pro make_circle,sources,fluxes1,fluxes2,circle1,circle2,help=help,algorithm=algorithm,plot=plot

if keyword_set(help) then begin
	print,'PURPOSE:'
	print,'	Mirror selected Clean regions onto MEM maps.'
	print,'	Also calculates the fluxes of the regions.'
	print,'HELP:'
	print,'	pro make_circle,sources,fluxes1,fluxes2,circle1,circle2'
	print,'	inputs:	sources (number of sources,integer)'
	print,'		algorithm (default: mem, keyword must be specified if other)'
	print,'		i.e. by default, circles are drawn onto MEM maps, but user may'
	print,'		want circles to be drawn on Clean maps'
	print,'		plot (default:no, /plot if user wants circle regions to be '
	print,'			visualized)'
	print,'	outputs:fluxes1'
	print,'		fluxes2 (fluxes of selected regions,double arrays)'
	print,'		circle1 (perimeter of circle regions for display on maps)
	print,'		circle2 (    "     "    "       "     "     "    "    " )
	print,'AUTHOR:'
	print,'	Rick Pernak, Goddard Space Flight Center'
	return
endif

if keyword_set(algorithm) then algorithm=algorithm else algorithm='mem'

;restore all of the maps and flare info/parameters
restore,'/home/pernak/RHESSI/paper/figures/Clean/normal/idl_plots/clean_figs_legend.sav'
restore,'/home/pernak/RHESSI/paper/figures/MEM/mem_figs_legend.sav'
restore,'/home/pernak/RHESSI/paper/figures/Pixon/pixon_figs_legend.sav'
;restore,'/home/pernak/RHESSI/MEM/mem.sav'
;can be changed to any MEM or Clean data/maps

print,'Making circles on ',algorithm,' maps'
; open up spreadsheet text file, take centroids columns and produce arrays
if (algorithm eq 'pixon') then begin
	stat_file ='/home/pernak/RHESSI/paper/spreadsheet/pixon.csv'
	area_file = '/home/pernak/RHESSI/paper/flux_data/pixon_areas_1118.txt'
	dummy = pixon
	data = data_pix
endif else if (algorithm eq 'mem') then begin
	stat_file ='/home/pernak/RHESSI/paper/spreadsheet/mem.csv'
	area_file = '/home/pernak/RHESSI/paper/flux_data/clean_areas_1107.txt'
	dummy = mem
	data = data_mem
endif else if (algorithm eq 'clean') then begin
	stat_file ='/home/pernak/RHESSI/paper/spreadsheet/clean.csv'
	area_file = '/home/pernak/RHESSI/paper/flux_data/clean_areas_1107.txt'
	dummy = clean
	data = data_mem
endif else begin
	print,'Not a valid method'
endelse

openr,1,stat_file
list = strarr(sources)
readf,1,list
close,1

x1 = dblarr(sources)
y1 = dblarr(sources)
x2 = dblarr(sources)
y2 = dblarr(sources)

for ctr=0,(sources-1) do begin
	xoff = data(ctr).xyoff[0]
	yoff = data(ctr).xyoff[1]
	xyint = data(ctr).xyint
	xydim = size(dummy,/dimensions)
	xdim = xydim[0]	;size of array in the x direction (equals y dimension size)
	split = str_sep(list(ctr),9B)

	x1[ctr] = split[5]
	y1[ctr] = split[6]
	x2[ctr] = split[7]
	y2[ctr] = split[8]
	
	x1[ctr] = ((x1[ctr]-xoff)/xyint) + (xdim/2.)
	y1[ctr] = ((y1[ctr]-yoff)/xyint) + (xdim/2.)
	x2[ctr] = ((x2[ctr]-xoff)/xyint) + (xdim/2.)
	y2[ctr] = ((y2[ctr]-yoff)/xyint) + (xdim/2.)
endfor
	
;centroids are arrays of x and y values, 1:bright source, 2:faint source
centroids1 = dblarr(2,sources)
centroids2 = dblarr(2,sources)
centroids1[0,*] = x1
centroids1[1,*] = y1
centroids2[0,*] = x2
centroids2[1,*] = y2

; take area file and break it up into 2 areas (1st and 2nd circle)
;	areas are copied from flux_output.txt file (default gui file name) and pasted
; into own file for convenience

openr,1,area_file
areas = dblarr(sources*2)
readf,1,areas
close,1

a1 = dblarr(sources) ;area of 1st(brightest) source
a2 = dblarr(sources) ;area of 2nd source

;take even numbered areas, which correspond to the brightest source
for ctr=0,(sources*2-1),2 do begin
	a1[ctr/2] = areas[ctr]
endfor

;take odd numbered areas, which correspond to the weakest source
for ctr=1,(sources*2-1),2 do begin
	a2[(ctr-1)/2] = areas[ctr]
endfor

;make a smaller second circle
;this is done to minimize overlapping of flux circles
for ctr=0,(sources-1) do begin
	if (ctr eq 9) or (ctr eq 13) then begin
		a2[ctr] = a2[ctr]/4. ;this will make a second circle with half the initial (gui) radius
		if (ctr eq 12) then a2[ctr] = a2[ctr]*1.5
	endif
endfor

; make circles in map

fluxes1 = dblarr(sources)
fluxes2 = dblarr(sources)

circle1 = dblarr(360,360,sources)
circle2 = dblarr(360,360,sources)

for ctr=0,sources-1 do begin
	r1 = double(sqrt(a1[ctr]/!pi)/xyint)
	r2 = double(sqrt(a2[ctr]/!pi)/xyint) ; radius of the two circle regions in CLEAN flux
	;pick appropriate square region around the origin
	xr1 = [-ceil(r1),ceil(r1)]
	yr1 = [-ceil(r1),ceil(r1)]
	xr2 = [-ceil(r2),ceil(r2)]
	yr2 = [-ceil(r2),ceil(r2)]
	map = dummy[*,*,ctr] ; output from save file

	;produce circles for maps
	phi = findgen(360) * !dtor
	circlex1 = r1*cos(phi)+centroids1[0,ctr]
	circley1 = r1*sin(phi)+centroids1[1,ctr]
	circle1[0,*,ctr] = circlex1
	circle1[1,*,ctr] = circley1
	circlex2 = r2*cos(phi)+centroids2[0,ctr]
	circley2 = r2*sin(phi)+centroids2[1,ctr]
	circle2[0,*,ctr] = circlex2
	circle2[1,*,ctr] = circley2

	;finds points in a circle of calculated radius in the square region

	;brightest region:
	temp_x1 = [0]
	temp_y1 = [0]
	for x1 = xr1[0],xr1[1] do begin
		for y1 = yr1[0],yr1[1] do begin
			if ((x1^2 + y1^2) le (r1^2)) then begin
				if (temp_x1[0] eq 0) then begin
					temp_x1 = [x1]
					temp_y1 = [y1]
				endif else begin
					temp_x1 = [temp_x1,x1]
					temp_y1 = [temp_y1,y1]
				endelse
			endif
		endfor
	endfor

	;move the regions from the origin to where the sources are on the map
	new_x1 = temp_x1 + centroids1[0,ctr]
	new_y1 = temp_y1 + centroids1[1,ctr]

	;make flux arrays
	arr1 = map[new_x1,new_y1]
	if (ctr eq 0) then fluxes1 = [total(arr1)] else fluxes1 = [fluxes1,total(arr1)]

	;faintest region:
	temp_x2 = [0]
	temp_y2 = [0]
	for x2 = xr2[0],xr2[1] do begin
		for y2 = yr2[0],yr2[1] do begin
			if ((x2^2 + y2^2) le (r2^2)) then begin
				if (temp_x2[0] eq 0) then begin
					temp_x2 = [x2]
					temp_y2 = [y2]
				endif else begin
					temp_x2 = [temp_x2,x2]
					temp_y2 = [temp_y2,y2]
				endelse
			endif
		endfor
	endfor

	;move the regions from the origin to where the sources are on the map
	new_x2 = temp_x2 + centroids2[0,ctr]
	new_y2 = temp_y2 + centroids2[1,ctr]

	;make flux arrays
	arr2 = map[new_x2,new_y2]
	if (ctr eq 0) then fluxes2 = [total(arr2)] else fluxes2 = [fluxes2,total(arr2)]

;plotting/visualization, plots 25 maps with circle regions for user
if keyword_set(plot) then begin
		lev = [.1,.2,.3,.4,.5,.6,.7,.8,.9,1.]
		window,1
		contour,map,xrange=[0,128],yrange=[0,128],/xst,/yst,/iso,lev=lev*max(map)
		oplot,new_x1,new_y1,psym=3
		oplot,new_x2,new_y2,psym=3
		stop
endif

endfor
return
end