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 by Michal. ;22.5.97 adapted for penumbra data by P. Brandt ;15.10.97 size criterium added by Michal on_error,1 ;*********** PARAMETERS ******************************* sx = 200 ;x-size of frames processed by OTRACK4 vmax=3 ;2 ;maximum velocity (km/s) dtim=52. ;time difference between frames (s) scal=0.0835 ;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