; $Id: graph.pro,v 1.40 2009/07/30 20:09:22 cooper Exp $
; $Log: graph.pro,v $
; Revision 1.40  2009/07/30 20:09:22  cooper
; put coordinates in carrmaps
;
; Revision 1.39  2009/07/30 15:44:56  cooper
; improved carrmap options
;
; Revision 1.38  2009/07/28 15:31:45  cooper
; fixed error in finding carrmap
;
; Revision 1.37  2009/07/27 20:03:54  cooper
; choose which euvi carrmap to see
;
; Revision 1.36  2009/07/24 19:55:01  cooper
; fixed carrmap window size
;
; Revision 1.35  2009/07/24 19:24:15  cooper
; explained input variables, added carrington map ability
;
; Revision 1.34  2008/07/30 22:08:56  sheeley
; added 2*pi to carr when it is negative
;
; Revision 1.33  2008/07/30 15:37:41  casto
; Changed title a bit
;
; Revision 1.32  2008/07/30 15:19:44  casto
; Added type to legend
;
; Revision 1.31  2008/07/30 15:12:09  casto
; Fixed title
;
; Revision 1.30  2008/07/30 14:58:26  casto
; fixed date bug
;
; Revision 1.29  2008/07/30 14:55:44  casto
; Added real date and type
;
; Revision 1.28  2008/07/29 19:59:09  sheeley
; reset car_rot for subsequent tries
;
; Revision 1.27  2008/07/29 15:59:16  casto
; graph.pro
;
; Revision 1.26  2008/07/22 22:10:24  sheeley
; fixed two bugs in the Carrington rotation edge effects procedure
;
; Revision 1.25  2008/07/22 19:39:29  sheeley
; adjusted box spacing
;
; Revision 1.24  2008/07/22 19:09:58  sheeley
; new CR notation
;
; Revision 1.23  2008/07/22 16:41:25  casto
; Added carrington rotation
;
; Revision 1.22  2008/07/22 16:16:55  sheeley
; chnaged < to lt and introduced phi0
;
; Revision 1.21  2008/07/22 14:56:12  casto
; Fixed carr long
;
; Revision 1.20  2008/07/21 19:09:46  casto
; Tweaked the legend
;
; Revision 1.19  2008/07/21 19:00:29  casto
; Added theta and phi
;
; Revision 1.18  2008/07/21 14:57:04  casto
; Fixed formatting
;
; Revision 1.17  2008/07/21 14:25:36  casto
; Changed acc size and added data points to r/t plot
;
; Revision 1.16  2008/07/18 19:43:45  sheeley
; changed (falpha) to (alpha) in one line
;
; Revision 1.15  2008/07/18 19:18:13  casto
; Fixed alpha bug
;
; Revision 1.14  2008/07/17 20:02:02  sheeley
; modified to display graph before saving it
;
; Revision 1.13  2008/07/17 15:34:24  casto
; Was in the middle of something
;
; Revision 1.12  2008/07/17 14:49:01  casto
; Don't have to exit display to print
;
; Revision 1.11  2008/07/17 13:03:16  casto
; Fixed trim namespace bug
;
; Revision 1.10  2008/07/16 21:52:00  sheeley
; added more legend 1 changes & moved print query to an earlier line
;
; Revision 1.9  2008/07/16 21:07:32  sheeley
; converted vf_i from R/hr to km/s
;
; Revision 1.8  2008/07/16 20:28:34  sheeley
; converted delta to degrees in first legend before putting it on the graph
;
; Revision 1.7  2008/07/16 20:09:56  sheeley
; added fast/slow option and changed the A/B parameter from (sp) to (side) which we have used in the past
;
; Revision 1.6  2008/07/16 19:11:08  casto
; Fixed small bugs/preferences
;
; Revision 1.5  2008/07/16 18:28:12  casto
; Calculates real distance to sun
;
; Revision 1.4  2008/07/15 20:21:32  nathan
; fixed ID tag
;
; Revision 1.3  2008/07/15 18:59:44  casto
; Made changes to interface and printing
;
; Revision 1.2  2008/07/15 18:01:14  casto
; Fixed what I was working on when I committed it
;
; Revision 1.1  2008/07/15 15:04:27  casto
; Beginning of curvefit
;

common share, type, radius

PRO fit, X, A, F, pder
common share, type, radius

curvef = obj_new(type, A)
F = curvef->r(X)
if n_params() GE 4 then begin
	pder(*,0) = curvef->dr0(X)
	pder(*,1) = curvef->dvi(X)
	pder(*,2) = curvef->dvf(X)
	pder(*,3) = curvef->dtau(X)
	pder(*,4) = -curvef->v(X)
ENDIF
end

PRO hif, X, A, F, pder
common share, type, radius

del = A[0]
r1 = A[1]
vf = A[2]

rho = (r1+vf*X)/radius
h = 1.-2.*rho*sin(del)+rho^2
F = acos((1.-rho*sin(del))/sqrt(h))

if n_params() ge 4 then begin
	pder(*,0) = rho*(rho-sin(del))/h
	pder(*,1) = cos(del)/(radius*h)
	pder(*,2) = X*cos(del)/(radius*h)
endif
end
pickfile = 1
while(pickfile) do begin
	filename = dialog_pickfile(filter='*.ht',path='/sjmaps08/ht_data')
	print, filename
	read_ht_file,filename,tai,ht,angle,date
	date = anytim2utc(date)
	if angle lt 180 then sat ="A" else sat = "B"
	angle *= !pi/180
        ;ans_kind = ''
        ;read,'Which kind of fit do you want? (fast(f) or slow(s)) ', ans_kind
        ;if (ans_kind eq 's') then begin
	    radius = (get_stereo_lonlat(date, sat))[0]/(695500.)
        ;endif else begin
            ;radius = 215.
        ;endelse
	b0 = (get_stereo_lonlat(date, sat, system='HCI'))[2]
	carr = (get_stereo_lonlat(date, sat, system='Carrington'))[1]
        if (carr lt 0.) then carr = carr+2.*!pi  ; nrs 7/30/08

	car_rot0 = get_stereo_carr_rot(date, sat) ; need an un-incremented value nrs 7/29/08
	print, 'carr = ',carr*180./!pi
	t = (tai-tai[0])/3600
	alpha = ht*!pi/180
	WINDOW, 0, XSIZE=600, YSIZE=400
	!p.multi = [0,1,1,0,0]
	plot, t, alpha*180./!pi, psym=2, charsize=2.5, $
	xtitle='t (hrs)', ytitle='!4a!3 (deg)', title=filename
	redo = ''
	read, 'Want to fit these data? (y/n) ', redo
	pickfile = (redo eq 'n')
endwhile

unsatisfied = 1
while(unsatisfied) do begin
	WINDOW, 0, XSIZE=700, YSIZE=900
	!p.multi = [0,1,2,0,0]
	plot, t, alpha*180./!pi, psym=2, charsize=2.5, $
	xtitle='t (hrs)', ytitle='!4a!3 (deg)', title=filename
	read, 't_boundary (beginning of upper fit in hrs) = ', t1
	hirange = where(t gt t1)
	hit = t[hirange]
	hialpha = alpha[hirange]
	read, 'delta (angle out of sky plane) = ', deltag
	deltag *= !pi/180.
	read, 'r1 (y-intercept of upper fit) = ', r1g
	read, 'vf = ', vfg
	vfg *= 1.2/232.
	params1 = double([deltag, r1g, vfg])
	yfit = CURVEFIT(hit, hialpha, 1.+0.*hit, params1, $
		FUNCTION_NAME='hif', /DOUBLE)
	delta = params1[0]
	r1 = params1[1]
	vf = params1[2]
	print, "delta = ", delta*180./!pi
	print, "r1 = ", r1
	print, "vf (upper fit) = ", vf*232./1.2
	oplot, hit, yfit*180./!pi
	r = radius*sin(alpha)/(cos(alpha-delta))
	plot, t, r, psym=2, charsize=3, xtitle='t (hrs)', ytitle='r/R_sun'
	oplot, t, r1+vf*t

	read, 't0 (time of max accel) = ', t0g
	read, 'tau (+/- timescale of accel) = ', taug
	read, 'r0 (radius at t0) = ', r0g
	read, 'vi = ', vig
	vig *= 1.2/232.
	params2 = double([r0g, vig, vfg, taug, t0g])
	type = ''
	read, "What type of function? ", type
	yfit2 = CURVEFIT(t, r, 1.+0.*t, params2,FUNCTION_NAME='fit',/DOUBLE)
	curvef = obj_new(type, params2)
	curvef->printout
	theta = acos(cos(b0)*cos(delta)*cos(angle)+sin(b0)*sin(delta))
	phi0 = atan(sin(angle)*cos(delta), $ 
	cos(b0)*sin(delta) - cos(delta)*cos(angle)*sin(b0))
        print,'phi0 = ',phi0*180./!pi
	phi = carr - phi0
	car_rot = car_rot0
	while (phi lt 0) do begin
		phi += 2*!pi
		car_rot += 1
	end
	car_rot -= fix(phi/(2*!pi))
	phi = phi mod (2*!pi) 

	plot, t, alpha*180./!pi, psym=2, charsize=2.5, $
	xtitle='t (hrs)', ytitle='!4a!3 (deg)',title=filename
	rho = curvef->r(t)/radius
	falpha = acos((1.-rho*sin(delta))/sqrt(1.-2.*rho*sin(delta)+rho^2))
	oplot, t, falpha*180./!pi
	plot, t, r, psym=2, charsize=3, xtitle='t (hrs)', ytitle='r/R_sun'
	oplot, t, r1+vf*t
	oplot, t, curvef->r(t)

	try = ''
	read, "Want to try again (y/n)? ", try
	unsatisfied = (try eq 'y')
endwhile

!p.multi = [0,2,2,0,0]
adjust = 1
while (adjust) do begin
	read, "What time should you cut it down to? ", cut
	WINDOW, 1, XSIZE=900, YSIZE=600
	ft = findgen(300)*cut/300.
	fr = curvef->r(ft)
	fv = curvef->v(ft)
	fa = curvef->a(ft)
	plot, ft, fv*232./1.2, charsize=1.5, $
		xtitle='t (hrs)', ytitle='v (km/s)'
	plot, ft, fa*232./(1.2*3.6), charsize=1.5, $
		xtitle='t (hrs)', ytitle= 'a (m/s^2)'
	plot, fr, fv*232./1.2, charsize=1.5, $ 
		xtitle='r/R_sun', ytitle='v (km/s)'
	plot, fr, fa*232./(1.2*3.6), charsize=1.5, $ 
		xtitle='r/R_sun', ytitle= 'a (m/s^2)'
	re = ''
	read, 'Want to adjust the time? (y/n) ', re
	adjust = (re eq 'y')
endwhile

set_plot, 'ps'
device, ysize=9.5,/in
device, yoffset=0.75,/in
device, /landscape

!p.position = [0.05,0.55,0.45,0.95]
plot, t, alpha*180./!pi, psym=2, charsize=1.5, $ 
	xtitle='t (hrs)', ytitle='!4a!3 (deg)'
rho = curvef->r(t)/radius
falpha = acos((1.-rho*sin(delta))/sqrt(1.-2.*rho*sin(delta)+rho^2))
oplot, t, falpha*180./!pi
 
arr = ['t1 = ' + strmid(strtrim(string(t1),1),0,4)+' hr', $
	"vf = " + strmid(strtrim(string(vf*232./1.2),1),0,5)+' km/s',$
	"r1 = " + strmid(strtrim(string(r1),1),0,5)+' R_s', $
	"!4d!3 = " + strmid(strtrim(string(delta*180./!pi),1),0,4), $
	'pa = ' + strmid(strtrim(string(angle*180./!pi),1),0,4), $
	'!4h!3 = ' + strmid(strtrim(string(theta*180./!pi),1),0,4), $
	'!4u!3 = ' + strmid(strtrim(string(phi*180./!pi),1),0,4), $
	'CR'+strmid(strtrim(string(car_rot),1),0,4)]

legend, arr, /top, /left, spacing=1.2, box=0, charsize=1.2

!p.position = [0.05,0.05,0.45,0.45]
plot, t, r, psym=2, charsize=1.5, xtitle='t (hrs)', ytitle='r/R_sun' 
oplot, t, curvef->r(t)

legend, curvef->toString(type),/top,/left,spacing=1.2,box=0,charsize=1.2
!p.position = [0.59,0.55,0.99,0.95]
plot, fr, fv*232./1.2, charsize=1.5, $ 
	xtitle='r/R_sun', ytitle='v (km/s)'

!p.position = [0.59,0.05,0.99,0.45]
plot, fr, fa*232./(1.2*3.6), charsize=1.5, $ 
	xtitle='r/R_sun', ytitle= 'a (m/s^2)'
datestr = '   '+strmid(utc2str(date),0,10)+' '+strmid(utc2str(date),11,5)
xyouts, 0.5, 1.0, filename+datestr+' UT', alignment=0.5, /normal
device, /close
set_plot, 'x'
spawn, 'display idl.ps' + '&'
!p.position = -1
!p.multi = -1

print, "phi (longitude): ", phi*180/!pi
print, "theta (latitude): ", theta*180./!pi
print, 'Carrington Rotation:: ',car_rot
viewcarr=''
read, 'Want to see the carrmap? (y/n): ', viewcarr
 
IF(viewcarr eq 'y') THEN BEGIN
	aorb=''
	notdone='y'
	cororeuvi=''
	criteria=''
	WHILE(notdone EQ 'y') DO BEGIN
	read, 'Cor2 or Euvi? (cor2/euvi): ', cororeuvi
	read, 'Which satellite? (a/b): ', aorb
	ew=''
	IF(cororeuvi EQ 'cor2') THEN BEGIN read, 'East or West limb? (E/W): ',ew
		ew=STRUPCASE(ew)
		read,'What solar radius? (3/5/7/10/13): ', criteria		
		CASE criteria OF
		'3':criteria=117  ;sid# in names Arnaud gave the files
		'5':criteria=113
		'7':criteria=121
		'10':criteria=125
		'13':criteria=129
		ELSE:print,'Not a valid solar radius'
		ENDCASE
		IF(datatype(criteria) NE 'STR') THEN BEGIN
			IF(aorb EQ 'b') THEN criteria++
			criteria='sid'+STRTRIM(STRING(criteria),2)
		ENDIF
	ENDIF ELSE IF(cororeuvi EQ 'euvi') THEN read, 'What wavelength? (171,195,284,304): ',criteria ELSE print,'Not a valid image'
	scar_rot=STRTRIM(STRING(FIX(car_rot)),2)
	map_img=findfile('/net/earth/secchi/carrmaps/'+aorb+'/'+cororeuvi+'/cm/png/*'+criteria+'*cr'+scar_rot+ew+'*')
	IF(map_img[0] EQ '') THEN print, 'Carrington Map for '+scar_rot+' not found: /net/earth/secchi/carrmaps/'+aorb+cororeuvi+'/cm/png/ ' $
	ELSE BEGIN map_img=map_img[n_elements(map_img)-1]
		map_img=read_png(map_img)
		WINDOW,2,xsize=880,ysize=500
		tvscl,map_img
		x=(phi*180./!pi)*2+79
		y=69-((theta*180./!pi)-180.)*2
		print,x,y
		xyouts,x-2,y-1,'(',/device
		xyouts,x+2,y-1,')',/device,color=1
		;;;Needs to be fixed, I used numbers rather than variables because I didn't know how to find values such as the lower left corner of the map
		 print,'Right Click to leave map reading mode'
		cursor,x,y,/device
		lat0='' & lon0=''
		WHILE(!MOUSE.button NE 4) DO BEGIN
			x=(x-79)/2. ;(80,70) is lower left corner of map in pixels with respect to the screen. Divide by two because
			y=ABS((y-69)/2.-180) ;map is (720x360 pixels). -180 to get 0 at north pole, instead of south
			lat='Latitude: '+STRTRIM(STRING(y),2)
			lon='Longitude: '+STRTRIM(STRING(x),2)
			IF(lat NE lat0 OR lon NE lon0) THEN BEGIN
				DEVICE,COPY=[600,450,250,50,500,10]
				xyouts,500,30,lon,/device,color=1,charsize=2
				xyouts,500,10,lat,/device,color=1,charsize=2
				lat0=lat & lon0=lon
			ENDIF
			cursor,x,y,/device
		ENDWHILE
	ENDELSE
	read,'Want to see another carrmap? (y/n): ', notdone
	ENDWHILE
ENDIF

tosave = ''
read, 'Want to save? (y/n) ', tosave
if (tosave eq 'y') then begin
	slashes = strsplit(filename, '/')
	name = dialog_pickfile(filter='*.ps', $
		path='/sjmaps08/ht_data/ps_files')
	spawn, 'mv idl.ps ' + name + '&'
endif else begin
	name='/sjmaps08/temp/printable.ps'
	spawn, 'mv idl.ps ' + name
	spawn, 'chmod a+w /sjmaps08/temp/printable.ps'
	spawn, 'rm -f idl.ps'
endelse

printer = ''
read, 'Want to print? (y/n) ', printer	
if (printer eq 'y') then spawn, 'lp -d hp4-136 '+name	
if (file_test('/sjmaps08/temp/printable.ps') NE -1) then spawn, 'rm -f /sjmaps08/temp/printable.ps'
end