pro slice2vel,file=file,fx=fx,fy=fy,min=min,max=max,save=save,spline=spline

; Explanation : This procedure is for displaying the slice-time plot and doing height-time
;               measurement for the top of the erupting feature, and then calculating the velocity
;               as the function of time.
; Inputs      : file    - .sav file containing the following data:
;                           2D slice array - each 1D slice lies along x direction, +x direction is
;                                            the radial direction; +y direction is the time
;                                            increasing direction.
;                           time           - corresponding to the slice array; in the format of the
;                                            UTC string.
;                         when file is not given, an example file 'slice.sav' will be used.
; Keywords    : fx,fy   - the factors of adjusting window size; when they are not set, the default
;                         is (1,1).
;               min,max - set the min and max of slice image intensity; when they are not set, the
;                         default is (0.2,2.5).
;               save    - when the keyword 'save' is set, the code will save time, height, velocity
;                         of selected points in a .sav file.
;               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 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     : png file displaying the slice-time plot and velocity-time plot;
;               .sav file containing time [xpoint_ut_str], height [height,height_err (km)], velocity
;               [vel,vel_err (km/s)] of selected points.
; Written     : Version 1, X. Cheng, Dec. 01, 2016
; Modified    : Version 2, X. Cheng, Mar. 17, 2017, adding keywords 'min', 'max'
;               Version 3, C. Xing,  Oct. 22, 2019, adding keywords 'spline'; optimizing plotting;
;                                                   using right mouse button to end the cursor


;--------------------------------------read data
if not keyword_set(file) then file='slice.sav'
restore,file,/ver
slice=rotate(slice,4)

;--------------------------------------adjust size of window
if not keyword_set(fx) then fx=1
if not keyword_set(fy) then fy=1

slice_tmp=congrid(slice,n_elements(slice[*,0])*fx,n_elements(slice[0,*])*fy)
ss=size(slice_tmp)
print,'X,Y dimension:',ss[1],ss[2]

;--------------------------------------plot slice
loadct,0
!p.background=255
!p.color=0
pos1=[0.06,0.13,0.46,0.93]
pos2=[0.56,0.13,0.96,0.93]
window,1,xs=ss[1]*2,ys=ss[2]

if not keyword_set(min) then min=0.2
if not keyword_set(max) then max=2.5

image_plot=bytscl(alog10(slice_tmp),min=min,max=max)
tvimage,image_plot,position=pos1

base_y0=0.0
unite=6.955e2/975.0                    ;Mm/arc
base_y1=n_elements(slice[0,*])*0.6*unite
utplot,time,[base_y0,base_y1],/nodata,position=pos1,xstyle=1,ystyle=1,/noerase,ytitle='Height (Mm)',charsize=2,charthick=1

;--------------------------------------cursor
print,'------------------------------------------------------------------------'
print,'Press the left mouse button to select points'
print,'Press the right mouse button to select the last point and end the cursor'
print,'------------------------------------------------------------------------'

loadct,34
cursor,x_s,y_s,/down
xpoint=x_s
ypoint=y_s
plots,x_s,y_s,psym=2
print,x_s,y_s

repeat begin
 h:
 cursor,x_s,y_s,/down
 if min(abs(xpoint-x_s)) eq 0 then begin
  goto,h
 endif else begin
  xpoint=[xpoint,x_s]
  ypoint=[ypoint,y_s]
  plots,x_s,y_s,psym=2
  print,x_s,y_s
 endelse
endrep until !Mouse.Button eq 4

loadct,13
tmp=sort(xpoint)
xpoint=xpoint[tmp]
ypoint=ypoint[tmp]
plots,xpoint,ypoint,color=65,symsize=2

;--------------------------------------calculate and plot velocity
height=ypoint*1e3
height_err=height
height_err[*]=unite*0.6*5*1e3          ;uncertainty is 5-6 pixel.
xpoint_err=fltarr(n_elements(xpoint))

if (keyword_set(spline)) then begin
 pp=imsl_cssmooth(xpoint,height)
 ypoint_spline=imsl_spvalue(xpoint,pp)
 vel=deriv(xpoint,ypoint_spline)
 vel_err=derivsig(xpoint,ypoint_spline,xpoint_err,height_err)
endif else begin
 vel=deriv(xpoint,height)
 vel_err=derivsig(xpoint,height,xpoint_err,height_err)
endelse

loadct,0
xpoint_ut=str2utc(anytim(xpoint+anytim(time[0]),/vms))
utplot,xpoint_ut,vel,ytitle='Velocity (km/s)',psym=-2,position=pos2,/noerase,charsize=2,charthick=1
eutplot,xpoint_ut,vel,vel_err,/unc

;--------------------------------------save
if (keyword_set(save)) then begin
 filestr=strmid(time[0],0,11)
 xpoint_ut_str=utc2str(xpoint_ut)
 save,filename=filestr+'.sav',slice,time,base_y0,base_y1,xpoint,ypoint,xpoint_ut_str,height,height_err,vel,vel_err
endif

device,decomposed=0
imager=tvrd(true=1)
write_png,filestr+'.png',imager,r,g,b

end