pro make_sigmas,sigma_flux1,sigma_flux2,sigma_fwhm1,sigma_fwhm2

;March 6, 2006
;make a sigma array for the flux and fwhm arrays that Ed made on March 3, 2006

restore,'/home/pernak/RHESSI/forward_fit/data/dates_fwhm_flux.sav'

nbands = (size(fwhm1,/dim))[0] ;number of bands
nflares = (size(fwhm1,/dim))[1] ;number of flares

openr,1,'/home/pernak/RHESSI/forward_fit/data/filenames.txt'
saved = strarr(nflares)
readf,1,saved
close,1

;sigma arrays to be used in plots
sigma_flux1 = fltarr(nbands,nflares)
sigma_flux2 = fltarr(nbands,nflares)
sigma_fwhm1 = fltarr(nbands,nflares)
sigma_fwhm2 = fltarr(nbands,nflares)

;loop over flares
for ctr=0,nflares-1 do begin
	restore,saved[ctr]
	bands = size(sigma,/n_elements)/2
	
	;assign zeroes to all bands of "bad flares"
	if (ctr eq 8) or (ctr eq 12) then begin
		sigma_flux1(*,ctr) = 0.
		sigma_flux2(*,ctr) = 0.
		sigma_fwhm1(*,ctr) = 0.
		sigma_fwhm2(*,ctr) = 0.
	endif else begin

		;loop over bands for first source
		for ctr2=0,bands-1 do begin
			sigma_flux1(ctr2,ctr) = sigma[ctr2*2].srcflux
			sigma_fwhm1(ctr2,ctr) = sigma[ctr2*2].srcfwhm
		endfor

		;loop over bands for second source
		for ctr3=0,bands-1 do begin
			sigma_flux2(ctr3,ctr) = sigma[(ctr3*2)+1].srcflux
			sigma_fwhm2(ctr3,ctr) = sigma[(ctr3*2)+1].srcfwhm
		endfor		
	endelse

endfor

;sf = sqrt((sigma_flux1/flux1)^2+(sigma_flux2/flux2)^2)

;make x axis (fwhm ratio) between .5 and 1
;flux ratio is kept the same because indices between fwhm and flux should match
;for ctr1 = 0,nbands-1 do begin ;loop over band
;	for ctr2 = 0,nflares-1 do begin ;loop over flare
;		if (fwhm1[ctr1,ctr2] gt fwhm2[ctr1,ctr2]) then $
;			size[ctr1,ctr2]=fwhm1[ctr1,ctr2]/(fwhm1[ctr1,ctr2]+fwhm2[ctr1,ctr2]) $
;			else size[ctr1,ctr2]=fwhm2[ctr1,ctr2]/(fwhm1[ctr1,ctr2]+fwhm2[ctr1,ctr2])
;		fluxes[ctr1,ctr2] = flux1[ctr1,ctr2]/(flux1[ctr1,ctr2]+flux2[ctr1,ctr2])
;	endfor
;endfor

;check for finite values
dummy1=fwhm1/fwhm2
dummy2=flux1/flux2

;weed out the bad sources *will be one dimensional array now*
good = where(finite(dummy1) and finite(dummy2) and $
	(dummy1 gt 0) and (dummy2 gt 0),ngood); and (sf lt 4.),ngood)

;calculate rows and columns
row = good/nbands
column = good mod nbands

;regenerate flux arrays (with finite, non-zero values) and make ratio arrays
new_flux1 = dblarr(ngood)
new_flux2 = dblarr(ngood)
fluxes = dblarr(ngood)

for ctr=0,ngood-1 do begin
	new_flux1[ctr] = flux1[column[ctr],row[ctr]]
	new_flux2[ctr] = flux2[column[ctr],row[ctr]]
	if (new_flux1[ctr] gt new_flux2[ctr]) then $
		fluxes[ctr] = new_flux1[ctr]/(new_flux1[ctr]+new_flux2[ctr]) else $
		fluxes[ctr] = new_flux2[ctr]/(new_flux1[ctr]+new_flux2[ctr])
endfor

fwhm1 = fwhm1[good]
fwhm2 = fwhm2[good]

set_plot,'ps'
device,filename='asymmetry_size_flux.ps'

;ploterr,size,fluxes,sf,psym=4,/xst,/yst,xr=[0,2],yr=[0,2],
;window,0
plot,fwhm1,fluxes,psym=4,/xst,/yst,xtitle='FWHM1',$
	ytitle='Flux1/(Flux1 + Flux2), Flux1 > Flux2'
for ctr=0,ngood-1 do begin
	coordinates=string(row[ctr],column[ctr],format='("(",1i,",",1i,")")')
	xyouts,fwhm1[ctr],fluxes[ctr],strcompress(coordinates,/rem)
endfor
legend,['(Flare Number, Band Number)'],box=0,/bottom,/right

;plot least squares line
x = findgen(30)/10
coeff1 = poly_fit(fwhm1,fluxes,1,/double)
oplot,coeff1[0]+x*coeff1[1]

;window,1
plot,fwhm2,fluxes,psym=4,xtitle='FWHM2',$
	ytitle='Flux1/(Flux1 + Flux2), Flux1 > Flux2'
for ctr=0,ngood-1 do begin
	coordinates=string(row[ctr],column[ctr],format='("(",1i,",",1i,")")')
	xyouts,fwhm2[ctr],fluxes[ctr],strcompress(coordinates,/rem)
endfor
legend,['(Flare Number, Band Number)'],box=0,/bottom,/right

;plot least squares line
x = findgen(30)/10
coeff2 = poly_fit(fwhm2,fluxes,1,/double)
oplot,coeff2[0]+x*coeff2[1]

device,/close
set_plot,'x'

;window,2
;plot,fwhm2,fluxes,psym=4,xtitle='FWHM',$
;	ytitle='Flux1/(Flux1 + Flux2), Flux1 > Flux2'
;oplot,fwhm1,fluxes,psym=2
;oplot,coeff1[0]+x*coeff1[1]
;oplot,coeff2[0]+x*coeff2[1],lines=2

return
end