PRO instvel4,tlim,vxylogi,vxylogo,new=new

; A program to obtain instantaneous velocities and trajectories
; of INW/OUT PGs from the object-tracking result array.
; Version 4 - adopted from ver.3 for Suetterlin's data
;	      (trajectory smoothed by SMSPLINE,
;	      simultaneous processing of INW/OUT PGs)
;	      semi-automatical classification of speed shapes
;		(function CLASSIF.PRO)
;             3 Apr 2000

;INPUTS: TLIM - minimum lifetime of objects (in integer minutes)
;	 inar - 3D array, result of object tracking (all quadrants),
;		composed od triml,trimt,trimr,trimb (RESTORE)
;	 ima  - underlying image--average of the series (RESTORE)
;        
;KEYWORD NEW  - if set, new output arrays and log files are created.
;		if not, old log files are restored and the work
;		continues with next objects

;OUTPUTS: VXYLOGI, VXYLOGO - 2-D arrays containing data
;	  for INW resp. OUT accepted objects:
;	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.
;	6. time-averaged intensity
;  These arrays are saved in log files vxy[tlim][i/o].log
;
;  Save-files  vxy[tlim][i/o].nn , one per accepted object.
;	Each file contains arrays xs,ys (smoothed trajectory)
;	and v (instantaneous velocity magnitude).

on_error,1

;Parameters --------------------------------------------------------

	reb=1			;  rebin factor for visualization
	sca=0.0815		;  arcsec per pixel
	lag=30		;  time between 2 frames (seconds)
	x0=220		;  x-position of umbral center (pix)
	y0=245		;  y-position of umbral center (pix)
	iph=1218.	;  mean quiet photosphere intensity

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

;Inputs, filenames, etc.

RESTORE,'avemov'
ima=avemov
avemov=0
im88=0
RESTORE,'trim'
inar=merge(triml,trimt)
inar=merge(inar,trimr)
inar=merge(inar,trimb)	; putting together all quadrants
triml=0
trimt=0
trimr=0
trimb=0

fili='vxy/vxy'+strtrim(tlim,2)+'i.'
filo='vxy/vxy'+strtrim(tlim,2)+'o.'

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

;Selection of objects with lifetimes ge tlim ----------------------

lf=intarr(siz(3))
for i=0,siz(3)-1 do begin
	w=where(inar(0,*,i),cnt)
	lf(i)=cnt		; no. of frames where object is present
endfor
lf=(lf-1)*lag/60.		; lifetimes in minutes (float)

pgno=where(lf ge tlim,sip)	; list of subscripts of selected objects

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

ixylogi=fltarr(7,sip)		; initialize output array (INW)
ixylogo=fltarr(7,sip)		; initialize output array (OUT)
orig=0

if (not keyword_set(new)) then begin
	restore,fili+'log'
	restore,filo+'log'	; here come old vxylogi, vxylogo
	sli=size(vxylogi)
	slo=size(vxylogo)
	ixylogi(*,0:sli(2)-1)=vxylogi	; insert old values
	ixylogo(*,0:slo(2)-1)=vxylogo
	lasti=vxylogi(0,sli(2)-1)	; id of last processed INW PG
	lasto=vxylogo(0,slo(2)-1)	; id of last processed OUT PG
   ; where to start new run
	last=max(lasti,lasto)		; id of last processed PG
	orig=where(pgno eq last)+1	; pos of first processed PG in pgno
	orig=orig(0)
	print,'Last PG already processed: ',last
	if orig gt sip-1 then message,'All PGs were already processed'
endif
vxylogi=ixylogi
vxylogo=ixylogo
ixylogi=0
ixylogo=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

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 (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,*,pgno(i))
pos=pos(where(pos))/iph
avin=total(pos)/c

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

if inw eq 1 then begin			;INW

	vxylogi(0,i)=pgno(i)		; object number
	vxylogi(1,i)=(c-1)*lag/60.	; lifetime (minutes)
	vxylogi(2,i)=dst		; travel distance
	vxylogi(3,i)=cdst		; distance from UC center
	vxylogi(4,i)=avvm		; average speed
	vxylogi(6,i)=avin		; average intensity

endif else begin			;OUT

	vxylogo(0,i)=pgno(i)		; object number
	vxylogo(1,i)=(c-1)*lag/60.	; lifetime (minutes)
	vxylogo(2,i)=dst		; travel distance
	vxylogo(3,i)=cdst		; distance from UC center
	vxylogo(4,i)=avvm		; average speed
	vxylogo(6,i)=avin		; average intensity

endelse

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

;Interactive smoothing of trajectory and calculation of speed

;Initial parameters for SMSPLINE
tens=5			;tension (decreases with value)
nods=4  		;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

if inw eq 1 then vxylogi(5,i)=clas else vxylogo(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 $
   if inw eq 1 then vxylogi(5,i)=nods else vxylogo(5,i)=nods

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

tplot:

if inw eq 1 then print,'Class:',vxylogi(5,i)  $
	else print,'Class:',vxylogo(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
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
	if inw eq 1 then vxylogi(1,i)=0 else vxylogo(1,i)=0
	print,'  ** Rejected **'
	end
  's':  begin			; destroy last PG and go to fin
	if inw eq 1 then vxylogi(1,i)=0 else vxylogo(1,i)=0
	goto,fin
	end
  else: begin
	if inw eq 1 then save,xs,ys,v,file=fili+strtrim(pgno(i),2) $
	else save,xs,ys,v,file=filo+strtrim(pgno(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:
wi=where(vxylogi(1,*),mi)
wo=where(vxylogo(1,*),mo)
print,'End of program.'
print,'Scanned: ',strtrim(i,2),' PGs.  Selected: ',strtrim(mi,2), $
      ' INW  and ',strtrim(mo,2),' OUT'
if mi ne 0 then begin
  vxylogi=vxylogi(*,wi)			; compress output
  save,vxylogi,file=fili+'log'
endif
if mo ne 0 then begin
  vxylogo=vxylogo(*,wo)			; compress output
  save,vxylogo,file=filo+'log'
endif
wdelete,0,1

end