pro vbe,ffile,result

cd,'/scratch/sob/lap93'

scale=0.06218*725    ; km/pixel

restore,'vy08'    ; logfile of 20 selected images

; read background files
backarr=intarr(4,20,3)  ;(seq.,x,y,int) x 20 frames x 3 backgr.points
cb=0	;counter of backgr. points
back:
read,'background ident.no.? (if 0, quit)',nb
if fix(nb) eq 0 then goto,further       ;without background
openr,1,'b'+strtrim(fix(nb),2)+'.dat'
readf,1,nn		;no. of real frames
dumb=intarr(4,nn)
readf,1,dumb
close,1
backarr(0,0,cb)=dumb	;array of background data, cb = 0 or 1 or 2
cb=cb+1
if cb gt 2 then goto,further
goto,back

further:
; read feature file
openr,1,ffile
readf,1,n		;no. of real frames
farr=intarr(4,n) 	;(seq.,x,y,intensity) for each frame
readf,1,farr
close,1

; calculations

; time differences in seconds (time [hours] in vy08(1,*))
tim=(vy08(1,farr(0,0):farr(0,n-1))-vy08(1,0))*3600.
dtim=tim-shift(tim,1) & dtim(0)=0.	;time differences

;arrays of displacements Dx, Dy, Dr (in km)
Dx=(farr(1,*)-shift(farr(1,*),0,1))*scale & Dx(0)=0
Dy=(farr(2,*)-shift(farr(2,*),0,1))*scale & Dy(0)=0
Dr=sqrt(Dx*Dx+Dy*Dy)

; velocities vx, vy, v
vx=Dx/dtim
vy=Dy/dtim


; FOURIER FILTRATION

;  interpolation to 512 equidistant samples in time
dti=(tim(n-1)-tim(1))/512. 	; step in time
tint=findgen(513)*dti+tim(1)    ; tint(0)=tim(1) !
vxi=fltarr(513)
vyi=vxi
vxi(0)=vx(1)
vyi(0)=vy(1)
for i=1,512 do begin		; interpolation of vx,vy
   j=max(where(tint(i) gt (tim+dti/10.)))  ; takes care of trunc.errors
   vxi(i)=vx(j)+(vx(j+1)-vx(j))/(tim(j+1)-tim(j))*(tint(i)-tim(j))
   vyi(i)=vy(j)+(vy(j+1)-vy(j))/(tim(j+1)-tim(j))*(tint(i)-tim(j))
endfor

filter:
plot,tim,vx,yra=[-2.,2.]
oplot,tim,vy,linestyle=2
read,'filter cutoff (pixel position):  ',cut

u=fltarr(513)
for i=0,cut do u(i)=(1.-cos(!pi*(cut-i)/cut))/2.

u=[u,reverse(u(1:511))]		; cosinus filter

fvxi=fft([vxi,reverse(vxi(1:511))],-1)*u
fvyi=fft([vyi,reverse(vyi(1:511))],-1)*u
vxif=fft(fvxi,1) & vxif=float(vxif(0:512))
vyif=fft(fvyi,1) & vyif=float(vyif(0:512))
oplot,tint,vxif
oplot,tint,vyif,linestyle=2

print,'smoothing ok? (y/n)'
if get_kbrd(1) eq 'n' then goto,filter

;  interpolation back to original sampling
vx1=vx & vy1=vy
vx1(1)=vxif(0)
vy1(1)=vyif(0)
vx1(0)=vx1(1)  ; replication of the first values (no physical sense)
vy1(0)=vy1(1)

for i=2,n-1 do begin
   j=max(where(tim(i) gt tint))
   vx1(i)=vxif(j)+(vxif(j+1)-vxif(j))/(tint(j+1)-tint(j))*(tim(i)-tint(j))
   vy1(i)=vyif(j)+(vyif(j+1)-vyif(j))/(tint(j+1)-tint(j))*(tim(i)-tint(j))
endfor

; computation of errors
evx=0. & evy=0.
for i=1,n-1 do begin
   evx=evx+(vx(i)-vx1(i))*(vx(i)-vx1(i))
   evy=evy+(vy(i)-vy1(i))*(vy(i)-vy1(i))
endfor
evx=evx/((n-1)*(n-2))
evy=evy/((n-1)*(n-2))

print,'errors in vx and vy:  ',sqrt(evx),sqrt(evy)

vx=vx1 & vy=vy1

oplot,tim,vx
oplot,tim,vy,linestyle=2

v=sqrt(vx*vx+vy*vy)

print,'error in v:'
print,sqrt((vx/v)^2*evx+(vy/v)^2*evy)

v1=Dr/dtim
window,2
plot,tim,v1
oplot,tim,v,linestyle=2
wset,0

;intensities
int=farr(3,*)

; P/B ratio
iback=fltarr(n)+1.
if cb eq 0 then goto,fin

back=intarr(4,n,cb)
for i=0,cb-1 do back(*,*,i)=backarr(*,farr(0,0):farr(0,n-1),i)

; computation of local background as weighted mean, where
;  the weights are the inverse of distances.

dd=fltarr(n)     ; distance
ddsum=fltarr(n)
iback=fltarr(n)

for i=0,cb-1 do begin
	dd=1./(sqrt((farr(1,*)-back(1,*,i))^2+(farr(2,*)-back(2,*,i))^2))
	iback=iback+back(3,*,i)*dd
	ddsum=ddsum+dd
endfor

iback=iback/ddsum

fin:
pb=int/iback

result=fltarr(12,n)
result(0,*)=farr(0,*)	; seq.no.
result(1,*)=farr(1,*)	; x-position
result(2,*)=farr(2,*)	; y-position
result(3,*)=Dx		; displacements
result(4,*)=Dy
result(5,*)=Dr
result(6,*)=tim/60.     ; time in minutes from the beginning of series
result(7,*)=vx		; velocities (filtered)
result(8,*)=vy
result(9,*)=v
result(10,*)=int	; intensity
result(11,*)=pb		; P/B ratio

end