PRO fcabriggs,noplot=noplot file='fca/93120802.fca' nlin=nlines(file) ; open the file OPENR,fcaunit,file,/GET_LUN ; read the number of fields of view used READF,fcaunit,nfov ; read in the field of view parameters fovparms=FLTARR(nfov+1,5) READF,fcaunit,fovparms fovparms=TRANSPOSE(fovparms(1:*,*)) ; read in the data data=FLTARR(nfov+1,nlin-6) READF,fcaunit,data ut=TRANSPOSE(data(0,*)) afov=TRANSPOSE(data(1:*,*)) year=1993 month=12 day=8 action=[' - accept data',' - reject data'] xx=REFORM(fovparms(3,*)) yy=REFORM(fovparms(4,*)) ; Note the difference between nfov and C(nfov,2)!!!! ; ncross=ROUND(GAMMA(nfov+1) $ ; no. of cross-correlations ; /(GAMMA(2+1)*GAMMA(nfov-2+1))) ncross=nfov ; since FOV pairs are to be taken in sequence xi=FLTARR(ncross) eta=xi good=BYTARR(ncross) pktau=FLTARR(ncross) ; Now identify the FOV pairs involved in the cross-correlations p=INDGEN(ncross) q=p ; k=0 ; FOR i=0,nfov-2 DO BEGIN ; FOR j=i+1,nfov-1 DO BEGIN ; p(k)=i ; q(k)=j ; k=k+1 ; ENDFOR ; ENDFOR p=INDGEN(nfov) q=(p+1) MOD nfov ; locate a UT gap and restrict data to before gap dut=(SHIFT(ut,-1)-ut)>0 gap=WHERE(dut GT 500,ngap) ut=ut(0:gap(0)) afov=afov(0:gap(0),*) ; 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) ; Branch back to here for repeated intervals DO_AGAIN: good(*)=0B ; Now identify the interval length and start time READ,'Enter the start time and interval in (fractional) UT hours: ' $ ,uthstart,uthlen utstart=3600.*uthstart nsamp=ROUND(uthlen*60.) IF (nsamp MOD 2) EQ 0 THEN nsamp=nsamp+1 utlen=3600.*uthlen PRINT,utstart,utlen,nsamp ; now identify data from the first full half-hour period 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) ; Now 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 MESSAGE,'Not enough data for '+ $ STRING(tau(j),FORMAT='(F6.2)')+ $ ' minute lag ' $ ,/INFORMATIONAL ENDFOR ; Now plot out the auto- and cross-correlations ; First the auto-correlations IF NOT KEYWORD_SET(noplot) THEN BEGIN !P.MULTI=[0,1,ncross+2,0,0] IF !D.NAME EQ 'X' THEN WINDOW,0,XSIZE=640,YSIZE=900 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 FOR i=1,nfov-1 DO OPLOT,tau,auto(i,*),LINE=i autofit=gaussfit(tau,REFORM(TOTAL(auto,1)/nfov),fp) OPLOT,tau,autofit,PSYM=1 w0=(WHERE(tau EQ 0))(0) PLOT,tau,TOTAL(auto,1)/nfov,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 corn=[[0,0],[!D.X_CH_SIZE,!D.Y_CH_SIZE]] szch=CONVERT_COORD(corn,/DEVICE,/TO_DATA) yspc=1.25*(szch(1,1)-szch(1,0)) zf=(xt-fp(1))/fp(2) OPLOT,xt,poly(xt,fp(3:5)) OPLOT,xt,fp(0)*EXP(-zf^2/2)+poly(xt,fp(3:5)),THICK=2 XYOUTS,xl+2,1.-yspc,STRING(fp(0),FORMAT='("Height: ",F6.2)'),ALIGNMENT=0 XYOUTS,xl+2,1.-2*yspc,STRING(fp(1),FORMAT='("Center: ",F6.2)'),ALIGNMENT=0 XYOUTS,xl+2,1.-3*yspc,STRING(fp(2),FORMAT='("Sigma: ",F6.2)'),ALIGNMENT=0 XYOUTS,xr-2,1.-yspc,STRING(fp(3),FORMAT='("Constant: ",F7.4)'),ALIGNMENT=1 XYOUTS,xr-2,1.-2*yspc,STRING(fp(4),FORMAT='("Linear: ",F7.4)'),ALIGNMENT=1 XYOUTS,xr-2,1.-3*yspc,STRING(fp(5),FORMAT='("Quadratic: ",F7.4)'),ALIGNMENT=1 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='(" 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 minimum 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 ; Fit a parabola to the peak pkfit=svdfit(tau(wpk),REFORM(cross(i,wpk)),3,yfit=yf) pkpq=-pkfit(1)/(2.*pkfit(2)) ; tau at peak PRINT,'Peak cross-correlation at tau = ',STRTRIM(pkpq,2),' min.' OPLOT,tau(wpk),yf,PSYM=1 algn=(wpk(2) LT nsamp/2) xv=xl+2+algn*(xr-xl-4) XYOUTS,xv,1.-0.5*yspc $ ,STRING(poly(pkpq,pkfit),pkpq $ ,FORMAT='("Peak ",F5.3," at ",F7.3," min")') $ ,ALIGNMENT=algn good(i)=1B pktau(i)=pkpq xi(i)=xx(q(i))-xx(p(i)) eta(i)=yy(q(i))-yy(p(i)) ENDELSE ENDELSE ENDFOR IF TOTAL(good) EQ ncross THEN BEGIN ntd=ABS(TOTAL(pktau))/TOTAL(ABS(pktau)) PRINT,STRING(ntd,action(ntd GE 0.2) $ ,FORMAT='("Normalized time discrepancy = ",F5.3,A)') ; Try to fit for -F/C and -G/C using eqn. (22) of Briggs, MAP handbook. ; Start by accumulating needed sums. s0=TOTAL(pktau*xi) s1=TOTAL(xi*xi) s2=TOTAL(xi*eta) s3=TOTAL(pktau*eta) s4=TOTAL(eta*eta) ; Now compute fit parameters: p1 = F/C, p2 = G/C. p1=(s3*s2-s0*s4)/(s1*s4-s2*s2) p2=(s0*s2-s1*s3)/(s1*s4-s2*s2) PRINT,p1,p2 $ ,FORMAT='("Results of fit: F/C=",E10.3,", G/C=",E10.3)' PRINT,'Measured tau Fitted tau' PRINT,[TRANSPOSE(pktau),-TRANSPOSE(p1*xi+p2*eta)] alt=-regress([TRANSPOSE(xi),TRANSPOSE(eta)] $ ,pktau,REPLICATE(1.,ncross),/relative_weight $ ,altfit,const,altsig,ftest,r,rmul,chisq) fbyc=[alt(0),altsig(0)] gbyc=[alt(1),altsig(1)] PRINT,'Fitting for F/C, G/C:' PRINT,alt,FORMAT='("Fitted coefficients: ",2F10.6)' PRINT,altsig,FORMAT='("Standard deviations: ",2F10.7)' PRINT,r,FORMAT='("Linear correlations: ",2F10.6)' PRINT,rmul,FORMAT='("Multiple linear correlation: ",F10.6)' PRINT,const,FORMAT='("Constant term in fit: ",F10.6)' ; Now compute A/C, B/C, and H/C using eqn. (23) of Briggs, MAP Handbook. alt=regress([TRANSPOSE(xi*xi),TRANSPOSE(eta*eta),TRANSPOSE(xi*eta)] $ ,pktau*pktau,REPLICATE(1.,ncross),/relative_weight $ ,altfit,const,altsig,ftest,r,rmul,chisq) abyc=[alt(0),altsig(0)] bbyc=[alt(1),altsig(1)] hbyc=0.5*[alt(2),altsig(2)] PRINT,'Fitting for A/C, B/C, and H/C:' PRINT,alt,FORMAT='("Fitted coefficients: ",3F10.6)' PRINT,altsig,FORMAT='("Standard deviations: ",3F10.7)' PRINT,r,FORMAT='("Linear correlations: ",3F10.6)' PRINT,rmul,FORMAT='("Multiple linear correlation: ",F10.6)' PRINT,const,FORMAT='("Constant term in fit: ",F10.6)' WSHOW,0,/ICONIC ENDIF ENDIF ENDELSE ; Do it again? resp='' READ,'Repeat for another time interval? ([Y]/N) ',resp IF STRUPCASE(resp) NE 'N' THEN GOTO,DO_AGAIN END