PRO fcameek,num,start,interval,finish

;IF N_PARAMS() EQ 0 THEN BEGIN
;    doc_library,'fcameek'
;    RETURN
;    ENDIF

patt=STRING(num,FORMAT='("*",I0,"*")')

IF N_ELEMENTS(start) EQ 0 THEN start=0		; hours UT
IF N_ELEMENTS(interval) EQ 0 THEN interval=0.5	; hours
IF N_ELEMENTS(finish) EQ 0 THEN finish=24	; hours UT

tmp=pickfile(/read,path='/jasper/cnsr3_data1/fca',filter=patt)
IF STRLEN(tmp) GT 0 THEN file=tmp ELSE BEGIN
    MESSAGE,'No file found; exiting'
    RETURN
    ENDELSE

; How big is the file?
nlin=nlines(file)

; open the file
OPENR,corrunit,file,/GET_LUN

; read the number of fields of view used
READF,corrunit,nfov
nfov=FIX(ROUND(nfov))
IF nfov NE 3 THEN BEGIN
    MESSAGE,'This analysis only applicable to 3-sensor data!'
    FREE_LUN,corrunit
    RETURN
    ENDIF

; read in the field of view parameters
fovparms=FLTARR(nfov+1,5)
READF,corrunit,fovparms
fovparms=TRANSPOSE(fovparms(1:*,*))

; read in the data
data=FLTARR(nfov+1,nlin-6)
READF,corrunit,data
FREE_LUN,corrunit
ut=TRANSPOSE(data(0,*))
afov=TRANSPOSE(data(1:*,*))

; Determine date parameters from file name
inpt=STRPOS(file,'/fca/')+5
year=1900+FIX(STRMID(file,inpt,2))
month=FIX(STRMID(file,inpt+2,2))
day=FIX(STRMID(file,inpt+4,2))

action=[' - accept data',' - reject data']

xx=REFORM(fovparms(3,*))
yy=REFORM(fovparms(4,*))
PRINT,"'Sensor' coordinates:"
fm='(A,'+STRTRIM(nfov,2)+'F10.2)'
PRINT,"xx:  ",xx,FORMAT=fm
PRINT,"yy:  ",yy,FORMAT=fm

ncross=nfov
good=BYTARR(ncross)	; data quality flags
pktau=FLTARR(ncross)	; lags for maximum cross-correlation
pkrho=pktau		; values of maximum cross-correlation

; Now identify the FOV pairs involved in the cross-correlations
p=INDGEN(ncross)
q=FIX((p+1) MOD nfov)
xi=xx(q)-xx(p)
eta=yy(q)-yy(p)
d=SQRT(xi*xi+eta*eta)			; inter-sensor distances
psi=(ATAN(xi,eta)+2.*!PI) MOD (2.*!PI)	; inter-sensor azimuths
PRINT,'xi:  ',xi,FORMAT=fm
PRINT,'eta: ',eta,FORMAT=fm
PRINT,'d:   ',d,FORMAT=fm
PRINT,'psi: ',psi*!RADEG,FORMAT=fm
PRINT,'Press any key to continue'
kkk=GET_KBRD(1)

; identify first and last UT values
mxut=MAX(ut,MIN=mnut)

; set up uniform UT array for interpolation of observations
nmin=CEIL(mxut/60.)-FLOOR(mnut/60.)+1
uut=(FINDGEN(nmin)+FLOOR(mnut/60.))*60

; interpolate to 1-minute intervals
iafov=FLTARR(nmin,nfov)
FOR i=0,nfov-1 DO iafov(*,i)=interpol(afov(*,i),ut,uut)

; Set up an array to hold the derived parameters
input=FLTARR(9,48)	; Up to 24 hrs of data in half-hour
output=FLTARR(8,48)	; intervals

; Set up a loop through all time intervals between specified limits
first=CEIL(mnut/(interval*3600.)+0.5)>ROUND(start/interval)
last=FLOOR(mxut/(interval*3600.)-1.5)<ROUND(finish/interval)-1
PRINT,first*interval,last*interval
uthlen=interval

; Now do the loop!
FOR hrut=first,last DO BEGIN
good(*)=0B
badcorr=0B

; Specify the interval length and start time
uthstart=hrut*interval
utstart=3600.*uthstart
nsamp=ROUND(uthlen*60.)
IF (nsamp MOD 2) EQ 0 THEN nsamp=nsamp+1
utlen=3600.*uthlen
PRINT,FORMAT='("---------------------------------------------")'
PRINT,utstart,utlen,nsamp $
     ,FORMAT='("Start ",F6.0," s UT; length ",F5.0," s; ",I0," samples")'

; Identify data from the specified interval
w=WHERE((utstart LE uut) AND (uut LE utstart+utlen),nw)

; Check for shallow fading
nofade=0B
FOR i=0,nfov-1 DO nofade=(nofade OR (stdev(iafov(w,i),m) LE 0.02*m))
IF nofade THEN PRINT,'Insufficient fading over this interval!' ELSE BEGIN

    tau=FINDGEN(nsamp)-(nsamp/2)
    auto=FLTARR(nfov,nsamp)
    cross=FLTARR(ncross,nsamp)

; Compute auto- and cross-correlations over this period
    FOR j=0,N_ELEMENTS(tau)-1 DO BEGIN
        utm=utstart+60*tau(j)
        wt=WHERE((utm LE uut) AND (uut LE utm+utlen),nwt)
        IF nw EQ nwt THEN BEGIN
	    FOR i=0,nfov-1 DO auto(i,j)=correlate(iafov(w,i),iafov(wt,i))
	    FOR i=0,ncross-1 DO $
		cross(i,j)=correlate(iafov(w,p(i)),iafov(wt,q(i)))
	    ENDIF ELSE BEGIN
	    MESSAGE,'Not enough data for ' $
		   +STRING(tau(j),FORMAT='(F6.2)')+' minute lag ' $
		   ,/INFORMATIONAL
	    badcorr=1B
	    ENDELSE
        ENDFOR

; Plot out the auto- and cross-correlations
; First the auto-correlations
    IF NOT badcorr THEN BEGIN
        !P.MULTI=[0,1,ncross+2,0,0]
        IF !D.NAME EQ 'X' THEN WINDOW,0,XSIZE=640,YSIZE=860 ELSE ERASE
        tstring=STRING(ymd2date(year,month,day,format="d$ n$ y$") $
		      ,uthstart,uthlen $
                      ,FORMAT='(A,"  630 nm airglow FCA: start "' $
			     +',F5.2,", len ",F4.2)')
        PLOT,tau,auto(0,*),XTITLE='Time lag (minutes)' $
	                  ,YTITLE='Autocorrelation' $
		          ,YRANGE=[0,1.0],/YSTYLE,TITLE=tstring $
			  ,LINESTYLE=1
        FOR i=1,nfov-1 DO OPLOT,tau,auto(i,*),LINE=i+1
	mauto=REFORM(TOTAL(auto,1)/nfov)
	OPLOT,tau,mauto
        autofit=gaussfnb(tau,mauto,fp)
        OPLOT,tau,autofit,PSYM=1

        w0=(WHERE(tau EQ 0))(0)

        PLOT,tau,mauto,PSYM=1 $
            ,XTITLE='Time lag (minutes)' $
            ,YTITLE='Gaussian Fit to Mean Autocorrelation' $
	    ,YRANGE=[-1,1],/YSTYLE
        xl=!X.CRANGE(0)
        xr=!X.CRANGE(1)
        xt=FINDGEN(2*(xr-xl)+1)*0.5+xl

; Determine the y-axis offset in data coordinates between lines of text
	corner=[[0,0],[!D.X_CH_SIZE,!D.Y_CH_SIZE]]
	sizech=CONVERT_COORD(corner,/DEVICE,/TO_DATA)
	yspace=1.25*(sizech(1,1)-sizech(1,0))

; Plot and label the Gaussian fit to the mean autocorrelation function
        zf=(xt-fp(1))/fp(2)
        OPLOT,xt,fp(0)*EXP(-zf^2/2),THICK=2
        XYOUTS,xl+2,1.-yspace $
	      ,STRING(fp(0),FORMAT='("Height: ",F6.2)'),ALIGNMENT=0
        XYOUTS,xl+2,1.-2*yspace $
	      ,STRING(fp(1),FORMAT='("Center: ",F6.2)'),ALIGNMENT=0
        XYOUTS,xl+2,1.-3*yspace $
	      ,STRING(fp(2),FORMAT='("Sigma:  ",F6.2)'),ALIGNMENT=0

; Plot and label each ofthe cross-correlation functions
        FOR i=0,ncross-1 DO BEGIN
	    yt=STRING(p(i),q(i) $
		     ,FORMAT='("Cross-correlation(fov ",I1,"; fov ",I1,")")')
	    tt=STRING(p(i),q(i) $
		     ,FORMAT='("<f(fov ",I1,",t) * f(fov ",I1,",t+!7s!x)")')
	    PLOT,tau,cross(i,*),XTITLE='Time lag !7s!x (minutes)',YTITLE=yt $
		               ,YRANGE=[0,1.0],/YSTYLE,TITLE=tt
	    pkpq=extremes(cross(i,*),1)	; find maxima in cross-correlation
	    IF pkpq(0) EQ -1 $
		THEN MESSAGE,'No peaks in '+STRTRIM(p(i),2) $
		            +'-'+STRTRIM(q(i),2)+' cross-correlation' $
		            ,/INFORMATIONAL $
                ELSE BEGIN
	        wgthalf=WHERE(cross(i,pkpq) GE 0.5,nwgthalf)	; > 0.5
	        IF nwgthalf EQ 0 $
	            THEN MESSAGE,'No significant peaks in '+STRTRIM(p(i),2) $
			        +'-'+STRTRIM(q(i),2)+' cross-correlation' $
		                ,/INFORMATIONAL $
	            ELSE BEGIN
	            pkpq=pkpq(wgthalf)	; select peaks with cross-corr > 0.5
	            IF N_ELEMENTS(pkpq) GT 1 THEN BEGIN	; multiple peaks:
	                mintau=MIN(ABS(pkpq-w0),wmin)	; choose peak
	                pkpq=pkpq(wmin)			; with min. lag
	                ENDIF ELSE pkpq=pkpq(0)
	            IF (pkpq LT 2) OR (nsamp-3 LT pkpq) THEN BEGIN
	                PRINT,'Lag for peak cross-correlation too close ' $
			     +'to maximum |lag|.'
	                PRINT,'No quadratic fit attempted.'
	                ENDIF ELSE BEGIN
	                wpk=LINDGEN(5)+pkpq-2	; 5 pts centered on peak
	                yf=1
	                pkfit=svdfit(tau(wpk) $	; Fit a parabola to the peak
				    ,REFORM(cross(i,wpk)),3,yfit=yf)
	                pkpq=-pkfit(1)/(2.*pkfit(2)); lag at peak
		        pktau(i)=pkpq
		        pkrho(i)=poly(pkpq,pkfit)	; peak cross-correlation
		        good(i)=1B
	                PRINT,'Peak '+STRTRIM(p(i),2)+'-' $
			      +STRTRIM(q(i),2) $
			      +' cross-correlation at ' $
			     ,STRTRIM(pkpq,2),' min. lag'
                        OPLOT,tau(wpk),yf,PSYM=1
	                algn=(wpk(2) LT nsamp/2)
	                xv=xl+2+algn*(xr-xl-4)
		        pktxt=STRING(pkrho(i),pkpq $
		             ,FORMAT='("Peak ",F5.3," at ",F7.3," min")')
	                XYOUTS,xv,1.-0.5*yspace,pktxt,ALIGNMENT=algn
	                ENDELSE
		    ENDELSE
	        ENDELSE
            ENDFOR

	IF TOTAL(good) EQ ncross THEN BEGIN
	    PRINT,'pktau: ',pktau
	    PRINT,'pkrho: ',pkrho
	    ntd=ABS(TOTAL(pktau))/TOTAL(ABS(pktau))
	    PRINT,STRING(ntd,action(ntd GE 0.2) $
	         ,FORMAT='("Normalized time discrepancy = ",F5.3,A)')

; Lags for peak correlation fairly consistent - carry out analysis
            IF ntd LT 0.2 THEN BEGIN
; Determine 'width' of mean autocorrelation from average |lag| for
; 60.7% correlation.  First find the 60.7% correlation points.
		over=REPLICATE(-1,nsamp)
		over(WHERE(mauto GT 0.5))=1	; > 60.7% correlation
		cross=WHERE(over*SHIFT(over,-1) EQ -1)
		IF N_ELEMENTS(cross) GT 2 THEN BEGIN	; keep only 2!
		    abslag=ABS(cross-nsamp/2)	; |lag|
		    sabslag=SORT(abslag)	; sorting array
		    cross=cross(sabslag(0:1))	; 2 crosses at min |lag|
		    IF cross(0) GT cross(1) THEN cross=reverse(cross)
		    ENDIF
; Interpolate to find lags at 60.7% autocorrelation
		slope=(mauto(cross+1)-mauto(cross))/(tau(cross+1)-tau(cross))
		tau61=(REPLICATE(EXP(-0.5),2)-mauto(cross))/slope+tau(cross)
                tauarg1=mean(ABS(tau61))	; use mean of both |lags|
		PRINT,tau61 $
		     ,FORMAT='("60.7% correlation lags: ",2F6.2," min.")'
		PRINT,tauarg1 $
		     ,FORMAT='("Mean 60.7% correlation lag: ",F6.2," min.")'
		FOR i=0,1 DO PLOTS,[tau61(i),tau61(i)],[4.8,5.8]
		input(*,hrut)=[utstart,utlen,tauarg1,pktau,pkrho]

; Fit to apparent velocity in order to redefine lags for peak
; cross-correlation (Meek, 1980, eqn. (5))
	        opktau=pktau		; save old peak lags
	        alt=regress([TRANSPOSE(xi),TRANSPOSE(eta)] $
		           ,pktau,pkrho*pkrho $
		           ,altfit,const,altsig,ftest,r,rmul,chisq)
                pktau=alt(0)*xi+alt(1)*eta	; redefine peak lags
	        FOR i=0,ncross-1 DO PLOTS,[opktau(i),opktau(i)] $
					 ,[0.,3.4-1.2*i]
	        FOR i=0,ncross-1 DO PLOTS,[pktau(i),pktau(i)] $
					 ,[0.,3.4-1.2*i] $
				         ,THICK=2

; Now solve for pattern parameters using eqns (6) - (17) from Meek
                bigQ=-2*ALOG(0.5)/(tauarg1*tauarg1*3600)	; eqn (9)
		mi=(-2*ALOG(pkrho)+bigQ*pktau*pktau*3600) $	; eqn (13a)
		   /(d*d)
                mij=mi(q)-mi(p)
		tan2theta=-(mij(1)*COS(2*psi(0)) $	; eqn (14a)
			   +mij(2)*COS(2*psi(1)) $
			   +mij(0)*COS(2*psi(2))) $
			  /(mij(1)*SIN(2*psi(0)) $
			   +mij(2)*SIN(2*psi(1)) $
			   +mij(0)*SIN(2*psi(2)))
		theta=ATAN(tan2theta)/2.		; in radians
		tilt=(!RADEG*theta+180.) MOD 180.	; in degrees
		rf=(mi(0)*COS(psi(1)-theta)^2 $		; eqn (14b)
		   -mi(1)*COS(psi(0)-theta)^2)/mij(0)
                r2=rf/(1+rf)				; must be > 0
		IF r2 LE 0 THEN BEGIN
		    MESSAGE,'Negative squared axial ratio!' $
			   ,/INFORMATIONAL
		    PRINT,r2,FORMAT='("r^2: ",E11.4)'
		    ENDIF ELSE BEGIN
		    b2=1.5*(1./r2+1)/TOTAL(mi)		; eqn (14c)
		    a2=b2*r2				; eqn (14d)
		    IF (a2 LE 0) OR (b2 LE 0) THEN BEGIN
			MESSAGE,'Negative a^2 or b^2!' $
			       ,/INFORMATIONAL
			PRINT,a2,FORMAT='("a^2: ",E11.4)'
			PRINT,b2,FORMAT='("b^2: ",E11.4)'
                        ENDIF ELSE BEGIN
			tanphimth=((pktau(0)*d(1)*COS(psi(1)-theta) $
				   -pktau(1)*d(0)*COS(psi(0)-theta)) $
				  /(pktau(1)*d(0)*SIN(psi(0)-theta) $
				   -pktau(0)*d(1)*SIN(psi(1)-theta))) $
				  *(b2/a2)		; eqn (15)
			phimth=ATAN(tanphimth)
			phi=phimth+theta
			dir=(!RADEG*phi+360.) MOD 360.
			bigV=(bigQ*pktau(0)*60/d(0)) $	; eqn (16)
			    /((COS(phimth)*COS(psi(0)-theta)/a2) $
			     +(SIN(phimth)*SIN(psi(0)-theta)/b2))
			IF bigV LT 0 THEN BEGIN		; V < 0 ?
			    bigV=-bigV			; change sign
			    dir=(dir+180.) MOD 360.	; and direction
			    ENDIF
			invc2=bigQ-bigV*bigV $		; eqn (17)
				  *(((COS(phimth)^2)/a2) $
				   +((SIN(phimth)^2)/b2))
			IF invc2 LE 0 THEN BEGIN
			    MESSAGE,'Negative c^2!' $
				   ,/INFORMATIONAL
			    PRINT,1./invc2,FORMAT='("c^2: ",E11.4)'
			    ENDIF ELSE BEGIN
			    c2=1./invc2
			    a=SQRT(a2)
			    b=SQRT(b2)
			    IF a LT b THEN BEGIN	; A < B ?
				swap,a,b
				tilt=(tilt+90.) MOD 180.
				ENDIF
			    c=SQRT(c2)
			    PRINT,a,FORMAT='("Major axis: ",F6.1," km; ",$)'
			    PRINT,b,FORMAT='("Minor axis: ",F6.1," km; ",$)'
			    PRINT,tilt $
				 ,FORMAT='("Long axis azimuth: ",F6.1," deg")'
			    PRINT,bigV,FORMAT='("Velocity: ",F6.3," km/s ",$)'
			    PRINT,dir,FORMAT='("toward azimuth ",F6.1," deg")'
			    PRINT,c/60. $
				 ,FORMAT='("Characteristic time: ",F5.1," min")'
			    output(*,hrut)=[utstart,utlen,a,b,tilt $
					   ,bigV,dir,c/60.]

; Now check the previously computed values of various intermediate
; parameters by recomputing them form the subsequently derived
; pattern parameters.
;			    PRINT,FORMAT='()'
;			    phimth=!DTOR*(dir-tilt)
;			    theta=!DTOR*tilt
;			    a2=a*a
;			    b2=b*b
;			    r2=a2/b2
;			    Qchk=bigV*bigV $
;				*(((COS(phimth)^2)/a2) $
;				 +((SIN(phimth)^2)/b2))+invc2
;			    PRINT,bigQ,Qchk $
;				 ,FORMAT='("Q: ",E11.4,"; check: ",E11.4)'
;			    tmaxi=(bigV/bigQ)*d $
;				 *((COS(phimth)*COS(psi-theta)/a2) $
;				  +(SIN(phimth)*SIN(psi-theta)/b2))
;			    PRINT,pktau*60 $
;				 ,FORMAT='("Peak lags: ",3F7.1)'
;			    PRINT,tmaxi $
;				 ,FORMAT='("Check:     ",3F7.1)'
;			    tmp=d*d $
;			       *(((COS(psi-theta)^2)/a2) $
;				+((SIN(psi-theta)^2)/b2)) $
;			       -bigQ*(pktau*60)^2
;			    rhomaxi=EXP(-0.5*tmp)
;			    PRINT,pkrho $
;				 ,FORMAT='("Peak correlations: ",3F6.3)'
;			    PRINT,rhomaxi $
;				 ,FORMAT='("Check:             ",3F6.3)'
;			    michk=((COS(psi-theta)^2) $
;				  +r2*(SIN(psi-theta)^2)) $
;				 /(r2*b2)
;			    PRINT,mi $
;				 ,FORMAT='("m parameters: ",3E12.4)'
;			    PRINT,michk $
;				 ,FORMAT='("Check:        ",3E12.4)'
			    ENDELSE
			ENDELSE
		    ENDELSE
                ENDIF
	    ENDIF
	ENDIF
    ENDELSE
ENDFOR

; Now save the analysis input data to a file
datfile=STRMID(file,0,STRLEN(file)-4)+'dat'
OPENW,datunit,datfile,/GET_LUN
PRINTF,datunit,input,FORMAT='(2F7.0,4F7.2,3F7.3)'
FREE_LUN,datunit

; Now save the analysis results to a file
fcafile=STRMID(file,0,STRLEN(file)-4)+'fca'
OPENW,fcaunit,fcafile,/GET_LUN
PRINTF,fcaunit,output,FORMAT='(5F7.0,F7.2,2F7.0)'
FREE_LUN,fcaunit

;IF !D.NAME EQ 'X' THEN WDELETE,0
END