PRO fcacompare,date,mtg_type=mtg_type,colour=colour

IF N_PARAMS() EQ 0 THEN BEGIN
    MESSAGE,'Usage: fcacompare,date,mtg_type=mtg_type,colour=colour' $
	   ,/INFORMATIONAL
    RETURN
    ENDIF

IF KEYWORD_SET(colour) THEN BEGIN
    lstyles=REPLICATE(0,4)
    cols=[0.4,0.6,0.8,0.9]
    loadct,13
    whitetop
    ENDIF ELSE BEGIN
    lstyles=[1,2,3,4]
    cols=REPLICATE(1,4)
    loadct,0
    ; whitebot
    ; blacktop
    ENDELSE

yscl=1.5
eurlat=80.0533
eurlon=273.5839
; Convert the input LONG date to a string for file-name matching
ds=STRING(date,FORMAT='(I6.0)')
yy=date/10000L
mm=(date-10000L*yy)/100L
dd=date MOD 100
filepatt='/jasper/cnsr3_data1/fca/'+ds+'*.fca'
files=FINDFILE(filepatt,COUNT=nfiles)
IF nfiles EQ 0 THEN BEGIN
    MESSAGE,'No files found for '+ds,/INFORMATIONAL
    RETURN
    ENDIF ELSE BEGIN
    parms=FLTARR(8,48,nfiles)
    parm1=FLTARR(8,48)
    FOR i=0,nfiles-1 DO BEGIN
	OPENR,fcaunit,files(i),/GET_LUN
	READF,fcaunit,parm1
	FREE_LUN,fcaunit
	parms(0,0,i)=parm1
	ENDFOR

    IF !D.NAME EQ 'PS' THEN psres=100 ELSE WINDOW,0,XSIZE=904,YSIZE=660
    multiplot,[8,6],/init,/verbose
    xleft=!X.MARGIN(0)*!D.X_CH_SIZE
    xright=!D.X_SIZE-!X.MARGIN(1)*!D.X_CH_SIZE
    xcenter=0.5*(xleft+xright)
    ybottom=!Y.MARGIN(0)*!D.Y_CH_SIZE
    ytop=!D.Y_SIZE-!Y.MARGIN(1)*!D.Y_CH_SIZE
    ycenter=0.5*(ybottom+ytop)
    XYOUTS,xcenter,ytop+0.5*!D.Y_CH_SIZE $
	  ,'F-LAYER PATCHES  FULL CORRELATION ANALYSIS  '+ds $
	  ,CHARSIZE=1.8,ALIGNMENT=0.5,/DEVICE
    XYOUTS,xcenter,ybottom-3.5*!D.Y_CH_SIZE $
	  ,'    TOP 3 ROWS: Dusk-dawn drift velocity (m/s)   ' $
	  +'       BOTTOM 3 ROWS: Dusk-dawn patch scale (km)    ' $
	  ,/DEVICE,ALIGNMENT=0.5
    XYOUTS,xleft-6.5*!D.X_CH_SIZE,ycenter $
	  ,'Noon-midnight patch scale (km)    ' $
	  +'        Antisunward drift velocity (m/s)' $
	  ,/DEVICE,ALIGNMENT=0.5,ORIENTATION=90.

    ; First compare velocities
    ;
    vxfca=FLTARR(4)
    vyfca=vxfca

    ; Read in drift velocities derived from individual patch
    ; analysis, and average over half-hour intervals
    ; The following code was largely chopped out of mkmtg.pro.
    ; ---------------------------------------------------------------------
    IF N_ELEMENTS(mtg_type) EQ 0 THEN BEGIN
	mtg_type='g' 
	MESSAGE,'Default montage type is Geographic',/INFORMATIONAL
	ENDIF ELSE BEGIN
	CASE STRLOWCASE(mtg_type) OF
       'g': mtg_type=STRLOWCASE(mtg_type)
       'p': mtg_type=STRLOWCASE(mtg_type)
      ELSE: REPEAT BEGIN
		READ,'Enter a montage type (Geographic/Pace): ',mtg_type
		mtg_type=STRLOWCASE(mtg_type)
		ENDREP UNTIL (mtg_type EQ 'g') OR (mtg_type EQ 'p')
	    ENDCASE
	ENDELSE
    path='/jasper/cnsr3_data1/pgm/'+ds+'02.mt'+mtg_type
    mtg=qgetfile(path)
    n=N_ELEMENTS(mtg)
    filine=0
    WHILE STRMID(mtg(filine),0,3) NE 'IM#' DO filine=filine+1
    READS,mtg(filine+1),firstn,FORMAT='(I3)'
    READS,mtg(n-1),lastn,FORMAT='(I3)'
    npix2=lastn-firstn+1
    PRINT,npix2,' images analysed',FORMAT='(I0,A)'
    info0=n-npix2
    times=fltarr(npix2)
    vdusk=times
    vasun=times
    p1=''
    FOR i=0,npix2-1 DO BEGIN
        j=i+info0
        READS,mtg(j),p0,p1,p2,p3,p4,p5,p6,p7,p8 $
             ,FORMAT='(I3,2X,A12,2X,F8.2,2X,I3,2X,F8.4,2(2X,I3),2(2X,F8.5))'
        times(i)=p2
        vasun(i)=p7
        vdusk(i)=p8
        ENDFOR
    novels=BYTARR(npix2)
    gap=WHERE((vasun EQ 0.) and (vdusk EQ 0.),ngap)
    IF ngap GT 0 THEN BEGIN
        novels(gap)=1
        good=REPLICATE(1B,N_ELEMENTS(novels))-novels
        FOR i=0,ngap-1 DO BEGIN
	    tmp=extrac(vasun,gap(i)-3,7)
	    vasun(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
	    tmp=extrac(vdusk,gap(i)-3,7)
	    vdusk(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
            ENDFOR
        ENDIF
    vasun=vasun*111660.		; convert from deg/s to m/s;
    vdusk=vdusk*111660.		; assume 111.66 km per deg of latitude
    
    ymd=date
    dst=ymd2date(1900+ymd/10000,(ymd/100) MOD 100,ymd MOD 100 $
	        ,format='d$ n$ Y$')
    PRINT,dst
    mint=MIN(times/3600.,MAX=maxt)
    PRINT,'Patch velocities span ',mint,' to ',maxt,' hours UT' $
	 ,FORMAT='(2(A,F5.2),A)'

    vavg=FLTARR(5,24)		; # samples, <V antisun>, <V 2dusk>,
				; s.d. of Vas, s.d. of V2d
    first=FLOOR(MIN(times/1800.))	; first half-hour with vels
    last=FLOOR(MAX(times/1800.))	; last half-hour with vels
    FOR j=first,last DO BEGIN
	win=WHERE((j LE times/1800.) AND (times/1800. LT j+1),nwin)
	IF nwin GT 1 THEN BEGIN
	    vavg(0,j)=nwin
	    vavg(3,j)=stdev(vasun(win),mvasun)
	    vavg(1,j)=mvasun
	    vavg(4,j)=stdev(vdusk(win),mvdusk)
	    vavg(2,j)=mvdusk
	    ENDIF
	ENDFOR
    ; ---------------------------------------------------------------------

    flag=0B
    FOR i=0,23 DO BEGIN
	doit=(3 LE i) AND (i LE 20)
	multiplot
	set_isoxy,-1000,1000,2000,0
        PLOT,[0,1],/NODATA $
	    ,XRANGE=[!X.CRANGE(1),!X.CRANGE(0)],/XSTYLE,XTICKS=4 $
	    ,XTICKV=[-500,0.0,500],YTICKV=[0.0,500,1000,1500] $
	    ,YRANGE=[!Y.CRANGE(0),!Y.CRANGE(1)],/YSTYLE,YTICKS=4 $
	    ,CHARSIZE=0.75
	FOR j=0,nfiles-1 DO BEGIN
	    IF TOTAL(parms(*,i,j)) GT 0 THEN BEGIN
		p=parms(*,i,j)
		IF NOT flag THEN BEGIN
		    utmid=p(0)+p(1)/2
		    h=FLOOR(utmid/3600)
		    m=FLOOR((utmid-3600*h)/60)
		    s=utmid-(h*60+m)*60
		    sunazel,yy,mm,dd,h,m,s,saz
		    IF doit THEN XYOUTS,0,1900. $
			  ,STRMID(strsec(p(0)),0,5)+'-' $
			  +STRMID(strsec(p(0)+p(1)),0,5) $
			  ,ALIGNMENT=0.5,CHARSIZE=0.75
		    ENDIF
		angle=p(6)-(saz-180.)			; + if drift duskward
		arrx=-p(5)*SIN(angle*!DTOR)*1000.	; + if drift dawnward
		arry=p(5)*COS(angle*!DTOR)*1000.
		vxfca(j)=arrx
		vyfca(j)=arry
		IF doit THEN myarrow,0.0,0.0,arrx,arry,/data $
		       ,hsize=-0.05 $;,col=cols(j)*(!D.N_COLORS-1) $
		       ,linestyle=lstyles(j)
		flag=flag+1
                ENDIF ELSE BEGIN
		vxfca(j)=0.
		vyfca(j)=0.
		ENDELSE
            ENDFOR

	IF flag GE 3 THEN BEGIN		; show the mean FCA velocity
	    wv=WHERE(vxfca*vxfca+vyfca*vyfca)
	    sdvx=stdev(vxfca(wv),mvx)
	    sdvy=stdev(vyfca(wv),mvy)
	    IF doit THEN ellipse,sdvx,sdvy,0.,0.,360.,mvx,mvy $
		   ,/data $;,col=!D.N_COLORS-2 $
		   ,linestyle=0
	    ENDIF

	IF vavg(0,i) GT 0 THEN BEGIN	; show mean patch velocity
	    IF doit THEN myarrow,0.0,0.0,-vavg(2,i),vavg(1,i),/data $
		   ,hsize=-0.1 $;,col=!D.N_COLORS-1 $
		   ,linestyle=0
	    IF doit THEN ellipse,vavg(4,i),vavg(3,i),0.,0.,360.,-vavg(2,i),vavg(1,i) $
		   ,/data $;,col=!D.N_COLORS-1 $
		   ,linestyle=0
	    ENDIF
	XYOUTS,900.,300.,STRING(FIX(flag),FORMAT='(I1)') $
	      ,ALIGNMENT=1 $
	      ,CHARSIZE=0.75
	flag=0B
        ENDFOR

    ; Now compare patch shapes
    flag=0B
    FOR i=0,23 DO BEGIN
	doit=(3 LE i) AND (i LE 20)
	multiplot
        set_isoxy,-1000,1000,-1000,1000
	PLOT,[0,1],/NODATA $
	    ,XTICKS=4,XTICKV=[-500,500] $
	    ,YTICKS=4,YTICKV=[-500,500] $
	    ,CHARSIZE=0.75
	FOR j=0,nfiles-1 DO BEGIN
	    IF TOTAL(parms(*,i,j)) GT 0 THEN BEGIN
		p=parms(*,i,j)
		IF NOT FLAG THEN BEGIN
		    IF doit THEN XYOUTS,0,-900 $
			  ,STRMID(strsec(p(0)),0,5)+'-' $
			  +STRMID(strsec(p(0)+p(1)),0,5) $
			  ,ALIGNMENT=0.5 $
			  ,CHARSIZE=0.75
		    utmid=p(0)+p(1)/2
		    h=FLOOR(utmid/3600)
		    m=FLOOR((utmid-3600*h)/60)
		    s=utmid-(h*60+m)*60
		    sunazel,yy,mm,dd,h,m,s,saz
		    ENDIF
		IF doit THEN ellipse,p(2),p(3),90.-(p(4)-saz),0.,360.,0.,0. $
		        $;,col=cols(j)*(!D.N_COLORS-1) $
		       ,linestyle=lstyles(j)
		flag=flag+1
                ENDIF
            ENDFOR
	IF doit THEN XYOUTS,900,700,STRING(FIX(flag),FORMAT='(I1)') $
	      ,ALIGNMENT=1 $
	      ,CHARSIZE=0.75
	flag=0B
        ENDFOR

    multiplot,[1,1],/init,/verbose
    ENDELSE

RETURN
END