;+
; NAME:
;       FCAMEEK
; PURPOSE:
;       Full correlation analysis of Polar Camera image data following
;       Meek, JATP, 1980.
; CATEGORY:
;       Image analysis
; CALLING SEQUENCE:
;       FCAMEEK,YEAR,MONTH,DAY,CAM,FILT,NUM,START,INTERVAL,STEP,FINISH $
;              , Wait = wait, PS = ps, Smooth = smooth, UseFiles = usefiles
; INPUTS:
;       YEAR, MONTH: Specify the (start) year and month of the
;               observations.  If only these two parameters are supplied
;               the routine sets up for multiple-day analysis.
; OPTIONAL INPUTS:
;       DAY:    Specifies the single UT day of the observations.
;       CAM, FILT: Specify the camera and filter numbers.
;       NUM:    The number of the file specifying the data to be
;               analyzed.
;       START:  The start time (UT hours). Default = 0.
;       INTERVAL: The length of the analysis interval (hours). Default =
;               0.5.
;       STEP:   The time between start times for successive analyses.
;               Default = INTERVAL.
;       FINISH: The end time (UT hours). Default = 24.  Use values > 24
;               if multi-day data are to be analyzed.
; KEYWORD PARAMETERS:
;       /WAIT:  If set, pauses after each plot until a key is pressed.
;       PS:     If set, specifies the name of the PostScript file to
;               which output is directed to for subsequent plotting.  
;               Note: this keyword and the WAIT keyword are mutually 
;               exclusive; if PS has a value, WAIT is set to 0.
;       Smooth: If set, specifies the width (in minutes) of the running 
;               average filter to be applied to the data for smoothing 
;               purposes before the FCA is applied.  If set as /Smooth,
;               a value of 5 minutes is used.
;       UseFiles: Should be set equal to a StrArr of fully-qualified path 
;               names, in order to use those files in the subsequent 
;               analysis.  This keyword allows the routine to be run in 
;               'batch' mode, rather than interactively.
; OUTPUTS:
;
; OPTIONAL OUTPUTS:
;       None.
; COMMON BLOCKS:
;       None.
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; EXAMPLE:
;
; SEE ALSO:
;
; MODIFICATION HISTORY:
;       Written by:     DPSteele, ISR, October 1993.
;       Documented by:  DPSteele, ISAS, May 1996.
;-
PRO fcameek,year,month,day,cam,filt,num,start,interval,step,finish $
	   ,wait=wait,ps=ps,smooth=sm,usefiles=usf

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

; Set up OS-dependent parameters.
@isitdos

; ON_ERROR,2

IF N_Elements(ps) GT 0 THEN BEGIN
    IF DataType(ps) NE 'STR' THEN BEGIN
	Help,ps
	Print,DataType(ps)
	Read,'Enter the name of the PostScript file to write: ',ps
	ENDIF
    psname=StrTrim(ps,2)
    wait=0
    ENDIF
    
IF N_Elements(sm) EQ 0 THEN sm=0
IF sm EQ 1 THEN sm=5

IF N_Elements(usf) GT 0 THEN BEGIN      ; validate supplied file names
    nusf=N_Elements(usf)
    valid=BytArr(nusf)
    FOR i=0,nusf-1 DO valid(i)=(rfindfile(usf(i)) NE '')
    IF Total(valid) NE nusf THEN BEGIN
	FOR i=0,nusf-1 DO Print,i,valid(i),usf(i),Format='(I0,I3,2X,A)'
	Message,'One or more invalid files supplied',/Informational
	Return
	ENDIF ELSE BEGIN
	file=StrUpCase(usf)
	cam=IntArr(nusf)
	filt=cam
	FOR i=0,nusf-1 DO BEGIN
	    cfnum=Fix(StrMid(file(i),StrPos(file(i),'.')-2,1))
	    cam(i)=cfnum/5
	    filt(i)=cfnum MOD 5
	    ENDFOR
	IF (Total(cam-cam(0)) NE 0) OR (Total(filt-filt(0)) NE 0) THEN BEGIN
	    Message,'Camera/filter discrepancy!'
	    Return
	    ENDIF ELSE BEGIN
	    cam=cam(0)
	    filt=filt(0)
	    ENDELSE
	ENDELSE
    ENDIF

IF N_Elements(day) EQ 0 THEN day=0

IF (N_Elements(file) EQ 0) THEN BEGIN
   IF ((N_Params() EQ 2) OR (day EQ 0)) THEN BEGIN   ; year, month only => multi-day!
	patt=STRING(year MOD 100,month,fca_suffix $
		   ,FORMAT='(2I2.2,"????",A)')
	REPEAT BEGIN
	    tmp=PickFile(/Read,Path=fcaroot,Filter=patt)
	    IF StrLen(tmp) GT 0 THEN BEGIN
		IF N_Elements(file) EQ 0 THEN BEGIN
		    cfnum=Fix(StrMid(tmp,StrPos(tmp,'.')-2,1))
		    fcam=cfnum/5
		    ffilt=cfnum MOD 5
		    IF (N_Elements(cam) GT 0) AND $
		       (N_Elements(filt) GT 0) THEN BEGIN
			discrep=(cam NE fcam) OR (filt NE ffilt)
			IF discrep THEN Message,'Camera/filter discrepancy!'
			ENDIF ELSE BEGIN
			cam=fcam
			filt=ffilt
			ENDELSE
		    file=StrUpCase(tmp)
		    ENDIF ELSE BEGIN
		    cfnum=Fix(StrMid(tmp,StrPos(tmp,'.')-2,1))
		    tcam=cfnum/5
		    tfilt=cfnum MOD 5
		    IF (tcam NE cam) OR (tfilt NE filt) $
			THEN Message,'Camera/filter discrepancy!'
		    file=[file,StrUpCase(tmp)]
		    ENDELSE
		ENDIF
	    ENDREP UNTIL StrLen(tmp) EQ 0
    
	ENDIF ELSE BEGIN

	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(step) EQ 0 THEN step=interval         ; hours
	IF N_ELEMENTS(finish) EQ 0 THEN finish=24           ; hours UT

	patt=STRING(year MOD 100,month,day,5*cam+filt,num,fca_suffix $
		   ,FORMAT='(3I2.2,I1,I1,A)')

	tmp=PickFile(/Read,Path=fcaroot,Filter=patt)
	IF STRLEN(tmp) GT 0 THEN file=STRUPCASE(tmp) ELSE BEGIN
	    MESSAGE,'No file found; exiting'
	    RETURN
	    ENDELSE
	ENDELSE

    ENDIF

IF N_Elements(file) EQ 0 THEN BEGIN
    Message,'No files selected',/Informational
    Return
    ENDIF ELSE BEGIN
    IF NOT dos THEN file=StrLowCase(file)
    Print,'Files read:'
    FOR i=0,N_Elements(file)-1 DO Print,i,file(i),Format='(I0,2X,A)'
    ENDELSE

IF N_Elements(ps) GT 0 THEN PSOpen,psname,/Long

nfiles=N_Elements(file)
nlin=INTARR(nfiles)
nfov=nlin
year=nlin
month=nlin
day=nlin
jday=LONARR(nfiles)
allfovparms=FLTARR(4,5*nfiles)  ; Data in correct format will fit here
match=1B        ; Flag to test parameter agreement among multiple files

FOR i=0,nfiles-1 DO BEGIN

    ; How big is the file?
    nlin(i)=N_Lines(file(i))

    ; open the file
    OpenR,corrunit,file(i),/GET_LUN

    ; read the number of fields of view used
    READF,corrunit,numfov
    nfov(i)=FIX(ROUND(numfov))
    IF nfov(i) 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(i)+1,5)
    READF,corrunit,fovparms
    allfovparms(0,5*i)=fovparms         ; copy parms for check below
    fovparms=TRANSPOSE(fovparms(1:*,*))

    ; Check whether fovparms agree with first file loaded
    IF i GT 0 THEN BEGIN
	thesematch=allfovparms(*,5*i:5*i+4) EQ allfovparms(*,0:4)
	
	match=match AND Min(thesematch)
	IF NOT match THEN BEGIN
	    MESSAGE,'Parameter mismatch in file '+STRTRIM(i,2) $
		   ,/INFORMATIONAL
	    RETURN
	    ENDIF
	ENDIF

    ; read in the data
    data=FLTARR(nfov(i)+1,nlin(i)-6)
    READF,corrunit,data
    FREE_LUN,corrunit
    
    ; Determine date parameters from file name
    IF dos THEN patt='FCA' ELSE patt='fca'
    inpt=STRPOS(file(i),dd+patt+dd)+5
    IF inpt LT 5 THEN Message,'Failed to decode filename for date'
    year(i)=1900+FIX(STRMID(file(i),inpt,2))
    month(i)=FIX(STRMID(file(i),inpt+2,2))
    day(i)=FIX(STRMID(file(i),inpt+4,2))
    jday(i)=JulDay(month(i),day(i),year(i))

    ; Concatenate data from several days
    IF i EQ 0 THEN BEGIN
	ut=TRANSPOSE(data(0,*))
	afov=TRANSPOSE(data(1:*,*))
	ENDIF ELSE BEGIN
	ut=[ut,TRANSPOSE(data(0,*))+(jday(i)-jday(0))*86400.]
	afov=[afov,TRANSPOSE(data(1:*,*))]
	ENDELSE

    ENDFOR

; Check that data are in correct order
timesort=Sort(ut)
ut=ut(timesort)
FOR j=0,(Size(afov))(2)-1 DO afov(*,j)=afov(timesort,j)

; Check data for discontinuities
ntimes=(Size(ut))(1)
dtime=ut(1:ntimes-1)-ut(0:ntimes-2)
resp=''
Print,String(Max(dtime,loc),ut(loc)/3600. $
	    ,Format='("Max. time gap = ",F7.0," sec at ",F5.1," UT")')
IF (N_Elements(usf) EQ 0) THEN BEGIN
    Read,'Continue with analysis? ([Y]/N) ',resp
    IF StrUpCase(resp) EQ 'N' THEN Return
    ENDIF

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

; identify first and last UT values
mxut=MAX(ut,MIN=mnut)
Print,String(mnut,mxut,mnut/3600.,mxut/3600. $
     ,Format='("Data span ",F7.0," - ",F7.0," s (",F6.1," - ",F6.1,' $
	    +'" h UT)")')

IF N_Elements(start) EQ 0 THEN Read,'Enter start time (hours): ',start
IF N_Elements(interval) EQ 0 $
    THEN Read,'Enter analysis interval (hours): ',interval
IF N_Elements(step) EQ 0 THEN Read,'Enter time step (hours): ',step
IF N_Elements(finish) EQ 0 THEN Read,'Enter finish time (hours): ',finish

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

ncross=nfov(0)
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(0))
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
IF (N_Elements(usf) EQ 0) THEN BEGIN
    PRINT,'Press any key to continue'
    kkk=GET_KBRD(1)
    ENDIF

; 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(0))
FOR i=0,nfov(0)-1 DO iafov(*,i)=interpol(afov(*,i),ut,uut)

; Smooth the data if desired
IF sm GT 0 THEN FOR i=0,nfov(0)-1 DO iafov(*,i)=Smooth(iafov(*,i),sm)

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

; Set up a loop through all time intervals between specified limits

; first=FLOOR(mnut/(interval*3600.)+0.5)>ROUND(start/interval)
first=FLOAT(start)
; last=CEIL(mxut/(interval*3600.)-0.5)<ROUND(finish/interval)
last=FLOAT(finish)
PRINT,first,last $
     ,FORMAT='("Start at ",I0," h UT; finish at ",I0," h UT")'
uthlen=interval
utlen=3600.*uthlen

nsamp=ROUND(uthlen*60.)
IF (nsamp MOD 2) EQ 0 THEN nsamp=nsamp+1

; Now do the loop!
index=-1
FOR hrut=first,last-interval,step DO BEGIN
    good(*)=0B
    badcorr=0B
    index=index+1

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

; Specify the interval length and start time
    uthstart=hrut
    utstart=3600.*uthstart
    PRINT,FORMAT='("---------------------------------------------")'
    PRINT,utstart,utlen,nsamp $
	 ,FORMAT='("Start ",F7.0," s UT; length ",F7.0," s; ",I0' $
		+'," samples")'

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

; Check for shallow fading
    nofade=0B
    FOR i=0,nfov(0)-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

; 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 the interval lies at either end of the data set, there won't be
; enough data points for lags in one direction.  Settle for lags up to
; one quarter of the total interval length, and redefine the unshifted
; array to work with the reduced time.

;            IF nw EQ nwt THEN BEGIN

	    IF nwt GT 0.5*nw THEN BEGIN

; The following section is new.  DPS, 960516

		IF nwt LT nw THEN BEGIN
		    ndiff2=(nw-nwt)/2
		    nutstart=uut(w(0))+ndiff2*60
		    mw=WHERE((nutstart LE uut) AND $
			     (uut LT (nutstart+nwt*60)),mwt)
;                   HELP,nw,nwt,ndiff2,nutstart,utstart,mwt
;                   STOP
		    ENDIF ELSE BEGIN
		    mw=w
		    wt=wt(0:nw-1)
		    ENDELSE

		IF N_ELEMENTS(mw) NE N_ELEMENTS(wt) THEN STOP
		    
		FOR i=0,nfov(0)-1 DO $
		    auto(i,j)=correlate(iafov(mw,i),iafov(wt,i))
		FOR i=0,ncross-1 DO $
		    cross(i,j)=correlate(iafov(mw,p(i)),iafov(wt,q(i)))
		ENDIF ELSE BEGIN
		MESSAGE,'Not enough data for ' $
		       +STRING(tau(j),FORMAT='(F7.2)')+' minute lag ' $
		       ,/INFORMATIONAL
		ENDELSE
	    ENDFOR
;           hlp,auto,cross
	    badcorr=(TOTAL(auto) EQ 0) AND (TOTAL(cross) EQ 0)

; 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') OR (!D.NAME EQ 'WIN') THEN BEGIN
		ss=Screen_Size()
		goldrat=(1+Sqrt(5))/2
		WINDOW,0,XSIZE=ss(1)/goldrat,YSIZE=ss(1)
		ENDIF ELSE ERASE
	    jdstart=JulDay(month(0),day(0),year(0))+uthstart/24
	    tstring=STRING(JD2Date(jdstart,Form="d$ n$ y$") $
			  ,filters(5*cam+filt),uthstart,uthlen $
			  ,FORMAT='(A,"  ",A," FCA: start "' $
				 +',F5.2,", len ",F5.2)')
	    IF sm GT 0 $
		THEN tstring=tstring $
			    +String(sm,Format='("  (",I0," min. smoothing)")')

	    PLOT,tau,auto(0,*),XTITLE='Time lag (minutes)' $
			      ,YTITLE='Autocorrelation' $
			      ,YRANGE=[0,1.0],/YSTYLE,TITLE=tstring $
			      ,LINESTYLE=1
	    FOR i=1,nfov(0)-1 DO OPLOT,tau,auto(i,*),LINE=i+1
	    mauto=REFORM(TOTAL(auto,1)/nfov(0))
	    OPLOT,tau,mauto

	    wgood=WHERE(mauto NE 0.)

;           STOP
	    
	    symnum=1+2*(nsamp GT 61)
	    autofit=gaussfnb(tau(wgood),mauto(wgood),fp)
	    OPLOT,tau(wgood),autofit,PSYM=symnum

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

	    PLOT,tau,mauto,PSYM=symnum $
		,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))
	    yspace=sizech(1,1)-sizech(1,0)

; Plot and label the Gaussian fit to the mean autocorrelation function
	    zf=((xt-fp(1))/fp(2))<10>(-10)
	    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: ",F7.2)'),ALIGNMENT=0
	    XYOUTS,xl+2,1.-3*yspace $
		  ,STRING(fp(2),FORMAT='("Sigma:  ",F7.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)")')

;               OPLOT,tau,SMOOTH(cross(i,*),3)
		cross(i,*)=SMOOTH(cross(i,*),3)
;               cross(i,*)=CONVOL(cross(i,*),[[0.25],[0.5],[0.25]] $
;                                ,/EDGE_TRUNCATE)

		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
	    timestmp

	    IF TOTAL(good) EQ ncross THEN BEGIN
		PRINT,'Lag at peak cross-correlation: ',pktau
		PRINT,'Value of pk cross-correlation: ',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
		    over(0)=-1          ; force -1 at
		    over(nsamp-1)=-1    ; array end-points
		    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: ",2F7.2," min.")'
		    PRINT,tauarg1 $
			 ,FORMAT='("Mean 60.7% correlation lag: ",F7.2," min.")'
		    FOR i=0,1 DO PLOTS,[tau61(i),tau61(i)],[4.8,5.8]
		    input(*,index)=[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(*,index)=[utstart,utlen,a,b,tilt $
					       ,bigV,dir,c/60.]

				ENDELSE
			    ENDELSE
			ENDELSE
		    ENDIF
		ENDIF
	    ENDIF
	ENDELSE
    IF Keyword_Set(wait) THEN kkk=Get_Kbrd(1)
    ENDFOR

; Now save the analysis input data to a file
datfile=STRMID(file(0),0,STRLEN(file(0))-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),0,STRLEN(file(0))-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
IF Keyword_Set(ps) THEN PSClose

Return
END