PRO instvel2,inar,pgno,ima,name,vxylog

; A program to obtain instantaneous velocities and trajectories
; from the object-tracking result array.
; Version 2 - trajectory smoothed by SMSPLINE,
;	      classification and more data in vxylog

;INPUTS: inar - 3D array, result of object tracking (trim)
;        pgno - array, list of subscripts of objects
;	 ima  - underlying image--one of the series
;	 name - string containing a part of filename

;OUTPUTS: vxylog - 2-D array containing for each accepted object
;	0. the identification number nn (referred to inar),
;	1. lifetime in minutes (derived from number of frames),
;	2. direct travel distance in arcsec,
;	3. distance from umbral center in arcsec,
;	4. average velocity magnitude (km/s),
;	5. classification of speed shape:
;	   -2: hill, -1: decrease, 0: const., 1: increase, 2: valley.
;		This array is saved in file vxy[name].log
;
;	  save-files  vxy[name].nn , one per accepted object.
;	        Each file contains arrays xs,ys (smoothed trajectory)
;		and v (instantaneous velocity magnitude).

;Parameters
	reb=5		;  rebin factor for visualization
	sca=0.125	;  arcsec per pixel
	lag=44.5265	;  seconds
	x0=123		;  x-position of umbral center (pix)
	y0=54		;  y-position of umbral center (pix)

;==================================================================

scale=sca*725.		;  km per pixel
si=size(ima)
sx=si(1)		;  parameter sx determined automatically
sip=size(pgno)
vxylog=fltarr(6,sip(1))		; define output array

loadct,3
window,0,xsiz=512,ysiz=382
window,1,xsiz=512
window,2,xsiz=si(1)*reb,ysiz=si(2)*reb
tvscl,rebin(ima,si(1)*reb,si(2)*reb)

for i=0, sip(1)-1 do begin	; BEGIN LOOP OVER OBJECTS

print,'Hit "s" to STOP, any other key to continue'
if get_kbrd(1) eq 's' then goto,fin

pos=inar(1,*,pgno(i))	;  positions in 1-D form
pos=pos(where(pos,c))	;  c = number of frames with object
print,' '
print,'Seq. number:',i,'    Object index:',pgno(i)
print,'No. of frames =',c,'     Lifetime (mnts.) =',(c-1)*lag/60.
x=pos mod sx
y=pos/sx

;Travel distance

dst=sqrt(float((x(c-1)-x(0))^2+(y(c-1)-y(0))^2))*sca
print,'Direct travel distance =',dst,'"',dst*725,' km'

; Average distance from umbral center

avx=round(total(x)/c)	; average position
avy=round(total(y)/c)
cdst=sqrt((avx-x0)^2+(avy-y0)^2)*sca
print,'Distance from umbral center =',cdst,'"'

; Average speed

  t=dindgen(c)*lag		;time coordinate in seconds
  sut=total(t)			;linear least squares fit...
  sut2=total(t*t)
  sux=total(x)
  suy=total(y)
  sutx=total(t*x)
  suty=total(t*y)
  ssx=sutx-sut*sux/c
  ssy=suty-sut*suy/c
  sst=sut2-sut*sut/c
  vx=scale*ssx/sst			;slopes (velocity km/s)
  vy=scale*ssy/sst
  avvm=sqrt(vx*vx+vy*vy)		;avg.velocity magnitude
  print,'Average speed =',avvm,' km/s'

; Output stuff

fil='vxy'+name+'.'+strtrim(pgno(i),2)
vxylog(0,i)=pgno(i)		; object number
vxylog(1,i)=(c-1)*lag/60.	; lifetime (min.)
vxylog(2,i)=dst			; travel distance
vxylog(3,i)=cdst		; distance from UC center
vxylog(4,i)=avvm		; average speed

;Interactive smoothing of trajectory and calculation of speed

;Initial parameters for SMSPLINE
tens=5			;tension (decreases with value)
nods=4  		;number of nodes

changepar:

wset,1
!p.multi=[0,0,2,0,0]
plot,x,ytitle='x-trajectory',/yst
xs=SMSPLINE(x,tens,nods)
oplot,xs,col=150
plot,y,ytitle='y-trajectory',/yst
ys=SMSPLINE(y,tens,nods)
oplot,ys,col=150
!p.multi=0

;Speed
vx=(xs-shift(xs,1))*scale/lag
vy=(ys-shift(ys,1))*scale/lag
vx=vx(1:c-1)
vy=vy(1:c-1)
v=sqrt(vx*vx+vy*vy)

wset,0
plot,v,ytitle='velocity magnitude (km/s)'

print,'Smoothing parameter =',nods
read,'New smoothing parameter (ge 3) or class (-2,-1,0,1,2): ',nods
; set new nods or store the class and continue
if nods ge 3 then goto,changepar else vxylog(5,i)=nods

;Trajectory plot

wset,2
plots,x*reb,y*reb,/device,col=255,thick=1
plots,xs*reb,ys*reb,/device,col=255,thick=2

print,'Discard this object ? (y/n)'
if get_kbrd(1) eq 'y' then vxylog(1,i)=0 else save,xs,ys,v,file=fil

plots,x*reb,y*reb,/device,col=0,thick=1
plots,xs*reb,ys*reb,/device,col=0,thick=2

endfor			; END LOOP OVER OBJECTS

fin:
w=where(vxylog(1,*),m)
print,'End of program. Objects selected:',m
if m ne 0 then begin
  vxylog=vxylog(*,w)
  save,vxylog,file='vxy'+name+'.log'
endif
wdelete,1
wset,0
end