PRO mkmtg,name,ps=ps,fillgaps=fillgaps,debug=debug

; IDL Version 3.1.0 (RISC/os mipseb)
; Journal File for cnsr3@jasper
; Working directory: /jasper/cnsr3_data1
; Date: Fri Jun 17 08:21:03 1994

IF N_PARAMS() EQ 0 THEN BEGIN
    PRINT,'Usage: mkmtg,name,ps=ps,fillgaps=fillgaps,debug=debug'
    RETURN
    ENDIF
 
path='/jasper/cnsr3_data1/pgm/'+name
type=STRMID(name,STRLEN(name)-1,1)	; g for geographic, p for PACE
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)
colshift=INTARR(npix2)
rowshift=colshift
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
    colshift(i)=p5
    rowshift(i)=p6
    vasun(i)=p7
    vdusk(i)=p8
    ENDFOR

dateline=0
WHILE (isnumber(getwrd(mtg(dateline),0)) EQ 0) DO dateline=dateline+1
date2ymd,mtg(dateline),year,month,day
IF year LT 50 THEN year=year+2000 ELSE year=year+1900

htline=dateline
WHILE (getwrd(mtg(htline),0) NE 'Assumed') DO htline=htline+1
emis_ht=FLOAT(getwrd(mtg(htline),3))

pixres=FLOAT(getwrd(mtg(info0-3),7))

good=WHERE((colshift NE 0) OR (rowshift NE 0),ngood)	; meaningful shifts

IF STRUPCASE(type) EQ 'P' THEN BEGIN

    drift_az=FLTARR(npix2)
    daz=drift_az
    v=drift_az

    ; Following lines added 941110 by D P Steele to get actual distance
    ; travelled from colshift and rowshift; PGM coordinates have varying
    ; km/deg of latitude.
    npix2=MAX(WHERE((vasun NE 0) AND (vdusk NE 0)))	; limit of useful data
    rdpgm,300,plat,plong,glat,glong,/dblsize		; get PGM coords

    ; Get Eureka MLTs for every hour of the current winter
    ; First build name of appropriate MLT file
    IF month GE 10 THEN refyr=year ELSE refyr=year-1
    yryr=STRING(refyr MOD 100,(refyr+1) MOD 100,FORMAT='(2I2)')
    mlt_file='/jasper/cnsr3_data1/pace/eur_mlt_'+yryr+'.dat'
    ; Now read contents of file
    rdcols,mlt_file,nmlt,yr,mo,dy,hr,mi,se,eurmlt
    ; Set up array coordinates of PGM pole and Eureka
    eur_col=160.1
    eur_row=189.4
    pole_col=191
    pole_row=191
    pole_glat=glat(pole_col,pole_row)
    pole_glong=glong(pole_col,pole_row)
    unrotmlt=ATAN(eur_col-pole_col,pole_row-eur_row)*!RADEG/15.
    ; Now loop through all lines read from .mt[gp] file
    FOR i=0,ngood-1 DO BEGIN
        gi=good(i)
        ; Set up the time of the image
        ti=times(gi)
        hour=FLOOR(ti/3600.)
        minute=FLOOR(ti-3600.*hour)/60
        second=ti-(hour*60+minute)*60
        ; Generate the pointer into the Eureka MLT array
        thishr=WHERE((yr EQ year) AND (mo EQ month) AND $
	             (dy EQ day) AND (hr EQ hour),nthishr)
        IF nthishr LE 0 THEN MESSAGE,'No MLT for this hour' ELSE BEGIN
	    ; Interpolate in the Eureka MLT array
	    mlt0=eurmlt(thishr(0))
	    mlt1=eurmlt(thishr(0)+1)
	    IF mlt1 LT mlt0 THEN mlt1=mlt1+24.
	    fhour=(minute*60.+second)/3600.
	    eur_mlt=((1.-fhour)*mlt0+fhour*mlt1) MOD 24
	    ENDELSE
        IF KEYWORD_SET(debug) THEN STOP
        ; "rotangle" is the CCW angle to rotate the image; the plat,
        ; plong arrays implicitly rotate CW under the stationary image.
        rotangle=(eur_mlt-unrotmlt)*15.
        ; Now find the appropriate element in the coord arrays
        rotate_xy,pole_col+colshift(gi),pole_row+rowshift(gi) $
	         ,-rotangle,pole_col,pole_row,shiftcol,shiftrow,/degrees
        IF KEYWORD_SET(debug) THEN STOP
        shiftcol=ROUND(shiftcol)
        shiftrow=ROUND(shiftrow)
        ; Find the geographic coords of the place to which the PGM pole
        ; shifted between images.
        shift_glat=glat(shiftcol,shiftrow)
        shift_glong=glong(shiftcol,shiftrow)
        ; Now work out distance covered and azimuth of motion
        earthrad=relc((pole_glat+shift_glat)/2.)+emis_ht
        ll2rb,pole_glong,pole_glat,shift_glong,shift_glat,adist,draz
	drift_az(gi)=draz
        IF KEYWORD_SET(debug) THEN STOP
        dist=earthrad*adist
        ; Get solar azimuth and determine angle between drift azimuth
        ; and antisunward direction
        sunazel,year,month,day,hour,minute,second,saz $
	       ,loc=[pole_glat,pole_glong]
        IF KEYWORD_SET(debug) THEN STOP
        asaz=(saz+180) MOD 360

        ; Make sure we have some estimate of the time interval between
        ; images
        dt=0.
        IF gi+1 LT npix2-1 THEN dt=times(gi+1)-times(gi)
        edt=0.
        oldv=SQRT(vasun(gi)*vasun(gi)+vdusk(gi)*vdusk(gi))
        IF oldv GT 0 THEN BEGIN
	    csgi=colshift(gi)
	    rsgi=rowshift(gi)
	    edt=pixres*SQRT(csgi*csgi+rsgi*rsgi)/oldv
	    ENDIF
        IF (edt GT 0) AND (dt EQ 0) THEN dt=edt

        ; Now compute 'correct' drift speed and direction
        v(gi)=dist/dt
        daz(gi)=drift_az(gi)-asaz
        mtg(info0+gi)=mtg(info0+gi)+STRING(drift_az(gi),daz(gi),v(gi) $
			      ,FORMAT='(F8.1,F7.1,F7.3)')
        ENDFOR

    ENDIF ELSE BEGIN
    daz=!RADEG*ATAN(vdusk,vasun)
    drift_az=FLTARR(N_ELEMENTS(daz))
    earthrad=relc(80.0533)
    v=SQRT(vasun*vasun+vdusk*vdusk)*111.69*((earthrad+emis_ht)/earthrad)
    FOR i=0,ngood-1 DO BEGIN
	gi=good(i)
        sechms,times(gi),hour,minute,second
        sunazel,year,month,day,hour,minute,second,sunaz,sunel
        drift_az(gi)=(daz(gi)+sunaz+180.) MOD 360.
        mtg(info0+gi)=mtg(info0+gi)+STRING(drift_az(gi),daz(gi),v(gi) $
			                  ,FORMAT='(F8.1,F7.1,F7.3)')
        ENDFOR
    ENDELSE

mvasun=v*COS(daz*!DTOR)
mvdusk=v*SIN(daz*!DTOR)

mtg(info0-1)=mtg(info0-1)+"  Drft_az  Dev'n  Speed"
path=STRMID(path,0,STRLEN(path)-3)+'cmt'+STRLOWCASE(type)
mputfile,path,mtg

flag=BYTARR(npix2)
gap=WHERE((mvasun EQ 0.) and (mvdusk EQ 0.),ngap)
IF ngap GT 0 THEN BEGIN
    flag(gap)=1
    good=REPLICATE(1B,N_ELEMENTS(flag))-flag
    wgood=WHERE(good)
    IF KEYWORD_SET(fillgaps) THEN BEGIN
        FOR i=0,ngap-1 DO BEGIN
	    tmp=extrac(mvasun,gap(i)-3,7)
	    mvasun(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
	    tmp=extrac(mvdusk,gap(i)-3,7)
	    mvdusk(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
            ENDFOR
        ENDIF
    ENDIF

dst=ymd2date(year,month,day,format='d$ n$ Y$')
resp=''
PRINT,'The observations span ',MIN(times/3600.,MAX=maxt),' to ',maxt, $
     ' hours UT',FORMAT='(2(A,F5.2),A)'
READ,'Enter the boundary times for the time panels you wish to plot: ',resp
nplots=nwrds(resp)-1
bt=FLTARR(nplots+1)
READS,resp,bt
; FOR i=0,nplots DO bt(i)=FLOAT(getwrd(resp,i))
IF KEYWORD_SET(ps) $
    THEN psopen,path+'.ps',/long $
    ELSE WINDOW,2,XSIZE=800,YSIZE=1000
!P.MULTI=[0,1,nplots,0,0]
maxv=MAX(mvasun)>1.4
FOR i=0,nplots-1 DO BEGIN
    mset_isoxy,bt(i),bt(i+1),maxv,0. $
	      ,xlock='xmd',ylock='ymx',nxrange=[-0.021,1.021]
    PLOT,[0,1],/NODATA,CHARSIZE=1.6 $
	,TITLE='F-Layer Ionospheric Convection, Eureka, NWT, '+dst $
	,XTITLE='Universal Time (hours)' $
	,YTITLE='Anti-Sunward Convection Velocity (km/s)' $
	,XRANGE=[!X.CRANGE(0),!X.CRANGE(1)],XMINOR=12 $
	,YRANGE=[!Y.CRANGE(0),!Y.CRANGE(1)]
    thisplot=WHERE((bt(i) LE times/3600.) AND (times/3600. LT bt(i+1)))
    ARROW,times(thisplot)/3600. $
	 ,REPLICATE(0.,N_ELEMENTS(thisplot)) $
         ,times(thisplot)/3600.-good(thisplot)*mvdusk(thisplot) $
	 ,good(thisplot)*mvasun(thisplot) $
         ,/DATA,HSIZE=-0.025
    mset_isoxy,0,0,0,0
    ENDFOR

IF NOT KEYWORD_SET(ps) THEN k=GET_KBRD(1)
ERASE

; Now calculate the convection electric field from the convection
; velocity.

; btot=49798.E-9	; T (IGRF 90 magnitude 300 km above Eureka)
; bvert=49775.E-9	; T (vertical component)
; bcoeff=(btot^2)/bvert

; Get the geographic coordinates of the PACE pole (N. hemisphere)
; at the assumed emission height.

PACEpole,emis_ht,gltp,glnp

; Calculate IGRF magnetic induction at PACE geomagnetic pole at
; altitude emis_ht
igrf,gltp,glnp,emis_ht,year,bx,by,bz
bx=1E-4*bx
by=1E-4*by
bz=1E-4*bz
bcoeff=(bx*bx+by*by+bz*bz)/bz

PRINT,'bcoeff=',bcoeff,FORMAT='(A,E10.2)'

ef2dusk=1000.*bcoeff*mvasun*1000.	; mV per meter
ef2noon=1000.*bcoeff*mvdusk*1000.
efmag=SQRT(ef2dusk^2+ef2noon^2)		; field magnitude
efaz=REPLICATE(270.,(N_ELEMENTS(efmag))); field azimuth
efaz(wgood)=90.-ATAN(ef2noon(wgood),-ef2dusk(wgood))*!RADEG
wn=WHERE(efaz LT 0.)
efaz(wn)=efaz(wn)+360.

!P.MULTI=[0,1,nplots,0,0]
ylims=[-30,70]
fac=maxv/(ylims(1)-ylims(0))
isoylims=fac*ylims
FOR i=0,nplots-1 DO BEGIN
    thisplot=WHERE((bt(i) LE times/3600.) AND (times/3600. LT bt(i+1)))
    mset_isoxy,bt(i),bt(i+1),isoylims(0),isoylims(1),nxrange=[-0.021,1.021]
    PLOT,[0,1],/NODATA,CHARSIZE=1.6 $
	,TITLE='Inferred Convection Electric Field, Eureka, NWT, '+dst $
	,XTITLE='Universal Time (hours)' $
	,YTITLE='Electric Field Magnitude (mV/m)' $
	,XRANGE=[!X.CRANGE(0),!X.CRANGE(1)],XMINOR=12 $
	,YRANGE=[!Y.CRANGE(0),!Y.CRANGE(1)]/fac
    OPLOT,times(thisplot)/3600.,efmag(thisplot),THICK=3
    OPLOT,times(thisplot)/3600.,efaz(thisplot)-270.,PSYM=-1
    PLOTS,!X.CRANGE,[0.,0.],LINESTYLE=1
    PLOTS,[0.92,0.96]*(!X.CRANGE(1)-!X.CRANGE(0))+!X.CRANGE(0) $
	 ,[0.93,0.93]*(!Y.CRANGE(1)-!Y.CRANGE(0))+!Y.CRANGE(0) $
	 ,PSYM=-1,/DATA
    XYOUTS,0.90*(!X.CRANGE(1)-!X.CRANGE(0))+!X.CRANGE(0) $
	  ,0.91*(!Y.CRANGE(1)-!Y.CRANGE(0))+!Y.CRANGE(0) $
	  ,'Angle (deg) of !17E!X toward noon from dusk' $
	  ,ALIGNMENT=1,/DATA
    ENDFOR

IF KEYWORD_SET(ps) THEN psclose ELSE BEGIN
    k=GET_KBRD(1)
    WDELETE,0
    WDELETE,2
    ENDELSE

RETURN
END