pro fit2h,file=file,spline=spline

; Explanation : This procedure is for fitting height-time profiles using four functions:
;		exponential, exponential+linear, power law, power law+linear.
; Inputs      : file   - .sav file containing the following data:
;                          height     - the unit is km
;                          height_err - the unit is km
;                          time       - in the format of the UTC string
;                        when file is not given, an example file 'tmp.sav' will be used.
; Keywords    : spline - when the keyword 'spline' is set, the code computes a smooth cubic spline
;                        approximation to height, which, instead of real height, is used to derive
;                        the (non-fitted) velocity. It should be noted that the keyword 'spline'
;                        will call some IMSL functions. For IDL 8.2-8.7, IMSL is available on
;                        Windows and Linux, while unavailable on OS X; for IDL 8.0-8.1, IMSL is
;                        available on Linux, 32-bit Windows, and 32-bit OS X.
; Outputs     : Chi^2 derived from fitting;
;               onset time; velocity and height at onset;
;               images showing fitting results
; Written     : Version 1, X. Cheng, Dec. 01, 2018
; Modified    : Version 2, X. Cheng, Oct. 14, 2019
;               Version 3, C. Xing,  Oct. 22, 2019, adding keywords 'file', 'spline'; optimizing
;                                                   plotting


;--------------------------------------read data
if not keyword_set(file) then file='ht.sav'
restore,file,/ver
nstart=0
nend=n_elements(height)-5              ;;adjustable data range, here 5 is an arbitrary value

xtime=utc2tai(str2utc(anytim(time[nstart:nend],/vms)))
xtime=xtime-xtime[0]
time=time[nstart:nend]
height=height[nstart:nend]
height_err=height_err[nstart:nend]
xtime_err=fltarr(n_elements(height))   

;--------------------------------------fitting with function:a*exp(bt)+h0
x=xtime
y=height
y_err=height_err
expr='p[0]*exp(p[1]*x)+p[2]'
myfit1=mpfitexpr(expr,x,y,y_err,perror=sigma1,bestnorm=bestnorm,yfit=yfit1,dof=dof)
chisq1=bestnorm/dof

;--------------------------------------fitting with function:a*exp(b*(t-c))+dt+h0
expr='p[0]*exp(p[1]*(x-p[2]))+p[3]*x+p[4]'
myfit2=mpfitexpr(expr,x,y,y_err,perror=sigma2,bestnorm=bestnorm,yfit=yfit2,dof=dof)
chisq2=bestnorm/dof

;--------------------------------------fitting with function:a*t^b+h0
expr='p[0]*(x^p[1])+p[2]'
parinfo=replicate({value:0.D,fixed:0,limited:[0,0],limits:[0.D,0]},3)
parinfo.value=[0.,3.,10.]
parinfo.fixed=[0,0,0]
parinfo[1].limited=[1,1]
parinfo[1].limits=[0.,10.]             ;;adjustable constraints
myfit3=mpfitexpr(expr,x,y,y_err,perror=sigma3,bestnorm=bestnorm,yfit=yfit3,parinfo=parinfo,dof=dof)
chisq3=bestnorm/dof

;--------------------------------------fitting with function:a*(t-b)^c+dt+h0
expr='p[0]*((x-p[1])^p[2])+p[3]*x+p[4]'
parinfo=replicate({value:0.D,fixed:0,limited:[0,0],limits:[0.D,0]},5)
parinfo.value=[0.,0.,3.,10.,0.]
parinfo.fixed=[0,1,0,0,0]
parinfo[1].limited=[0,0]
parinfo[1].limits=[-1000.,1000.]
parinfo[2].limited=[1,1]
parinfo[2].limits=[0.,10.]
parinfo[3].limited=[1,1]
parinfo[3].limits=[0.,100.]
myfit4=mpfitexpr(expr,x,y,y_err,perror=sigma4,bestnorm=bestnorm,yfit=yfit4,parinfo=parinfo,dof=dof)
chisq4=bestnorm/dof

;--------------------------------------calculate velocity
if (keyword_set(spline)) then begin
 pp=imsl_cssmooth(xtime,height)
 y_spline=imsl_spvalue(xtime,pp)
 vel=deriv(xtime,y_spline)
 vel_err=derivsig(xtime,y_spline,xtime_err,height_err)
endif else begin
 vel=deriv(xtime,height)
 vel_err=derivsig(xtime,height,xtime_err,height_err)
endelse

vel1=myfit1[0]*myfit1[1]*exp(myfit1[1]*x)
vel2=myfit2[0]*myfit2[1]*exp(myfit2[1]*(x-myfit2[2]))+myfit2[3]
vel3=myfit3[0]*myfit3[1]*x^(myfit3[1]-1)
vel4=myfit4[0]*myfit4[2]*(x-myfit4[1])^(myfit4[2]-1)+myfit4[3]

;--------------------------------------plot
;pos11=[0.1,0.10,0.48,0.9]
;pos12=[0.6,0.10,0.98,0.9]

window,1,xs=1200,ys=600
loadct,0
!p.color=0
!p.background=255
!p.multi=[0,2,1,0,0]
loadct,12

averinter=float(xtime[-1])/float(n_elements(xtime)-1)
timemin=utc2str(tai2utc(utc2tai(str2utc(anytim(time[nstart],/vms)))-averinter))
timemax=utc2str(tai2utc(utc2tai(str2utc(anytim(time[nend],/vms)))+averinter))
timerange=[timemin,timemax]

hplerr=height+height_err
hmierr=height-height_err
hpl=[hplerr,yfit1,yfit2,yfit3,yfit4]
hmi=[hmierr,yfit1,yfit2,yfit3,yfit4]
yrange=[min(hmi)-(max(hpl)-min(hmi))/7.,max(hpl)+(max(hpl)-min(hmi))/7.]

utplot,time,height,timerange=timerange,yrange=yrange,psym=1,charsize=1.8,symsize=1.5,thick=3,xstyle=1, ystyle=1,ytitle='Height [km]',charthick=1
uterrplot,time,hplerr,hmierr,color=0,thick=3

outplot,time,yfit1,psym=-3,color=200,thick=3
outplot,time,yfit2,psym=-3,color=100,thick=3
outplot,time,yfit3,psym=-3,color=150,thick=3
outplot,time,yfit4,psym=-3,color=30,thick=3

onset=alog(myfit2[3]/myfit2[1]/myfit2[0])/myfit2[1]+myfit2[2]
onset_str=utc2str(tai2utc(utc2tai(str2utc(anytim(time[0],/vms)))+onset))
outplot,[onset_str,onset_str],yrange,psym=-3,color=100,thick=5,linestyle=1

al_legend,['ae!Ubt!N+h!D0!N	        '+cggreek('chi')+'!U2!N!X='+strmid(strtrim(string(chisq1),1),0,3)+' '$
,'ae!Ub(t-t!I0!U)!N+ct+h!D0!N	  '+cggreek('chi')+'!U2!N!X='+strmid(strtrim(string(chisq2),1),0,3)+' '$
,'at!Ub!N+h!D0!N	         '+cggreek('chi')+'!U2!N!X='+strmid(strtrim(string(chisq3),1),0,3)+' '$
,'a(t-t!D0!N)!Ub!N+ct+h!D0!N	 '+cggreek('chi')+'!U2!N!X='+strmid(strtrim(string(chisq4),1),0,3)+' ']$
,textcolor=[200,100,150,30],box=1,charsize=1.5,charthick=1,/left,spacing=2.

vplerr=vel+vel_err
vmierr=vel-vel_err
vpl=[vplerr,vel1,vel2,vel3,vel4]
vmi=[vmierr,vel1,vel2,vel3,vel4]
vrange=[min(vmi)-(max(vpl)-min(vmi))/7.,max(vpl)+(max(vpl)-min(vmi))/7.]

utplot,time,vel,timerange=timerange,yrange=vrange,psym=1,charsize=1.8,symsize=1.5,thick=3,ystyle=1,xstyle=1,ytitle='Velocity (km/s)',charthick=1,/noerase
uterrplot,time,vplerr,vmierr,color=0,thick=3

outplot,time,vel1,psym=-3,color=200,thick=3
outplot,time,vel2,psym=-3,color=100,thick=3
outplot,time,vel3,psym=-3,color=150,thick=3
outplot,time,vel4,psym=-3,color=30,thick=3

;--------------------------------------outputs
print,'-----------------------------------------------------------------------'
print,'Chi^2(exponential, exponential+linear, power law, power law+linear):'
print,chisq1,chisq2,chisq3,chisq4
print,'Exp Onset time (t>0):               ',onset_str
print,'Onset linear velocity (km/s):',2*myfit2[3],' @ '+onset_str
print,'Onset height (Mm):           ',(myfit2[0]*exp(myfit2[1]*(onset-myfit2[2]))+myfit2[3]*onset+myfit2[4])/1e3,' @ '+onset_str
print,'----------end ---------------------------------------------------------'

end