pro sep_ang,dates,coordinates,help=help

if keyword_set(help) then begin
	print,'PURPOSE:'
	print,'	Convert centroid arrays, which are in units of arcsec, to latitude '
	print,'	and longitude, then finds separation angles and orientation angles '
	print,'	(alpha) for flares in "Analysis of RHESSI Flares Using Radio '
	print,'	Astronomical Techniques".'
	print,'HELP:'
	print,'pro sep_ang,dates,coordinates'
	print,'	inputs: none'
	print,'	outputs:coordinates structure, which contains: '
	print,'		s*: separation angle in radians (dbl array, for all 3 methods)'
	print,'		a*: orientation angle (alpha) in radians (dblarr for 3 methods)'
	print,'		dates: dates of flares (string array)'
	print,'		new_s: new values for old s, converted to Mm (dblarr 3 methods)'
	print,'AUTHOR:'
	print,'	Rick Pernak, Goddard Space Flight Center'
return
endif

openr,1,'/home/pernak/RHESSI/paper/spreadsheet/mem.csv'
mem_arr = strarr(25)
readf,1,mem_arr
close,1

openr,1,'/home/pernak/RHESSI/paper/spreadsheet/clean.csv'
clean_arr = strarr(25)
readf,1,clean_arr
close,1

openr,1,'/home/pernak/RHESSI/paper/spreadsheet/pixon.csv'
pixon_arr = strarr(25)
readf,1,pixon_arr
close,1

;dates of flares only used in paper, can be changed
dates=['2002-02-20','2002-07-23','2002-08-20','2002-09-08','2002-09-20',$
'2002-11-09','2003-03-18','2003-06-02','2003-06-09','2003-06-10',$
'2003-10-19','2003-10-24','2003-11-19','2004-10-30','2004-10-31',$
'2004-10-31','2004-11-06','2004-12-01','2005-01-15','2005-01-20',$
'2005-07-13','2005-07-30','2005-08-22','2005-09-10','2005-09-17']

centroid_plot,25,bm,fm,bc,fc,bp,fp
spawn,'rm centroid*.ps'

mem_s = dblarr(25)
clean_s = dblarr(25)
pixon_s = dblarr(25)

;convert old separation s to Mm
;MEM old s
for ctr=0,24 do begin
	split = str_sep(mem_arr(ctr),'	')
	old_s = double(split[9])
	;call program for determining solar radius in arcsec on a given day
	pbr = pb0r(date,/arcsec)
	r = pbr(2)
	mem_s[ctr] = old_s/r*700 ;Sun is 700 Mm in radius
endfor

;Clean old s
for ctr=0,24 do begin
	split = str_sep(clean_arr(ctr),'	')
	old_s = double(split[9])
	;call program for determining solar radius in arcsec on a given day
	pbr = pb0r(date,/arcsec)
	r = pbr(2)
	clean_s[ctr] = old_s/r*700 ;Sun is 700 Mm in radius
endfor

;Pixon old s
for ctr=0,24 do begin
	split = str_sep(pixon_arr(ctr),'	')
	old_s = double(split[9])
	;call program for determining solar radius in arcsec on a given day
	pbr = pb0r(date,/arcsec)
	r = pbr(2)
	pixon_s[ctr] = old_s/r*700 ;Sun is 700 Mm in radius
endfor

new_s = dblarr(3,25)
new_s[0,*] = mem_s
new_s[1,*] = clean_s
new_s[2,*] = pixon_s

;convert to latitude and longitude
ll_bm = xy2lonlat(bm,dates)
ll_fm = xy2lonlat(fm,dates)
ll_bc = xy2lonlat(bc,dates)
ll_fc = xy2lonlat(fc,dates)
ll_bp = xy2lonlat(bp,dates)
ll_fp = xy2lonlat(fp,dates)

;assign phi
phi_bm = ll_bm[0,*]*!dtor
phi_bc = ll_bc[0,*]*!dtor
phi_bp = ll_bp[0,*]*!dtor
phi_fm = ll_fm[0,*]*!dtor
phi_fc = ll_fc[0,*]*!dtor
phi_fp = ll_fp[0,*]*!dtor

;assign co-latitude: 90 deg - latitude
theta_bm = (90-ll_bm[1,*])*!dtor
theta_bc = (90-ll_bc[1,*])*!dtor
theta_bp = (90-ll_bp[1,*])*!dtor
theta_fm = (90-ll_fm[1,*])*!dtor
theta_fc = (90-ll_fc[1,*])*!dtor
theta_fp = (90-ll_fp[1,*])*!dtor

sep_mem = dblarr(25)
sep_clean = dblarr(25)
sep_pixon = dblarr(25)

;calculate coordinates in MEM
x_bm = sin(theta_bm)*cos(phi_bm)
y_bm = sin(theta_bm)*sin(phi_bm)
z_bm = cos(theta_bm)
x_fm = sin(theta_fm)*cos(phi_fm)
y_fm = sin(theta_fm)*sin(phi_fm)
z_fm = cos(theta_fm)

;calculate separation angle in MEM using dot product
for ctr=0,24 do begin
	sep_mem[ctr]=asin(sqrt(1-((x_bm[ctr]*x_fm[ctr])+(y_bm[ctr]*y_fm[ctr])+(z_bm[ctr]*z_fm[ctr]))^2))
endfor

;calculate coordinates in Clean
x_bc = sin(theta_bc)*cos(phi_bc)
y_bc = sin(theta_bc)*sin(phi_bc)
z_bc = cos(theta_bc)
x_fc = sin(theta_fc)*cos(phi_fc)
y_fc = sin(theta_fc)*sin(phi_fc)
z_fc = cos(theta_fc)

;calculate separation angle in Clean using dot product
for ctr=0,24 do begin
	sep_clean[ctr]=asin(sqrt(1-((x_bc[ctr]*x_fc[ctr])+(y_bc[ctr]*y_fc[ctr])+(z_bc[ctr]*z_fc[ctr]))^2))
endfor

;calculate coordinates in Pixon
x_bp = sin(theta_bp)*cos(phi_bp)
y_bp = sin(theta_bp)*sin(phi_bp)
z_bp = cos(theta_bp)
x_fp = sin(theta_fp)*cos(phi_fp)
y_fp = sin(theta_fp)*sin(phi_fp)
z_fp = cos(theta_fp)

;calculate separation angle in Pixon using dot product
for ctr=0,24 do begin
	sep_pixon[ctr]=asin(sqrt(1-((x_bp[ctr]*x_fp[ctr])+(y_bp[ctr]*y_fp[ctr])+(z_bp[ctr]*z_fp[ctr]))^2))
endfor

;calculate length of R2-R1, the connecting vector between two components
r2r1_mem=sqrt(2*(1-sin(theta_bm)*sin(theta_fm)*cos(phi_bm-phi_fm)-cos(theta_bm)*cos(theta_fm)))
r2r1_clean=sqrt(2*(1-sin(theta_bc)*sin(theta_fc)*cos(phi_bc-phi_fc)-cos(theta_bc)*cos(theta_fc)))
r2r1_pixon=sqrt(2*(1-sin(theta_bp)*sin(theta_fp)*cos(phi_bp-phi_fp)-cos(theta_bp)*cos(theta_fp)))

;calculate alpha
alpha_mem = asin(sqrt(1-((sin(theta_bm)+sin(theta_fm))*sin((phi_fm-phi_bm)/2)/r2r1_mem)^2))
alpha_clean = asin(sqrt(1-((sin(theta_bc)+sin(theta_fc))*sin((phi_fc-phi_bc)/2)/r2r1_clean)^2))
alpha_pixon = asin(sqrt(1-((sin(theta_bp)+sin(theta_fp))*sin((phi_fp-phi_bp)/2)/r2r1_pixon)^2))

coordinates = $
	{smem:sep_mem,sclean:sep_clean,spixon:sep_pixon,$
	amem:alpha_mem,aclean:alpha_clean,apixon:alpha_pixon,$
	new_s:new_s}

;{llbm:ll_bm,llfm:ll_fm,llbc:ll_bc,llfc:ll_fc,llbp:ll_bp,llfp:ll_fp

$alias rm rm
$rm centroid*.ps
$alias rm 'rm -i'

return
end