PRO flow,date

;Calculation of LCT flow maps based on destretched, subsonic
;filtered and registered WL and IR movies. WL movies have to
;be degraded to the resolution of IR ones.
;Basic parameters: LCT Gaussian window FWHM = 9 pix (1.49")
;		   integration length 29 frames (10 minutes)
;		   lag=1 frame (21 s), 1pix/lag = 5.724 km/s
;INPUT: date = 23 or 25
;	files irdata/ifi*.bin, wldata/wsi*.bin
;OUTPUT: SAVE file containing 3-D arrays. The 3-rd dimension
;	is used to store multiple 2-D arrays corresponding to
;	different tracking intervals. Apodized edges are cut away.
;	ARRAYS: vx, vy - x,y velocitry components for IR (km/s)
;		wx, wy - x,y velocitry components for WL (km/s)
;		vm, wm - speed magnitudes for IR resp. WL (km/s)
;		ave, awe - images averaged over tracking
;			intervals (IR and WL)
;CALLED: flowmakr and its subroutines

on_error,1
if (date ne 23) and (date ne 25) then  $
   message,'Sorry, I do not know how to handle this dataset'
dat=strtrim(date,2)

;PARAMETERS

scale=5.724	;velocity scale km/s per pix/lag

if date eq 23 then begin

  si=[3,408,320,118]
  n=4	;number of tracking intervals (in 1 movie)
	;limits of tracking intervals (frames)
  bo=[[5,33],[31,59],[57,85],[84,112]]
  					
endif else begin

  si=[3,388,300,150]
  n=5	;number of tracking intervals (in 1 movie)
	;limits of tracking intervals (frames)
  bo=[[7,35],[34,62],[61,89],[87,115],[114,142]]

endelse

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

;ARRAYS

vx=fltarr(si(1)-20,si(2)-20,2*n)	;cut away apodized edges
vy=fltarr(si(1)-20,si(2)-20,2*n)	; (10 pix each side)
wx=fltarr(si(1)-20,si(2)-20,2*n)
wy=fltarr(si(1)-20,si(2)-20,2*n)
vm=fltarr(si(1)-20,si(2)-20,2*n)
wm=fltarr(si(1)-20,si(2)-20,2*n)
ave=fltarr(si(1)-20,si(2)-20,2*n)
awe=fltarr(si(1)-20,si(2)-20,2*n)
mov=intarr(si(1),si(2),si(3))

;RUNS

; IR flowmaps -----------------------

openr,1,'../irdata/ifi'+dat+'a.bin'
readu,1,mov
close,1
print,' Processing 1st part, IR'

for i=0,n-1 do begin

	sub=mov(10:si(1)-11,10:si(2)-11,bo(0,i):bo(1,i))
	av=total(sub,3)/29.
	ave(0,0,i)=av			;average image

	flowmakr,sub,1,9,1,vx1,vy1,/qfit2
	vm1=sqrt(vx1*vx1+vy1*vy1)
	vx(0,0,i)=vx1
	vy(0,0,i)=vy1
	vm(0,0,i)=vm1

endfor

openr,1,'../irdata/ifi'+dat+'b.bin'
readu,1,mov
close,1
print,' Processing 2nd part, IR'

for i=0,n-1 do begin

	sub=mov(10:si(1)-11,10:si(2)-11,bo(0,i):bo(1,i))
	av=total(sub,3)/29.
	ave(0,0,i+n)=av			;average image

	flowmakr,sub,1,9,1,vx1,vy1,/qfit2
	vm1=sqrt(vx1*vx1+vy1*vy1)
	vx(0,0,i+n)=vx1
	vy(0,0,i+n)=vy1
	vm(0,0,i+n)=vm1

endfor

; WL flowmaps -----------------------

openr,1,'../wldata/wsi'+dat+'a.bin'
readu,1,mov
close,1
print,' Processing 1st part, WL'

for i=0,n-1 do begin

	sub=mov(10:si(1)-11,10:si(2)-11,bo(0,i):bo(1,i))
	av=total(sub,3)/29.
	awe(0,0,i)=av			;average image

	flowmakr,sub,1,9,1,vx1,vy1,/qfit2
	vm1=sqrt(vx1*vx1+vy1*vy1)
	wx(0,0,i)=vx1
	wy(0,0,i)=vy1
	wm(0,0,i)=vm1

endfor

openr,1,'../wldata/wsi'+dat+'b.bin'
readu,1,mov
close,1
print,' Processing 2nd part, WL'

for i=0,n-1 do begin

	sub=mov(10:si(1)-11,10:si(2)-11,bo(0,i):bo(1,i))
	av=total(sub,3)/29.
	awe(0,0,i+n)=av			;average image

	flowmakr,sub,1,9,1,vx1,vy1,/qfit2
	vm1=sqrt(vx1*vx1+vy1*vy1)
	wx(0,0,i+n)=vx1
	wy(0,0,i+n)=vy1
	wm(0,0,i+n)=vm1

endfor

;SCALING to km/s and SAVING

vx=vx*scale
vy=vy*scale
vm=vm*scale
wx=wx*scale
wy=wy*scale
wm=wm*scale

print,' Saving results into: ','flows'+dat
save,vx,vy,vm,ave,wx,wy,wm,awe, file='flows'+dat

END