;  $Id: fitall.pro,v 1.2 2024/10/08 19:35:48 nathan Exp $
;  $Log: fitall.pro,v $
;  Revision 1.2  2024/10/08 19:35:48  nathan
;  contents of SATPLOT_2.3.tar.gz
;
;  Revision 2.3 2024/09/09 00:00:00 penteado
;  Updating satplot package to V2.3
;+
;
; Name:       fitall
; 
; Purpose:    fit elongation vs. time to infer constant direction, speed, launch and arrival times (for ICMEs, CIRs)
;   	      using the Fixed-Phi-Fitting method, assuming a point-like shape
;             (e.g. Sheeley et al. 1999 ApJ, Rouillard et al. 2008 GRL, Moestl et al. 2009,2010,2011)
; 	      to be used within the SolarSoft SATPLOT package	
; 
; Parameters:  file_in... the name of the satplot .ht file which contains the elongation time data
;
; Keywords:  	cut...  two integers indicating how many data points of the CME 
;   	    	    	are not used at the beginning and end of the tracking data;    	    	
;    	    	    	default=[0,0] 
;               dir initial direction, defaults to +60 (for STEREO B), -60 (for A)
;               _extra anyo other keyword are passed on to utplot
;    	    	     
; Calling sequence: results=fitfp_satplot(filename, cut=[cutfront, cutback])       
;
; Example: results=fitfp_satplot('satplot__20110213_120000__20110218_120000__pa270_d10_B.ht')       	
;  
; Side effects: uses spice to get Earth and STEREO positions
;               uses fitj.pro 
;                  
; History:    12 Feb 2012: started rewriting to be used with satplot
;             7 Mar 2012: continued, now works with .ht files from satplot
;             12 Mar 2012: streamlined the input process
;             09 Sep 2024: updated documentation (penteado)
;
; Author:     Christian Moestl, SSL Berkeley and University of Graz, Austria
;             Tanja Rollett - University of Graz, Austria
;
;-

function fitall, file_in, cut=cut, dir=dir, _extra = extra_keywords

;---------------------------------------
;read in .ht files produced by satplot
str = STRARR(200)

OPENR, 10, file_in
   dummy = ''
   i = 0L
   WHILE NOT(EOF(10)) DO BEGIN
      READF, 10, dummy
      str[i] = dummy
      i = i + 1
   ENDWHILE
CLOSE, 10
str = str[0:i-1]
searchsc=strmid(str,13,1);
spacecraft= searchsc(3);
a = where( strmid(str,0,1) eq '#' )
a = max(a)+1
str = str[a:*]

track_date=strmid(str,9,19)
track_y=float(strmid(str,30,5))
inst=strmid(str,35,4)

;final structure containing the times, elongations, instrument, spacecraft
tr = {track_date: track_date, track_y: track_y, inst: inst, sc:spacecraft}
;---------------------------------------



radtodeg=180./!dpi;

time=tr.track_date
epsilon=tr.track_y  
scname=tr.sc;

; set parameters for the imaging observatory
if scname eq 'A' then begin 
    othersc='B'
    angle_signum=-1;
     IF KEYWORD_SET(dir) THEN BEGIN
      initialphi=dir; 
     ENDIF ELSE BEGIN 
      initialphi=-60;
     ENDELSE        
   end
    
if scname eq 'B' then begin
    othersc='A'
    angle_signum=1
    IF KEYWORD_SET(dir) THEN BEGIN
      initialphi=dir; 
     ENDIF ELSE BEGIN 
      initialphi=+60;
     ENDELSE  
  end
;;;___________________________________________-

;use correct sign for A (-1) and B (+1)  
epsilon=epsilon*angle_signum  

s=size(tr.track_date)

IF KEYWORD_SET(cut) THEN BEGIN    
    time2=strarr(s(1)-cut(0)-cut(1))    
    epsilon2=fltarr(s(1)-cut(0)-cut(1))
    time2=time(0+cut(0):s(1)-cut(1)-1)
    epsilon2=epsilon(0+cut(0):s(1)-cut(1)-1)
    time=0
    epsilon=0
    ;    print, time2, epsilon2, size(time2)   
    time=time2
    epsilon=epsilon2
    print, 'Data points cutted from track for fitting:', cut
    print, '   '
ENDIF    





;get STEREO positions
pos1=get_stereo_lonlat(tr.track_date(s(1)/2), scname, system='HEE')
pos2=get_stereo_lonlat(tr.track_date(s(1)/2), scname, system='HEEQ')

posother=get_stereo_lonlat(tr.track_date(s(1)/2), othersc, system='HEE')
posother2=get_stereo_lonlat(tr.track_date(s(1)/2), othersc, system='HEEQ')

;get position of the Earth
posearth=get_stereo_lonlat(tr.track_date(s(1)/2), 'earth', system='HEE')


print, '   '
print, 'IMAGING OBSERVATORY:'
print, 'Position of STEREO-', scname, ' in HEEQ long/lat:', pos1(1)*radtodeg, pos1(2)*radtodeg
print,'Position of STEREO-', scname, ' in HEE long/lat:',  pos2(1)*radtodeg, pos2(2)*radtodeg

print, 'on ', tr.track_date(s(1)/2) , ' (average date of the elongation values)'
print, '  '



print, 'IN SITU OBSERVATORY:'
print, 'Position of STEREO-', othersc, ' in HEEQ long/lat:',  posother(1)*radtodeg, posother(2)*radtodeg
print, 'Position of STEREO-', othersc, ' in HEE long/lat:',  posother2(1)*radtodeg, posother2(2)*radtodeg
print, 'on ', tr.track_date(s(1)/2)
print, '   '

tes=anytim(time)-anytim(time(0)); time to seconds for fitting

;convert time back from seconds to UTC
tzero=anytim(time(0)) 




;;;;;;;;;;;;;;;;;;;;;;;;;; SLOPE FITTING 

degtorad=!dpi/180;
H0=pos1(0);Imaging observatory to Sun in km
;starting points
phi=initialphi*degtorad; 




vr=500;
tinit=tes(0)-86400; first HI1 observation minus one day


X=[phi, vr, tinit];
common myfit,xueber,yueber, H0ueber
xueber=tes;
yueber=epsilon;
H0ueber=H0;
RES=0;
;Fixed Phi fitting
RES=AMOEBA(1.0e-5,FUNCTION_NAME='fitj',FUNCTION_VALUE=values,P0=X,scale=[1,100,1e4]);


X=[phi, vr, tinit];
common myfit2, xueber2, yueber2, dst, scnameueber
xueber2=tes;
yueber2=epsilon;
scnameueber=scname;
dst=H0;

;Harmonic Mean fitting
RESH=0;
RESH=AMOEBA(1.0e-5,FUNCTION_NAME='fith2',FUNCTION_VALUE=values,P0=X,scale=[1,100,1e4]);


print, 'STARTING POINT ', X(0)/degtorad, X(1), X(2)
print, 'Fixed- Phi Minimization: ', RES[0]/degtorad,'  ', RES[1],  '   ',RES[2]
print, 'Harmonic Mean Minimization: ', RESH[0]/degtorad,'  ', RESH[1],  '   ',RESH[2]







;**********************************************************
;------------- FP DIRECTIONS
earthangle=RES[0]*radtodeg+pos1(1)*radtodeg
otherangle=(pos1(1)*radtodeg-posother(1)*radtodeg )+RES[0]*radtodeg 
print, '          '
print, 'ICME direction: positive longitude angles -> solar west (right)'
print, 'FP:'
print, 'ICME Angle to Earth   ', num2str(earthangle, FORMAT='(F10.2)')
print, 'ICME Angle to STEREO-',othersc,'   ', num2str(otherangle,FORMAT='(F10.2)')
print, '  '

;------------- HM DIRECTIONS
earthangleh=RESH[0]*radtodeg+pos1(1)*radtodeg
otherangleh=(pos1(1)*radtodeg-posother(1)*radtodeg )+RESH[0]*radtodeg 
print, 'HM:'
print, 'ICME Angle to Earth   ', num2str(earthangleh, FORMAT='(F10.2)')
print, 'ICME Angle to STEREO-',othersc,'  ', num2str(otherangleh,FORMAT='(F10.2)')




;------FP Arrival time calculation
;launch time
tl0=RES(2)

;Arrival time: constant speed Sun to 1 AU, launch on Sun center
arrivaltime_st_sec=tzero+tl0+posother(0)/RES[1]; in seconds
arrst=strmid(anytim(arrivaltime_st_sec, /vms),0,17); as date

arrivaltime_earth_sec=tzero+tl0+posearth(0)/RES[1]
arrearth=strmid(anytim(arrivaltime_earth_sec, /vms),0,17)

arrivaltime_ace_sec=arrivaltime_earth_sec-1.5*1e6/RES[1];
arrace=strmid(anytim(arrivaltime_ace_sec, /vms),0,17)


;------HM Arrival time calculation
;launch time
tl0h=RESH(2)

;check if the spacecraft is hit at all: -90, +90 around direction
if cos(otherangleh/radtodeg) > 0 then begin 
 arrivaltime_st_sech=tzero+tl0h+posother(0)/(RESH[1]*cos(otherangleh/radtodeg)); in seconds
 arrsth=strmid(anytim(arrivaltime_st_sech, /vms),0,17); as date
endif else begin
 print, 'STEREO-',scname,' is not hit by the HM circle'
endelse

if cos(earthangleh/radtodeg) > 0 then begin
 arrivaltime_earth_sech=tzero+tl0h+posearth(0)/(RESH[1]*cos(earthangleh/radtodeg))
 arrearthh=strmid(anytim(arrivaltime_earth_sech, /vms),0,17)
 arrivaltime_ace_sech=arrivaltime_earth_sech-1.5*1e6/(RESH[1]*cos(earthangleh/radtodeg));
 arraceh=strmid(anytim(arrivaltime_ace_sech, /vms),0,17)
endif else begin
 print, 'Earth is not hit by the HM circle'
endelse


print, '         '

print, 'RESULTS SUMMARY'
print, 'FP:'
print, 'phi =   ', num2str(RES[0]*radtodeg,FORMAT='(F10.2)')
print, 'speed =   ', num2str(RES[1],FORMAT='(F11.2)')
launchtime=strmid(anytim(tzero+tl0, /vms),0,17)
print, 'launch time solar center:  ', launchtime
print, 'arrivaltime ACE: ', arrace
print, 'arrivaltime Earth: ', arrearth
print, 'arrivaltime STEREO-',othersc,': ', arrst


elements=size(yueber)
normalized_residue=fitj(RES)/elements(1)

print, 'fitting residue FP:', normalized_residue
print, '  '


print, 'HM:'
print, 'phi =   ', num2str(RESH[0]*radtodeg,FORMAT='(F10.2)')
print, 'speed =   ', num2str(RESH[1],FORMAT='(F11.2)')
launchtimeh=strmid(anytim(tzero+tl0h, /vms),0,17)
print, 'launch time solar center:  ', launchtimeh
;check if the Earth is hit by the HM circle
if cos(earthangleh/radtodeg) > 0 then begin
 print, 'arrivaltime ACE: ', arraceh
 print, 'arrivaltime Earth: ', arrearthh
 endif else begin
  arrearthh='0';
  arraceh='0';
endelse

;check if the other STEREO is hit by the HM circle
if cos(otherangleh/radtodeg) > 0 then begin
  print, 'arrivaltime STEREO-',othersc,': ', arrsth
  endif else begin
    arrsth='0';
  endelse
 
 
normalized_residueh=fith2(RESH)/elements(1)
print, 'fitting residue HM:', normalized_residueh
print, '  '



;FP: synthetic elongation time function from fits

;make interpolated time array so the curve is smooth
timeinterp=findgen(101)
tesinterp=findgen(101)

timestep=(anytim(time(elements(1)-1))-anytim(time(0)))/100;
 for i=0,100 do begin
  timeinterp(i)=anytim(time(0))+i*timestep;
  tesinterp(i)=i*timestep;  
 endfor

delta=90*degtorad-RES[0]     
;rho=RES[1]*(tes-RES[2])/H0;
rho=RES[1]*(tesinterp-RES[2])/H0;
fitmin=atan(rho*cos(delta)/(1-rho*sin(delta) ) )/degtorad


;HM: synthetic elongation time function from fits
b=RESH[0]
;because of the problems with the acosine one needs to switch the beta back to negative 
;so that its positive in fith.pro
if scname eq 'A' then b=-RESH[0]
v=RESH[1]
t=(tesinterp-RESH[2])
a=(2*dst)/(v*t)-cos(b)
c=sin(b)
fitminh=-acos( (-c^2+a*sqrt((c^2)*(-1+a^2+c^2)))/(c*(a^2+c^2))      ) /degtorad  

;--------------------------------------------------------------------------
;draw figure in window

set_plot,'X'

window, xsize=800,ysize=1000, retain=2
!p.multi=[0,1,2]
!p.background = long(16777215)
!p.color = 0



;********upper panel fixed phi
 ymax=max(tr.track_y)+10
 
 ;plotting time - elongation
 utplot, time, abs(epsilon), psym=4, ytitle='Elongation in degree', charsize=1.5,  $
 background=255,color=0,title=['STEREO-'+scname+' HI1/2 ICME FP track fitting'],yrange=[0,round(ymax)], $
 _extra = extra_keywords;
 ;plotting the fit
 outplot, anytim(timeinterp,/vms), abs(fitmin), color=0, thick=2, psym=0;, xrange=[xmin, xmax],yrange=[0,ymax]
 
 xyouts,tes(2),ymax-4,['CME angle to Earth = '$
 +num2str(earthangle,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2

 
xyouts,tes(2),ymax-8,['to STEREO-'+scname+' = '$
 +num2str(RES[0]/degtorad,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2
 
xyouts,tes(2),ymax-12,['to STEREO-'+othersc+' = '$
 +num2str(otherangle,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2

 
xyouts,tes(2),ymax-17,['V = '+num2str(RES[1],FORMAT='(I14)')+ $
' km/s'], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-23,['launch time: '+launchtime], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-26,['Arrival ACE: '+arrace], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-29,['Arrival ST'+othersc+': '+arrst], /DATA, Color=0, charsize=1.2
xyouts,tesinterp(n_elements(tesinterp)-30),ymax-3,['Fitting Residue: '+num2str(normalized_residue, format='(F5.3)')], /DATA, Color=0, charsize=1.2

;plot HI1 borders
hi12border=[18.7, 18.7, 18.7]
hi1border=[4, 4, 4]
hi1borderup=[24, 24, 24]
outplot, time(0:3), hi12border, linestyle=2, color=0
outplot, time(0:3), hi1border, linestyle=1, color=0
outplot, time(0:3), hi1borderup, linestyle=1, color=0



;********lower panel harmonic mean




 ;plotting time - elongation
 utplot, time, abs(epsilon), psym=4, ytitle='Elongation in degree', charsize=1.5,  $
 background=255,color=0,title=['STEREO-'+scname+' HI1/2 ICME HM track fitting'],yrange=[0,round(ymax)], $
 _extra = extra_keywords;
 ;plotting the fit
 outplot, anytim(timeinterp,/vms), abs(fitminh), color=0, thick=2, psym=0;, xrange=[xmin, xmax],yrange=[0,ymax]
 

 xyouts,tes(2),ymax-4,['CME angle to Earth = '$
 +num2str(earthangleh,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2

 
xyouts,tes(2),ymax-8,['to STEREO-'+scname+' = '$
 +num2str(RESH[0]/degtorad,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2
 
xyouts,tes(2),ymax-12,['to STEREO-'+othersc+' = '$
 +num2str(otherangleh,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2

 
xyouts,tes(2),ymax-17,['V = '+num2str(RESH[1],FORMAT='(I14)')+ $
' km/s'], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-23,['launch time: '+launchtimeh], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-26,['Arrival ACE: '+arraceh], /DATA, Color=0, charsize=1.2

xyouts,tes(2),ymax-29,['Arrival ST'+othersc+': '+arrsth], /DATA, Color=0, charsize=1.2
xyouts,tesinterp(n_elements(tesinterp)-30),ymax-3,['Fitting Residue: '+num2str(normalized_residueh, format='(F5.3)')], /DATA, Color=0, charsize=1.2

;plot HI1 borders
hi12border=[18.7, 18.7, 18.7]
hi1border=[4, 4, 4]
hi1borderup=[24, 24, 24]
outplot, time(0:3), hi12border, linestyle=2, color=0
outplot, time(0:3), hi1border, linestyle=1, color=0
outplot, time(0:3), hi1borderup, linestyle=1, color=0









;make jpeg output
filejpg='fpfit_'+file_in+'.jpg'
;x2jpeg, filejpg
 write_png, 'fpfit_'+file_in+'.png', tvrd()























;make similar eps output


set_plot,'PS'
 thickness=3
 device, /encapsulated, $
 filename='fpfit_'+file_in+'.eps',$
 xsize=22, ysize=27, /color, bits_per_pixel=8 

 ; the plot elongation vs. time
 utplot, time, abs(epsilon), psym=4, ytitle='Elongation in degree', charsize=1.5,  $
    background=255,color=0,title=['STEREO-'+scname+' HI1/2 ICME FP track fitting'],yrange=[0,round(ymax)], $
    charthick=thickness, thick=thickness;, xrange=[launchtime, time(n_elements(time)-1)]
 ;the plot for the fit 
 outplot, anytim(timeinterp,/vms), abs(fitmin), color=0, thick=thickness, psym=0;, xrange=[xmin, xmax],yrange=[0,ymax]
 
 xyouts,tes(2),ymax-4,['CME angle to Earth = '$
  +num2str(earthangle,FORMAT='(F10.1)')+' deg'], $
  /DATA, Color=0, charsize=1.2, charthick=thickness
 
 xyouts,tes(2),ymax-8,['to STEREO-'+scname+' = '$
  +num2str(RES[0]/degtorad,FORMAT='(F10.1)')+' deg'], $
  /DATA, Color=0, charsize=1.2, charthick=thickness
 
 xyouts,tes(2),ymax-12,['to STEREO-'+othersc+' = '$
  +num2str(otherangle,FORMAT='(F10.1)')+' deg'], $
  /DATA, Color=0, charsize=1.2, charthick=thickness
 
 xyouts,tes(2),ymax-17,['V = '+num2str(RES[1],FORMAT='(I14)')+ $
 ' km/s'], /DATA, Color=0, charsize=1.2, charthick=thickness

 xyouts,tes(2),ymax-23,['launch time: '+launchtime], /DATA, Color=0, charsize=1.2, charthick=thickness
 xyouts,tes(2),ymax-26,['Arrival ACE: '+arrace], /DATA, Color=0, charsize=1.2, charthick=thickness
 xyouts,tes(2),ymax-29,['Arrival ST'+othersc+': '+arrst], /DATA, Color=0, charsize=1.2, charthick=thickness
 xyouts,tesinterp(n_elements(tesinterp)-30),ymax-2,['Fitting Residue: '+num2str(normalized_residue, format='(F5.3)')], /DATA, Color=0, charsize=1.2, charthick=thickness

;indicate borders of HI1/HI2 field of view
hi12border=[ 18.7, 18.7, 18.7]
hi1border=[4, 4, 4]
hi1borderup=[24, 24, 24]
outplot, time(0:3), hi12border, linestyle=2, color=0, thick=thickness
outplot, time(0:3), hi1border, linestyle=1, color=0, thick=thickness
outplot, time(0:3), hi1borderup, linestyle=1, color=0, thick=thickness




;********lower panel harmonic mean




 ;plotting time - elongation
 utplot, time, abs(epsilon), psym=4, ytitle='Elongation in degree', charsize=1.5,  $
 background=255,color=0,title=['STEREO-'+scname+' HI1/2 ICME HM track fitting'],yrange=[0,round(ymax)], $
 _extra = extra_keywords, charthick=thickness, thick=thickness
 ;plotting the fit
 outplot, anytim(timeinterp,/vms), abs(fitminh), color=0, thick=thickness, psym=0;, xrange=[xmin, xmax],yrange=[0,ymax]
 

 xyouts,tes(2),ymax-4,['CME angle to Earth = '$
 +num2str(earthangleh,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2, charthick=thickness

 
xyouts,tes(2),ymax-8,['to STEREO-'+scname+' = '$
 +num2str(RESH[0]/degtorad,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2, charthick=thickness
 
xyouts,tes(2),ymax-12,['to STEREO-'+othersc+' = '$
 +num2str(otherangleh,FORMAT='(F10.1)')+' deg'], $
 /DATA, Color=0, charsize=1.2,charthick=thickness

 
xyouts,tes(2),ymax-17,['V = '+num2str(RESH[1],FORMAT='(I14)')+ $
' km/s'], /DATA, Color=0, charsize=1.2, charthick=thickness

xyouts,tes(2),ymax-23,['launch time: '+launchtimeh], /DATA, Color=0, charsize=1.2, charthick=thickness

xyouts,tes(2),ymax-26,['Arrival ACE: '+arraceh], /DATA, Color=0, charsize=1.2,charthick=thickness

xyouts,tes(2),ymax-29,['Arrival ST'+othersc+': '+arrsth], /DATA, Color=0, charsize=1.2, charthick=thickness
xyouts,tesinterp(n_elements(tesinterp)-30),ymax-2,['Fitting Residue: '+num2str(normalized_residueh, $
	format='(F5.3)')], /DATA, Color=0, charsize=1.2, charthick=thickness

;plot HI1 borders
hi12border=[18.7, 18.7, 18.7]
hi1border=[4, 4, 4]
hi1borderup=[24, 24, 24]
outplot, time(0:3), hi12border, linestyle=2, color=0, thick=thickness
outplot, time(0:3), hi1border, linestyle=1, color=0, thick=thickness
outplot, time(0:3), hi1borderup, linestyle=1, color=0, thick=thickness


















device, /close      















;make a structure with the results
results={launchtime_fp:launchtime, arrival_earth_fp:arrearth, $
	arrival_ace_fp:arrace, arrival_stereo_fp:arrst, $
	angle_earth_fp: earthangle, angle_observer_fp: RES(0)*radtodeg, $
         angle_other_stereo_fp:otherangle,speed_fp:RES(1), residue_fp:normalized_residue,$
        launchtime_hm:launchtimeh, arrival_earth_hm:arrearthh, $
	arrival_ace_hm:arraceh, arrival_stereo_hm:arrsth, $
	angle_earth_hm: earthangleh, angle_observer_hm: RESH(0)*radtodeg, $
         angle_other_stereo_hm:otherangleh,speed_hm:RESH(1), residue_hm:normalized_residueh }

 
  
 
set_plot,'X'

;return the structure
return, results


end