PRO instvel,inar,ima,name,vxylog,new=new

; A program to obtain instantaneous velocities and trajectories
; from the object-tracking result array (type TRIM).
; Version 5 - basic, universal. All objects in INAR are analyzed.
;	      trajectory smoothed by SMSPLINE,
;	      classification and general data in vxylog
;	      INW and OUT in a common logfile
;	      semi-automatical classification of speed shapes
;		(function CLASSIF.PRO)
;	    29 March 2001
;
;INPUTS: inar - 3D array, result of object tracking (trim)
;	 ima  - underlying image--one of the series
;	 name - string containing a part of filename
;
;KEYWORD NEW  - If set, new output array and log file is created.
;		If not, old log file is restored and the work
;		continues with next objects.	
;
;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 (class):
;	   -2: hill, -1: decrease, 0: const., 1: increase, 2: valley.
;	6. INW/OUT ident.:  INW: 1, OUT: 0
;	7. time-averaged intensity
;	  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).

on_error,1

;Parameters =======================================================

	reb=5		;  rebin factor for visualization
	sca=0.125	;  arcsec per pixel
	lag=21.99	;  seconds
	x0=40		;  x-position of umbral center (pix)
	y0=42		;  y-position of umbral center (pix)
	iph=10000.	;  mean quiet photosphere intensity

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

filn='vxyud/vxy'+name+'.'

scale=sca*725.		;  km per pixel
si=size(ima)
sx=si(1)		;  parameter sx determined automatically
siz=size(inar)
if siz(0) ne 3 then message,'Input data array must be 3-D'
sip=siz(3)
print,'sx=',strtrim(sx,2),'  sca=',strtrim(sca,2),'  lag=',strtrim(lag,2), $
	'  x0=',strtrim(x0,2),'  y0=',strtrim(y0,2)
print,'names:  ',filn,'*'

;Check for previous work and prepare output arrays ----------------

ixylog=fltarr(8,sip)		; initialize output array
orig=0

if (not keyword_set(new)) then begin
	restore,filn+'log'	; here comes old vxylog
	sl=size(vxylog)
	ixylog(*,0:sl(2)-1)=vxylog	; insert old values
	last=vxylog(0,sl(2)-1)		; id of last processed object
   ; where to start new run
	orig=fix(last+1)		; pos of first newly processed object
	print,'Last PG already processed: ',last
	if orig gt sip-1 then message,'All PGs were already processed'
endif

vxylog=ixylog
ixylog=0

;Prepare graphics -------------------------------------------------

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


; BEGIN LOOP OVER OBJECTS -----------------------------------------

for i=orig, sip-1 do begin	

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

pos=inar(1,*,i)		;  positions in 1-D form
pos=pos(where(pos,c))	;  c = number of frames with object
print,' '
print,'Object seq. number:',i
print,'No. of frames =',c,'     Lifetime (minutes) =',(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'

;INW or OUT ?

dsco=sqrt((x(0)-x0)^2+(y(0)-y0)^2)	;orig. dist. from u.center
dscf=sqrt((x(c-1)-x0)^2+(y(c-1)-y0)^2)	;final dist. from u.center
if dsco gt dscf then inw=1 else inw=0	;inw=1..INW, inw=0..OUT or TANG
if inw eq 1 then print,'INWARD moving' else print,'OUTWARD moving'

; 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'

; Average peak intensity

pos=inar(0,*,i)
pos=pos(where(pos))/iph
avin=total(pos)/c

; Output stuff -----------------------------------------------

vxylog(0,i)=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
; vxylog(5,i)                     reserved for speed class
vxylog(6,i)=inw			; INW/OUT ident.
vxylog(7,i)=avin		; average intensity

;------------------------------------------------------------

;Interactive smoothing of trajectory and calculation of speed

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

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)

;Classification - Automatic mode

clas=CLASSIF(v)
wset,0
plot,v,ytitle='velocity magnitude (km/s)'
xyouts,250,300,strtrim(clas,2),charsiz=3,charthick=2,/device,col=150
vxylog(5,i)=clas
print,'Automatic classification OK ? (any_key=yes / no=n)'
if get_kbrd(1) ne 'n' then goto,tplot

;Manual mode

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)'
xyouts,200,300,'MANUAL',charsiz=3,charthick=2,/device,col=150

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

;Trajectory plot --------------------------------------------------

tplot:

print,'Class:',vxylog(5,i)

wset,2
plots,x*reb,y*reb,/device,col=255,thick=1
plots,xs*reb,ys*reb,/device,col=255,thick=2
;head at final position
xf=xs(c-1)*reb		; final positions
yf=ys(c-1)*reb
re2=reb/2
;plots,[xf-re2,xf+re2,xf+re2,xf-re2,xf-re2], $
;      [yf-re2,yf-re2,yf+re2,yf+re2,yf-re2],/device,col=255,thick=2

print,'Reject this object (y), accept (any_key), or STOP (s) ?'
dcs=get_kbrd(1)

case dcs of
  'y':	begin			;destroy it
	  vxylog(1,i)=0
	  print,'  ** Rejected **'
	end
  's':  begin			; destroy last PG and go to fin
	  vxylog(1,i)=0
	  goto,fin
	end
  else: begin
	  save,xs,ys,v,file=filn+strtrim(i,2)
	  print,'  {{ Saved }}'
	end
endcase

plots,x*reb,y*reb,/device,col=0,thick=1
plots,xs*reb,ys*reb,/device,col=0,thick=2
;head at final position
;plots,[xf-re2,xf+re2,xf+re2,xf-re2,xf-re2], $
;      [yf-re2,yf-re2,yf+re2,yf+re2,yf-re2],/device,col=0,thick=2

endfor

; END LOOP OVER OBJECTS --------------------------------------------

fin:
w=where(vxylog(1,*),m)
print,'End of program.'
print,'Scanned: ',strtrim(i+1,2),'   Selected: ',strtrim(m,2),'  objects.'

if m ne 0 then begin
  vxylog=vxylog(*,w)		; compress output
  save,vxylog,file=filn+'log'
endif
wdelete,0,1

end