;WAVE_DIR: estimate the true direction of a height-time path, assuming uniform radial motion in 3D
;vr
;delta
;range='hi1','hi2','hi1,2'


;list of fitting functions
;function parameters: A = [delta, tau, t0] = [angle, Dsun/Vr, start time]
;the data are weighted more heavily if the density of data points is greater

pro tan_ang, x, a, f, pder
	f = (x-a[2])*cos(a[0])/(a[1]-(x-a[2])*sin(a[0]))
	 if (N_PARAMS() ge 4) then begin
	  pder = fltarr(N_ELEMENTS(x),3)
	  denom = (a[1]-(x-a[2])*sin(a[0]))^2
	  pder(*,0) =(x-a[2])*(x-a[2]-a[1]*sin(a[0]))/denom
	  pder(*,1) =-1*(x-a[2])*cos(a[0])/denom
	  pder(*,2) =-1*a[1]*cos(a[0])/denom
	 endif
end

;function parameters: p = [da,dv,d0,xa,xv,x0,sunda,sundv,sund0,t0] = [angle acc,angle vel, angle initial, radial distanc acc, radial distance vel, radial distance initial, Dsun acc, Dsun,time]
pro tan_ang_t, t, p, f, pder
	da = p[0]
	dv = p[1]
	d0 = p[2]
	xa = p[3]
	xv = p[4]
	x0 = p[5]
	sunda = p[6]
	sundv = p[7]
	sund0 = p[8]
	t0 = p[9]
        delta=.5*da*((t-t0)^2)+dv*(t-t0)+d0
        x=.5*xa*((t-t0)^2)+xv*(t-t0)+x0
        sund=.5*sunda*((t-t0)^2)+sundv*(t-t0)+sund0
	f = x*cos(delta)/(sund-(x*sin(delta)))
	 if (N_PARAMS() ge 4) then begin
	  df_dx=(1/x)*(f+(f^2)*(tan(delta)))
	  df_ddelta= (f^2) - (f*tan(delta))
	  df_dsund=-1*(f^2)/(x*cos(delta))
	  dx_dt0 = (-1*xa*(t-t0))-xv
	  ddelta_dt0 = (-1*da*(t-t0))-dv
	  dsund_dt0 = (-1*sunda*(t-t0))-sundv
	  pder = fltarr(N_ELEMENTS(x),10)
	  pder(*,0) = .5*df_ddelta*(t-t0)^2
	  pder(*,1) = df_ddelta*(t-t0)
	  pder(*,2) = df_ddelta
	  pder(*,3) = .5*df_dx*(t-t0)^2
	  pder(*,4) = df_dx*(t-t0)
	  pder(*,5) = df_dx
	  pder(*,6) = .5*df_dsund*(t-t0)^2
	  pder(*,7) = df_dsund*(t-t0)
	  pder(*,8) = df_dsund
	  pder(*,9) = (df_dx*dx_dt0) + (df_ddelta*ddelta_dt0) +(df_dsund*dsund_dt0)
	 endif
end

;function parameters: pr = [delta,v_s,t0,acc_s,tacc] = [angle,vr/Dsun, initial time at sun, acc/Dsun, start time acc]
pro tan_ang_r, t, pr, f, pder
;this includes a constant acc starting at the first data point collected
	delta = pr[0]
	v_s = pr[1]
	t0 = pr[2]
	acc_s = pr[3]
        tacc = pr[4]
       x_s=(.5*acc_s*((t-tacc)^2))+v_s*(t-t0)
	f = (x_s)*cos(delta)/(1-((x_s)*sin(delta)))
	 if (N_PARAMS() ge 4) then begin
          delta = pr[0]
	  v_s = pr[1]
	  t0 = pr[2]
	  acc_s = pr[3]
          x_s=(.5*acc_s*((t-tacc)^2))+v_s*(t-t0)
          f = (x_s)*cos(delta)/(1-((x_s)*sin(delta)))
	  df_dx_s=(1/x_s)*(f+(f^2)*(tan(delta)))
	  df_ddelta= (f^2) - (f*tan(delta))
	  dx_s_dt0 = (-1*acc_s*(t-tacc))-v_s
	  pder = fltarr(N_ELEMENTS(t),5)
	  pder(*,0) = df_ddelta
	  pder(*,1) = df_dx_s*(t-t0)
	  pder(*,2) = (df_dx_s*dx_s_dt0)
	  pder(*,3) = .5*df_dx_s*(t-tacc)^2
          pder(*,4) = -1*df_dx_s*acc_s*(t-tacc)
	 endif
end

;function parameters: pc = [delta,v_s,t0] = [angle,vr/Dsun, initial time at sun]
pro tan_ang_c, t, pc, f, pder
	delta = pc[0]
	v_s = pc[1]
	t0 = pc[2]
        x_s=v_s*(t-t0)
	f = (x_s)*cos(delta)/(1-((x_s)*sin(delta)))
	 if (N_PARAMS() ge 4) then begin
          delta = pc[0]
	  v_s = pc[1]
	  t0 = pc[2]
          x_s=v_s*(t-t0)
          f = (x_s)*cos(delta)/(1-((x_s)*sin(delta)))
	  df_dx_s=(1/x_s)*(f+(f^2)*(tan(delta)))
	  df_ddelta= (f^2) - (f*tan(delta))
	  dx_s_dt0 = -v_s
	  pder = fltarr(N_ELEMENTS(t),3)
	  pder(*,0) = df_ddelta
 	  pder(*,1) = df_dx_s*(t-t0)
	  pder(*,2) = (df_dx_s*dx_s_dt0)
	 endif
end


function wave_dir2a,htlist,vr=vr,delta=delta,t0=t0,range=range,sat=sat,func=func,hassat=hassat
;func default is tan_ang
;fit can only be tan_ang, tan_ang_r...(those listed above)
;don't load in only cor2 ht files.
;the errors and chisquareds are calculated with the assumption that the movies from which the jmaps were made have all available frames in them
;to calculate reasonable values for errors, multiple ht data sets must be collected for a given path on jmap.



IF(datatype(htlist) EQ 'UND') THEN htlist = dialog_pickfile(filter = '*.ht', path = '/sjmaps07/ht_data/', GET_PATH=mypath,/MULTIPLE_FILES)
 for i=0,n_elements(htlist)-1 DO BEGIN
   read_ht_file, htlist[i], tai, ht, pa

beginning:
 if (datatype(func) EQ 'UND') then func='tan_ang' 
if (datatype(range) EQ 'UND') then range='undefined' 

;getting time vs. tan(alpha) from data for different alpha ranges
	 if (range EQ 'undefined') then begin
	  x = double(tai)
	  y = double(tan(!pi*ht/180))
       	 endif
	 if (range EQ 'hi1') then begin  
	  cut_hi1=where(ht GT 3.28); -1 if only cor2, first element is 0 if hi1 or hi2
	  cut_hi2=where(ht GT 23.28);-1 if only hi1 cor2, first element minus 1 is first alpha in hi2
	  if (cut_hi1[0] EQ -1) then element_start=0 else element_start=cut_hi1[0]
	  if (cut_hi2[0] EQ -1) then element_end=n_elements(ht)-1 else element_end=cut_hi2[0]-1
          x = double(tai[element_start:element_end])
	  y = double(tan(!pi*ht[element_start:element_end]/180))
	 endif
	 if (range EQ 'hi2') then begin
	  cut_hi2=where(ht GT 18.36)
	  if (cut_hi2[0] EQ -1) then element_start=0 else element_start=cut_hi2[0]
	  element_end=n_elements(ht)-1
	  x = double(tai[element_start:element_end])
	  y = double(tan(!pi*ht[element_start:element_end]/180))
	 endif
	 if (range EQ 'hi1,2') then begin
	  cut_hi1=where(ht GT 3.28)
	  if (cut_hi1[0] EQ -1) then element_start=0 else element_start=cut_hi1[0]
	  element_end=n_elements(ht)-1
	  x = double(tai[element_start:element_end])
	  y = double(tan(!pi*ht[element_start:element_end]/180))
	 endif

  if (i EQ 0) then time=x
  if ~(i EQ 0) then time=[time,x]
  if (i EQ 0) then tan_alpha=y 
  if ~(i EQ 0) then tan_alpha=[tan_alpha,y]


 endfor


;determining range
alpha_min=min(atan(tan_alpha)*180/!pi)
alpha_max=max(atan(tan_alpha)*180/!pi)
if ((alpha_min GT 3.28) and (alpha_max LT 23.28)) then range='hi1'
if (alpha_min GT 18.36) then range='hi2'
if ((alpha_min GT 3.28) and (alpha_min LT 18.36) and (alpha_max GT 23.28)) then range='hi1,2'
if ((alpha_min LT 3.28) and (alpha_max LT 23.28) and (alpha_max GT 3.28)) then range='cor2,hi1'
if ((alpha_min LT 3.28) and (alpha_max GT 23.28)) then range='all3'
if ((alpha_min LT 3.28) and (alpha_max LT 3.28)) then range='cor2'


;sort data for curvefit? sort wrt time
tta=double(fltarr(n_elements(time),2)) 
tta[*,0]=time
tta[*,1]=tan_alpha
time=tta[sort(tta[*,0]),0]
tan_alpha=tta[sort(tta[*,0]),1]

;get total days
;min_yr=strmid(utc2str(tai2utc(time[0])),0,4)
;min_mo=strmid(utc2str(tai2utc(time[0])),5,2)
;min_dy=strmid(utc2str(tai2utc(time[0])),8,2)
;min_date= min_yr + min_mo + min_dy
;max_yr=strmid(utc2str(tai2utc(time[n_elements[time]-1])),0,4)
;max_mo=strmid(utc2str(tai2utc(time[n_elements[time]-1])),5,2)
;max_dy=strmid(utc2str(tai2utc(time[n_elements[time]-1])),8,2)
;max_date= max_yr + max_mo + max_dy
;total_days=ceil((time[n_elements(time)-1]-time[0])/3600/24)
;total_days=max_date-min_date+1
;if (strlen(total_days) GT 2) then
;total_mo= strmid(total_days,strlen(total_days)-4,

;getting dates for which data points have been taken
nti=n_elements(time)
for i=0,nti-1 do begin
;stop
  temp_yr=strmid(utc2str(tai2utc(time[i])),0,4)
  temp_mo=strmid(utc2str(tai2utc(time[i])),5,2)
  temp_dy=strmid(utc2str(tai2utc(time[i])),8,2)
  temp_date=temp_yr + temp_mo + temp_dy
  if (i EQ 0) then dates=temp_date
  if ~(i EQ 0) then begin
    if (dates[n_elements(dates)-1] EQ temp_date) then dates=dates else dates=[dates,temp_date]
  endif
endfor
;getting realtimes which correspond to times when induvidual frames were taken. Also Dsun from file
ut0=tai2utc(time[0])
ut1=tai2utc(time[nti-1])

for k=0,2 do begin
  if (k EQ 0) then begin
     tel='cor2'
     itype='img'
  endif
  if (k EQ 1) then begin
     tel='hi_1' 
     itype='img'
  endif
  if (k EQ 2) then begin
     tel='hi_2' 
     itype='img'
  endif
  if (datatype(sat) EQ 'UND') then sat='a'
    sat=STRLOWCASE(sat)
    s=scc_read_summary(date=[ut0,ut1], type=itype, spacecraft=sat, telescope=tel)
    if (n_elements(s) GT 1) then begin
;      if ((datatype(Dsun_start) EQ 'UND')) then begin
;         dsun_start_file=l[0]
;         secchi_prep,dsun_start_file,header,img
;         wcs=fitshead2wcs(header)
;         Dsun_start=wcs.position.dsun_obs
;Dsun_start=Dsun-.00005
;      endif
;      if ((i EQ n_elements(dates)-1) and (datatype(Dsun_end) EQ 'UND')) then begin
;        dsun_end_file=l[n_elements(1)-1]
;        secchi_prep,dsun_end_file,header,img
;        wcs=fitshead2wcs(header)
;        Dsun_end=wcs.position.dsun_obs   
;Dsun_end=Dsun+.0000005    
;      endif
      if (tel EQ 'cor2') then realtime_cor2=anytim2tai(s.date_obs)
      if (tel EQ 'hi_1') then realtime_hi1=anytim2tai(s.date_obs)
      if (tel EQ 'hi_2') then realtime_hi2=anytim2tai(s.date_obs)
    endif
endfor

help,realtime_cor2,realtime_hi1,realtime_hi2

Dsun_start = 1.5e+11
Dsun_end = 1.5e+11

;match time to realtime
deg=atan(tan_alpha)*180/!pi
;IF datatype(realtime_hi1) EQ 'UND' THEN message,'REALTIME_HI1 is undefined. Check mount points of FITS data.'
for i=0,n_elements(time)-1 do begin
  if (deg[i] LT 3.28) then begin
    difference=realtime_cor2-time[i]  
    temp_centertime=realtime_cor2[min(where(min(difference^2) EQ (difference^2)))]
    if (i EQ 0) then centertime=temp_centertime else centertime=[centertime,temp_centertime]
  endif
  if ((deg[i] GT 3.28) and (deg[i] LT 23.28) and (range NE 'hi2')) then begin
    difference=realtime_hi1-time[i]  
    temp_centertime=realtime_hi1[min(where(min(difference^2) EQ (difference^2)))]
    if (i EQ 0) then centertime=temp_centertime else centertime=[centertime,temp_centertime]
  endif
  if ((deg[i] GT 23.28) or ((deg[i] GT 18.36) and (range EQ 'hi2'))) then begin
    difference=realtime_hi2-time[i]  
    temp_centertime=realtime_hi2[min(where(min(difference^2) EQ (difference^2)))]
    if (i EQ 0) then centertime=temp_centertime else centertime=[centertime,temp_centertime]
  endif
endfor


;take averages of deg's for each centertime since there may be multiple data point corrsp to a given centertime
avg_deg=double(fltarr(n_elements(centertime)))
sample_variance_deg=double(fltarr(n_elements(centertime)))
n_datapoints=intarr(n_elements(centertime))
variance_avg_deg=double(fltarr(n_elements(centertime)))

for i=0,n_elements(centertime)-1 do begin
  sames=where(centertime[i] EQ centertime)
  for j=0,n_elements(sames)-1 do begin;calculating average deg from data points which corrsp. to same centertime
    deg_part=deg[sames[j]]
    if (j EQ 0) then sum_deg=deg_part else sum_deg=sum_deg + deg_part
  endfor  
  avg_deg[i]=sum_deg/n_elements(sames)
  n_datapoints[i]=n_elements(sames)
  if (n_datapoints[i] GT 5) then begin
    for j=0,n_datapoints[i]-1 do begin;calc sample variance from data points which corrsp. to same centertime
      sigma_deg_part=(avg_deg[i]-deg[sames[j]])^2
      if (j EQ 0) then sum_sigma_deg=sigma_deg_part else sum_sigma_deg=sum_sigma_deg + sigma_deg_part
    endfor
    sample_variance_deg[i]=sum_sigma_deg/(n_datapoints[i]-1)
  endif
  if ~(n_datapoints[i] GT 5) then sample_variance_deg[i]=0 
  variance_avg_deg[i]=sample_variance_deg[i]/n_datapoints[i]

endfor

;taking out degnerate data points from centertime,avg_deg, and sample_variance_deg
for i=0,n_elements(centertime)-1 do begin
  temp_newtime=centertime[i]
  temp_newdeg=avg_deg[i]
  temp_newvariance_deg=variance_avg_deg[i]
  temp_n_datapoints=n_datapoints[i]
  if (i EQ 0) then begin
    newtime=temp_newtime 
    newdeg=temp_newdeg 
    newvariance_deg=temp_newvariance_deg
    new_n_datapoints=temp_n_datapoints
  endif
  if ((i NE 0) and (newtime[n_elements(newtime)-1] NE temp_newtime)) then begin
    newtime=[newtime,temp_newtime]
    newvariance_deg=[newvariance_deg,temp_newvariance_deg]
    newdeg=[newdeg,temp_newdeg]
    new_n_datapoints=[new_n_datapoints,temp_n_datapoints]
  endif
endfor


;assuming maximum variance for data points with too few data points to properly estimate variance
for i=0, n_elements(newvariance_deg)-1 do begin
  max_variance=max(newvariance_deg)
  if (newvariance_deg[i] EQ 0) then newvariance_deg[i]=max_variance/new_n_datapoints[i]
endfor

;calculating variables used in fit
newtime= double(newtime)
newtan_alpha= double(tan(!pi*newdeg/180))
if (min(newvariance_deg) EQ 0) then begin
  sigmas=0
  for i=0,n_elements(newvariance_deg)-1 do begin
     temp_nvta=double(1.0/new_n_datapoints[i])
     if (i EQ 0) then newvariance_tan_alpha=temp_nvta else newvariance_tan_alpha=[newvariance_tan_alpha,temp_nvta]
  endfor
endif
if (min(newvariance_deg) NE 0) then begin
  sigmas=1
  newvariance_tan_alpha=double(newvariance_deg*(((cos(!pi*newdeg/180)^(-2))*!pi/180)^2))
endif

  
itmax=30000
tol=1e-14
weights=double(fltarr(n_elements(newtime)))
for i=0, n_elements(newtime)-1 do begin
  weights[i]=double(1./newvariance_tan_alpha[i]);using gaussian weighting
endfor  




trial_delta=[40,30,20,50,60,70,10,0,80,0,-10,-20,-30,-40,-50]
trial_vr=[250,275,300,325,350,375,400,425,450,475,500,225,200,175,150]
i=0
j=0
    

initial:

;INITIAL GUESSES AND CONVERSIONS
for j=0, n_elements(trial_delta)-1 do begin 
  for i=0, n_elements(trial_vr)-1 do begin 
    Dsun=(Dsun_start+Dsun_end)/2000;km
    sigma_Dsun=((abs(Dsun_start-Dsun_end)/2) + (5e+3/sqrt(2)))/1000;km
    if (datatype(refit) EQ 'UND') then begin
      Vr=trial_vr[i];km/s
      delta=trial_delta[j]
    endif
    ;read,'Vr=',Vr
    ;read, 'delta=', delta
    ;delta=;degrees
    velocity = (newdeg[1]-newdeg[0])/(newtime[1]-newtime[0]);angular
    delta=!pi*(delta)/180
    tau=Dsun/Vr
    t0=newtime[0]-(newdeg[0]/velocity)
    a = double([delta,tau,t0])
    vr_initial=vr
    delta_initial=delta*180/!pi
    t0_initial=utc2str(tai2utc(t0)) 
    p = double([0.,0.,delta,0.,vr,6.96e+5,0.,0.,Dsun,t0])
    pr = double([delta,vr/Dsun,t0,1e-3,min(newtime)]);pr = [delta,v_s,t0,acc_s,tacc]
    pc =double([delta,vr/Dsun,t0])

    if (func EQ 'tan_ang') then begin
      fit=CURVEFIT(newtime,newtan_alpha,weights,a,sigma,chisq=chisq,/double,function_name='tan_ang',iter=iter,itmax=itmax,status=status,tol=tol,yerror=yerror)
   
      delta=180*a[0]/!pi
      Vr=Dsun/a[1];km/s
      t0=utc2str(tai2utc(a[2]))
      sigma_delta=180*sigma[0]/!pi
      sigma_vr=vr*sqrt((sigma[1]/a[1])^2+(sigma_Dsun/Dsun)^2)
      sigma_t0=sigma[2]/3600
    endif

    if (func EQ 'tan_ang_r') then begin
      fit=CURVEFIT(newtime,newtan_alpha,weights,pr,sigma,chisq=chisq,/double,function_name='tan_ang_r',iter=iter,itmax=itmax,status=status,tol=tol,yerror=yerror)
      delta=180*pr[0]/!pi
      Vr0=Dsun*pr[1];km/s
      acc=pr[3]*Dsun
      tacc=utc2str(tai2utc(pr[4]))
      vr=vr0 + acc*(max(newtime)-pr[4])
      t0=utc2str(tai2utc(pr[2]))
      sigma_delta=180*sigma[0]/!pi
      sigma_vr=vr*sqrt((sigma[1]/pr[1])^2+(sigma_Dsun/Dsun)^2)
      sigma_t0=sigma[2]/3600

    endif

    if (func EQ 'tan_ang_c') then begin
      fit=CURVEFIT(newtime,newtan_alpha,weights,pc,sigma,chisq=chisq,/double,function_name='tan_ang_c',iter=iter,itmax=itmax,status=status,tol=tol,yerror=yerror)
      delta=180*pc[0]/!pi
      Vr=Dsun*pc[1];km/s
      t0=utc2str(tai2utc(pc[2]))
      sigma_delta=180*sigma[0]/!pi
      sigma_vr=vr*sqrt((sigma[1]/pc[1])^2+(sigma_Dsun/Dsun)^2)
      sigma_t0=sigma[2]/3600
    endif

    if (sigmas EQ 1) then begin 
      if ((chisq LT 1000)and (delta LT 180) and (delta GT -180) and (vr LT 700) and (vr GT 100)) then begin
         
        goto, final    
      endif
    endif
    if (sigmas EQ 0) then begin 
      if ((delta LT 180) and (delta GT -180) and (vr LT 700) and (vr GT 100)) then begin
        goto, final    
      endif
    endif

;    undefine, a
;    undefine, pr
;    undefine,sigma
;    undefine, chisq
;    undefine, iter
;    undefine, status
;    undefine, yerror
;    undefine, delta
;    undefine, vr
;    undefine,t0
;    undefine,fit
;    undefine,acc

  endfor
endfor

final: 
  

;print,htlist
;print, 'PA:',pa
;print, 'range:'+string(range)
;print, 'number of data points:'+string(n_elements(time)) 
;print, 'number of data points after averaging:'+string(n_elements(newtime))
;print, 'status:'+ string(status)
;print,'iter:'+string(iter)
;if (sigmas EQ 0) then print,'chisq(nosigmas):'+string(chisq)
;if (sigmas EQ 1) then print,'chisq:'+string(chisq)
;print,'vr_initial:'+string(vr_initial)
;print,'delta_initial:'+string(delta_initial)
;print,'t0_initial:'+string(t0_initial)
;if (sigmas EQ 0) then begin
;  print,'vr:'+string(vr)
;  print,'delta:'+string(delta)
;  print,'t0:'+string(t0)
;endif
;if (sigmas EQ 1) then begin
;  print,'vr:'+string(vr) + ' +- ' + string(sigma_vr)
;  print,'delta:'+string(delta) + ' +- ' + string(sigma_delta)
;  print,'t0:'+string(t0) + ' +- ' + string(sigma_t0) +'hrs'
;endif
;print, 'yerror(deg)'+string(atan(yerror)*180/!pi)
;if (func EQ 'tan_ang_r') then begin
; print, 'vr0'+string(vr0)
; print, 'acc'+string(acc)
; print, 'tacc'+string(tacc)
;endif

item99= 'Sat: ' + strtrim(sat)
item0='PA: '+ strtrim(string(fix(pa)),2) ; strmid(string(pa),6,strlen(string(pa)))
item1='range: ' + string(range)
item2='number of data points:'+strmid(string(n_elements(time)),6,strlen(string(n_elements(time)))) 
item3='status:' ;+ string(status)
item4='iter:' ;+string(iter)
if (sigmas EQ 0) then item5='chisq(nosigmas):'+strmid(string(chisq),6,strlen(string(chisq)))
if (sigmas EQ 1) then item5='chisq:'+strmid(string(chisq),6,strlen(string(chisq)))
item6='vr_initial:'+string(vr_initial)
item7='delta_initial:'+string(delta_initial)
item8='t0_initial:'+string(t0_initial)
;if (sigmas EQ 0) then begin
  item19='vr:'+strmid(string(vr),6,strlen(string(vr)))
  item20='delta:'+strmid(string(delta),6,strlen(string(delta)))
  item21='t0: '+ string(t0)
;endif
;if (sigmas EQ 1) then begin
  item9='Vr: '+ strtrim(string(fix(vr)),2) ;strmid(string(vr),6,strlen(string(vr))) + ' +- ' + strmid(string(sigma_vr),6,strlen(string(sigma_vr)))
  item10='delta: '+ strtrim(string(fix(delta)),2) ;strmid(string(delta),6,strlen(string(delta)))+' +- '+strmid(string(sigma_delta),6,strlen(string(sigma_delta)))
  item11='t0: '+string(t0) ; + ' +- ' + strmid(string(sigma_t0),6,strlen(string(sigma_t0))) + 'hrs'
;endif
item12='alpha(degrees)'
item14='days since t0'
item13=string(htlist)
item15='number of data points after averaging:'+string(n_elements(newtime))
item16='yerror(deg):'+strmid(string(atan(yerror)*180/!pi),5)
 
        t0 = anytim2tai(t0)
  	del = delta*!pi/180
        pa = pa*!pi/180
        rs=30*.696e6  
        tss=t0+rs/vr
        place=get_stereo_lonlat(string(utc2str(tai2utc(tss))),string(sat),/degrees)
        b0 = place[2]*!pi/180
	theta = acos(sin(del)*sin(b0) + cos(b0)*cos(del)*cos(pa)) * 180/!pi
	phi = acos((sin(del)*cos(b0) - cos(del)*sin(b0)*cos(pa)) / sin(!pi*theta/180)) * 180/!pi
        ts=utc2str(tai2utc(tss))

item30='theta:'+strmid(string(theta),6,strlen(string(theta))-10)
item31='phi:'+strmid(string(phi),6,strlen(string(phi))-10)
item32='ts: '+string(ts)













indp=newtime/3600/24-newtime[0]/3600/24
dp=180*atan(newtan_alpha)/!pi
dpfit=180*atan(fit)/!pi

window,0
plot,indp,dpfit,xminor=12,xtickinterval=1,ytickinterval=5,xtitle='days since first data point',ytitle=item12,charsize=2
oplot,indp,dp,psym=4
;items=[item0,item9,item10,item11,item1,item5,item16,item13]
IF (keyword_set(hassat)) THEN items = [item99,item0,item31,item30,item9,item10]$ 
	ELSE items=[item0,item31,item30,item9,item10] ;removed item11, item32 from end of list (times)
;**********times are incorrect, since they are hardwired for COR1, HI1, HI2 formatted jmaps, so removed from legend ;tjc 6-22-09
legend,string(items),/bottom,/right,spacing=1.2,box=0,charsize=2

return,a

refit=''
;read, 'want to refit?(y/n)  ', refit
if (refit EQ 'y') then begin
   vr_guess=''
   read, 'initial guess for vr(km/s):  ',vr_guess
   vr=float(vr_guess)
   delta_guess=''
   read, 'initial guess for delta(degrees):  ',delta_guess
   delta=float(delta_guess)
   ;print, initial guess for vr (type 'vr=...', then .c)
   ;stop
   ;print, initial guess for delta (type 'delta=...', then .c)
   ;stop
   goto, initial
endif

indp=newtime/3600/24-a[2]/3600/24
dp=180*atan(newtan_alpha)/!pi
dpfit=180*atan(fit)/!pi


window,0
plot,indp,dpfit,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12
oplot,indp,dp,psym=4
items=[item0,item9,item10,item11,item1,item5,item16,item13]
legend,string(items),/bottom,/right,spacing=1.2,box=0,charsize=1






;save the results
		
	   ans=''
;           read, 'want to save graph as is(y/n)? ', ans
           if (ans eq 'y') then begin
                set_plot,'ps'
                ans2=''
                read,'enter filename= ',ans2
                device,filename=ans2
                device,ysize=9,/in
		device,xsize=7,/in
                device, yoffset=0.75,/in				 
	        plot,indp,dpfit,title=ans2,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12
		oplot,indp,dp,psym=4
		legend,string(items),/bottom,/right,spacing=1.2,box=0,charsize=1
                device,/close
                select_windows ;set_plot,'x'
 ;print the results
                ans3=''
                read, 'want to print(y/n)? ', ans3
                if (ans3 eq 'y') then begin
                   spawn, 'lp -d hp4-136 '+ string(ans2)
	        endif
	       
	   endif 



;save the results w/o extra stuff

window,0
plot,indp,dpfit,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12
oplot,indp,dp,psym=4
items=[item0,item19,item20,item21,item1]
legend,string(items),/bottom,/right,spacing=1.2,box=0,charsize=1
		
	   ans=''
;           read, 'want to save graph w/o chisq and files(y/n)? ', ans
           if (ans eq 'y') then begin
                set_plot,'ps'
                ans2=''
                read,'enter filename= ',ans2
                device,filename=ans2
                device,ysize=9,/in
		device,xsize=7,/in
                device, yoffset=0.75,/in				 
	        plot,indp,dpfit,title=ans2,xminor=12,xtickinterval=1,ytickinterval=5,xtitle=item14,ytitle=item12
		oplot,indp,dp,psym=4
		legend,string(items_less),/bottom,/right,spacing=1.2,box=0,charsize=1
                device,/close
                select_windows ;set_plot,'x'
 ;print the results
                ans3=''
                read, 'want to print(y/n)? ', ans3
                if (ans3 eq 'y') then begin
                   spawn, 'lp -d hp4-136 '+ string(ans2)
	        endif
	       
	   endif 


wdelete
return, a


END