;+
; NAME:
;	IMF_CONF
; PURPOSE:
;	Cross-correlation of IMF By and polar cap convection direction, 
;	December 8, 1993.
; CATEGORY:
;	Event specific data analysis
; CALLING SEQUENCE:
;	IMF_CONV, IMF_FILE, IMF_COL, CONV_FILE [, CONV_COL] $
;               , START=start, FINISH=finish, PS=ps, NPKS=npks $
;	        , TRIPLEX=triplex, SFIT=sfit, CORRSIGN=corrsign
; INPUTS:
;	IMF_FILE: Specifies the full pathname of the file containing
;		the IMF data (in either Cartesian or polar GSM coordinates).
;	IMF_COL: Specifies the column of the IMF_FILE containing the 
;		coordinate to be used for the cross-correlation.
;	CONV_FILE: Specifies the full pathname of the file containing
;		the convection parameters to be used.
; OPTIONAL INPUTS:
;	CONV_COL: If set, specifies the column of the CONV_FILE containing
;		the parameter to be used for the cross-correlation (D = 10).
; KEYWORD PARAMETERS:
;	/HELP:	If set, displays internal documentation.
;	START:	If specified, gives the UT start hour of the correlation 
;		interval.
;	FINISH: As above, except the finish hour.
;	MAXLAG:	If set, specifies the maximum lag between the IMF By and the 
;		convection direction for which a cross-correlation is to be 
;		computed.  The default is (finish - start).
;	/PS:	If set, send the plots to a PS file instead of a window.
;	NPKS:	Specifies the number of cross-correlation peaks to annotate.
;	/TRIPLEX: If set, all text is in Triplex Roman instead of Simplex Roman.
;	/SFIT:	If set, do separate fits of convection direction to IMF By for 
;		By > 0 and By < 0.
;	CORRSIGN: Determines whether positive correlations (CORRSIGN=1) or
;		negative correlations (CORRSIGN=-1) are sought.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	Two files are read.  If /PS is set, a PostScript plot file is created.
; RESTRICTIONS:
;	Hardwired to use two specific file names.  Could be generalized with
;	moderate effort.
; PROCEDURE:
;	Straightforward.  Data sets are interpolated to 1-min steps.
; EXAMPLE:
;	
; SEE ALSO:
;	
; MODIFICATION HISTORY:
; 	Written by:	DPS, ISR/ISAS, Sept. 9-18, 1995.
;-
PRO imf_conv,imf_file,imf_col,conv_file,conv_col $
            ,help=help,start=start,finish=finish,maxlag=maxlag,ps=ps,npks=npks $
            ,triplex=triplex,sfit=sfit,corrsign=corrsign

IF KEYWORD_SET(help) THEN BEGIN
    doc_library,'imf_conv'
    RETURN
    ENDIF

IF N_PARAMS() EQ 3 THEN conv_col=10
IF N_ELEMENTS(start) EQ 0 THEN start=0.
IF N_ELEMENTS(finish) EQ 0 THEN finish=24.
IF N_ELEMENTS(maxlag) EQ 0 THEN maxlag=finish-start
IF N_ELEMENTS(npks) EQ 0 THEN npks=3
IF KEYWORD_SET(triplex) THEN font='!17' ELSE font='!3'
IF N_ELEMENTS(sfit) EQ 0 THEN sfit=0
IF N_ELEMENTS(corrsign) EQ 0 THEN corrsign=1
IF ABS(corrsign) NE 1 THEN corrsign=1

; Read in the IMF data file

; imf=TRANSPOSE(ddread('\images\pgm\dec0893.txt'))
imf=TRANSPOSE(ddread(imf_file))
utimf=imf(*,0)
by=imf(*,imf_col)
butrange=[MIN(utimf,MAX=tmp),tmp]

; Read in the convection velocity file

; conv=qgetfile('\images\pgm\93120802.cmg')
conv=qgetfile(conv_file)
nl=N_ELEMENTS(conv)
lasttext=WHERE(STRPOS(conv,'IM#') EQ 0,nlt)
IF nlt LT 1 THEN BEGIN
    MESSAGE,'Couldn''t find end of text',/INFORMATIONAL
    RETURN
    ENDIF
ndata=nl-lasttext(0)-1
data=INDGEN(ndata)+lasttext(0)+1

; Set up convection parameter arrays

keep=conv(data)
utconv=FLTARR(ndata)
dev=utconv

; Extract convection parameters

len=STRLEN(keep)
maxlen=MAX(len)
FOR i=0,ndata-1 DO BEGIN
    ki=keep(i)
    utconv(i)=FLOAT(getwrd(ki,2))
    IF len(i) EQ maxlen THEN BEGIN
	dev(i)=FLOAT(getwrd(ki,conv_col))
	ENDIF
    ENDFOR

; Interpolate missing data

dutrange=[MIN(utconv,MAX=tmp),tmp]
full=WHERE(len EQ maxlen,nfull)
omit=WHERE(len LT maxlen,nomit)
dev(omit)=interpol(dev(full),utconv(full),utconv(omit))

; Extend both data sets artificially to a longer time span

datarange=[butrange(0)<dutrange(0),butrange(1)>dutrange(1)]
fillrange=datarange+[-10000.,+10000.]
dutimf=MEDIAN(utimf-SHIFT(utimf,1))
dutconv=MEDIAN(utconv-SHIFT(utconv,1))
utimf=[fillrange(0),utimf(0)-dutimf $
      ,utimf $
      ,utimf(N_ELEMENTS(utimf)-1)+dutimf,fillrange(1)]
by=[0.,0.,by,0.,0.]
utconv=[fillrange(0),utconv(0)-dutconv $
      ,utconv $
      ,utconv(N_ELEMENTS(utconv)-1)+dutconv,fillrange(1)]
dev=[0.,0.,dev,0.,0.]

; Interpolate both data sets to a uniform time grid

urange=ROUND(fillrange/60.)*60
utime=makex(FLOAT(urange(0)),FLOAT(urange(1)),60.)
uby=interpol(by,utimf,utime)
udev=interpol(dev,utconv,utime)

; Locate the interval over which to do the cross-correlation

first=3600.*start>butrange(0)>dutrange(0)
last=3600.*finish<butrange(1)<dutrange(1)
w=WHERE((first LE utime) AND (utime LE last))
uby=uby(w)

; Now do the cross-correlations

ncc=60*maxlag+1
crosscorr=FLTARR(ncc)
FOR i=0,ncc-1 DO BEGIN
    shudev=udev(w+i)
    crosscorr(i)=correlate(uby,shudev)
    ENDFOR

; Locate the peak cross-correlations and sort by magnitude

peaks=extremes(crosscorr,corrsign)
peaks=peaks(WHERE(crosscorr(peaks)/corrsign GT 0.,npeaks))
npks=npks<npeaks
ord=reverse(SORT(crosscorr(peaks)/corrsign))
top3=peaks(ord(0:npks-1))

; Plot the results

IF KEYWORD_SET(ps) THEN psopen,'imf_conv.ps',/long $
		   ELSE WINDOW,0,XSIZE=650,YSIZE=650
!P.MULTI=[0,1,2,0,0]
lag=FINDGEN(ncc)
PLOT,lag,crosscorr $
    ,XTITLE='Lag (minutes) between IMF B!iy!n and convection direction' $
    ,YTITLE='Cross-correlation coefficient' $
    ,YRANGE=[-1,1],/YSTYLE $
    ,TITLE=font+'Cross-Correlation of IMF B!iy!n and ' $
          +'Polar Cap Convection Direction' $
    ,LINESTYLE=0
PLOTS,!X.CRANGE,REPLICATE(0.0,2),LINESTYLE=1
charht=CONVERT_COORD([0,!D.X_CH_SIZE],[0,!D.Y_CH_SIZE],/DEVICE,/TO_DATA)
charht=charht(1,1)-charht(1,0)
; XYOUTS,0.75*!X.CRANGE(0)+0.25*!X.CRANGE(1),!Y.CRANGE(0)+2*charht $
;       ,'IMF B!iy!n lags convection direction',ALIGNMENT=0.5
; XYOUTS,0.25*!X.CRANGE(0)+0.75*!X.CRANGE(1),!Y.CRANGE(0)+2*charht $
;       ,'IMF B!iy!n leads convection direction',ALIGNMENT=0.5
XYOUTS,0.05*!X.CRANGE(0)+0.95*!X.CRANGE(1),!Y.CRANGE(0)+charht $
      ,'December 8, 1993, '+STRMID(strsec(first),0,5) $
       +' - '+STRMID(strsec(last),0,5)+' UT' $
      ,ALIGNMENT=1
FOR i=0,npks-1 DO BEGIN
    lagi=lag(top3(i))
    cci=crosscorr(top3(i))
    delht=(((!Y.CRANGE(1)-cci)-2.*charht)/3.)<(0.3*charht)
    HELP,delht
    PLOTS,REPLICATE(lagi,2),[!Y.CRANGE(0),cci]
    XYOUTS,lagi,cci+2*delht+charht $
	  ,STRING(top3(i),FORMAT='("!7s!x = ",I0," min")') $
	  ,ALIGNMENT=0.5
    XYOUTS,lagi,cci+delht $
	  ,STRING(cci,FORMAT='("r=",F5.2)') $
	  +STRING(ROUND(cci^2*100),FORMAT='(" (",I2,"%)")') $
	  ,ALIGNMENT=0.5
    ENDFOR

; Evaluate the linear regression coefficients between IMF By and drift angle

bestlagdev=udev(w+top3(0))
lagtime=utime(w+top3(0))
lagst=STRING(top3(0),FORMAT='("(!7s!x = ",I0," min)")')
PLOT,lagtime/3600.,bestlagdev $
    ,XTITLE='UT (hours), December 8, 1993' $
    ,YTITLE='Drift Angle (deg CCW from antisunward)' $
    ,YRANGE=[-30,50],/YSTYLE $
    ,TITLE='IMF B!iy!n Contribution to Drift Angle Variation '+lagst

; Get ready to show a legend

xdvpos=!P.CLIP(0)
ydvpos=!P.CLIP(3)
npos=CONVERT_COORD([xdvpos,ydvpos],/DEVICE,/TO_NORMAL)
npos=npos(0:1)
lform='("DA = ",F4.1," deg + ",F4.1," deg/nT")'
IF KEYWORD_SET(sfit) THEN BEGIN
    pfit=0.
    nfit=0.
    pby=WHERE(uby GE 0.,npby)
    nby=WHERE(uby LT 0.,nnby)
    pparms=svdfit(uby(pby),bestlagdev(pby),2,yfit=pfit)
    nparms=svdfit(uby(nby),bestlagdev(nby),2,yfit=nfit)
    PRINT,pparms
    PRINT,nparms
    OPLOT,lagtime(pby)/3600.,pfit,PSYM=1
    OPLOT,lagtime(nby)/3600.,nfit,PSYM=4
    labels=['B!iy!n positive;  '+STRING(pparms,FORMAT=lform) $
           ,'B!iy!n negative; '+STRING(nparms,FORMAT=lform)]
    syms=[1,4]
    ENDIF ELSE BEGIN
    bfit=0.
    bparms=svdfit(uby,bestlagdev,2,yfit=bfit)
    PRINT,bparms
    OPLOT,lagtime/3600.,bfit,PSYM=2
    labels='Best fit: '+STRING(bparms,FORMAT=lform)
    syms=2
    ENDELSE

legend,labels,psym=syms,position=npos

IF KEYWORD_SET(ps) THEN psclose

RETURN
END