PRO otrim,datin,inlf,lim,msz,datout

;OTRIM - a program to revise the results of OTRACK2:
;Objects with lifetime shorter than parameter LIM,
;with area smaller than parameter MSZ, and
;with velocities larger than VMAX are eliminated.
;
;INPUTS:  DATIN - 3-D data volume created by OTRACK2
;	   ([intens., pos., size] * frames * objects)
;	  INLF  - 2-D array with lifetimes produced by OTRACK2
;	  LIM - minimum accepted lifetime (in frame lags)
;	  MSZ - minimum accepted size (area) in pixels
;OUTPUT:  DATOUT - 3-D volume with revised data.
;
;18.6.96, checked with OTROK, OTROK2.
;12.7.96, changed velocity determination
;Michal.
;22.5.97 adapted for penumbra data by P. Brandt
;15.10.97 size criterion added

on_error,1

;*********** PARAMETERS *******************************
sx = 180	;x-size of frames processed by OTRACK2
vmax= 1  ;3 	;maximum velocity (km/s)
dtim=20.	;time difference between frames (s)
scal=0.062	;scale arcsec/pix
;******************************************************
print,'OTRIM: sx =',sx,'     minlife =',lim,'     minsize =',msz

kms=scal*725.			;conversion pixels to km

;Eliminate UDs with lifetimes below the limit

slf=size(inlf)
if slf(0) eq 2 then w=where(inlf(0,*) ge lim,n)
if slf(0) eq 1 then w=where(inlf ge lim,n)
if slf(0) gt 2 then message,'****wrong input array inlf***'	
					;n - new no. of objects
print,'Original size:',slf(0:2)
print,'Lifetime crit.: selected',n

datout=datin(*,*,w)

;Eliminate objects faster than vmax and smaller than msz

v=fltarr(n)			;velocity
asz=fltarr(n)			;area

;Loop over objects
for i=0,n-1 do begin
  w=where(datout(0,*,i),m)	;subscripts of existing data,
				;m-1 is the lifetime
  ;2-D position
  pos=datout(1,w,i)	;1-D coordinate
  x=pos mod sx		;2-D coords
  y=pos/sx

  ;Average velocity
  ;Linear least squares fit (formulae see Rektorys 1973, p.1065)
  ;Mean velocity is defined as a slope of x(t) and y(t)

  t=dindgen(m)*dtim		;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/m
  ssy=suty-sut*suy/m
  sst=sut2-sut*sut/m
  vx=float(ssx/sst)*kms
  vy=float(ssy/sst)*kms
  v(i)=sqrt(vx*vx+vy*vy)

  ;Average object area
  asz(i)=total(datout(2,w,i))/m

endfor

w=where(v le vmax,n)
print,'Velocity crit.: selected',n
w=where((v le vmax) and (asz ge msz),n)
datout=datout(*,*,w)
print,'Velocity and Size crit.: selected',n

END