PRO gridtest,name,update=update,fullscale=fullscale 
 
; Set up OS-dependent parameters. 
 
@isitdos 
 
IF N_ELEMENTS(fullscale) EQ 0 THEN fullscale=300. 
 
; Set the edge length of a single image from the .pgm file. 
 
edge=383 
 
; Load the .mtg file into a STRING array. 
 
path=pgmroot+dd+name 
mtg=qgetfile(path) 
 
; Determine how many images go into the montage. 
 
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)' 
 
; Get the .pgm image resolution 
info0=n-npix2 
READS,mtg(info0-3),pixres,FORMAT='(40X,F6.4)' 
 
; Get the assumed emission height, and work out the spatial scale 
; of the individual images. 
 
htline=0 
WHILE STRMID(mtg(htline),0,7) NE 'Assumed' DO htline=htline+1 
emis_ht=FLOAT(getwrd(mtg(htline),3)) 
 
; Set up the arrays to hold the essential .mtg file contents. 
 
times=FLTARR(npix2) 
cshift=INTARR(npix2) 
rshift=cshift 
vdusk=times 
vasun=times 
p1='' 
 
; Extract the file contents from the STRING array. 
 
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 
    cshift(i)=-p5       ; opposite sign to what is in .mtg file 
    rshift(i)=-p6       ; ditto 
    vasun(i)=p7 
    vdusk(i)=p8 
    ENDFOR 
 
; Identify images close to hours and half-hours 
 
; hh=ROUND(times)/1800                  ; whole # of halfhours since 00 
; close2hh=ABS((hh-SHIFT(hh,1)))<1      ; this is 1 where hh increments 
hh=times/1800. 
close2hh=(ABS(hh-ROUND(hh)) LT (1.5/30.)) 
 
; Identify and fill gaps in the velocity and shift arrays 
 
flag=BYTARR(npix2) 
gap=WHERE((vasun EQ 0.) and (vdusk EQ 0.),ngap) 
IF ngap GT 0 THEN BEGIN 
    flag(gap)=1                         ; flag gaps for later use 
    good=REPLICATE(1B,N_ELEMENTS(flag))-flag 
    FOR i=0,ngap-1 DO BEGIN 
        k=gap(i) 
        tmp=extrac(vasun,k-3,7)         ; first do antisunward vels 
        vasun(k)=TOTAL(tmp)/TOTAL(tmp GT 0)     ; mean of non-0 vels 
        IF k NE npix2-1 $               ; dt=interval between images 
            THEN dt=times(k+1)-times(k) $ 
            ELSE dt=times(k-1)-times(k-2) 
        rshift(k)=ROUND(dt*vasun(k)/pixres)     ; estimate rshift 
        tmp=extrac(vdusk,k-3,7)         ; repeat for duskward vels 
        vdusk(k)=TOTAL(tmp)/TOTAL(tmp GT 0) 
        cshift(k)=ROUND(dt*vdusk(k)/pixres) 
        PRINT,k,cshift(k),rshift(k),vasun(k),vdusk(k),FORMAT= $ 
             '(I3,2X,2(I4,2X),2(F8.5,2X))' 
        ENDFOR 
    ENDIF 
 
csz= ( cth = 2 ) 

ifile=path
CASE STRLOWCASE(STRMID(ifile,STRLEN(ifile)-1,1)) OF
'g': ext='geo'
'p': ext='pgm'
ELSE: MESSAGE,'unknown montage type'
    ENDCASE

; Decorate the montage with appropriate coordinates.  First, create 
; the window. 

; wedge=460 
wedge=400       ; try this to see how the plot looks 
xedge=3*wedge/2
yedge=5*wedge/2
IF KEYWORD_SET(update) $ 
    THEN WINDOW,/FREE,XSIZE=xedge,YSIZE=yedge,/PIXMAP $ 
    ELSE DEVICE,SET_RESOLUTION=[xedge,yedge],Z_BUFFERING=0 
!P.THICK=csz 

; Establish a DATA coordinate system. 

PLOT,INDGEN(wedge)-wedge/2,INDGEN(5*wedge/2)-wedge/2 $ 
    ,/DEVICE,/NODATA,POSITION=[0,0,wedge-1,5*wedge/2-1] $ 
    ,XSTYLE=5,YSTYLE=5,XMARGIN=[0,0],YMARGIN=[0,0] $ 
    ,XRANGE=[-wedge/2,wedge/2],YRANGE=[-wedge/2,2*wedge] 

IF ext EQ 'pgm' THEN BEGIN 

    ; Draw a circle at 85 deg PACE geomagnetic latitude. 

    arcs,5./pixres,0,360,0,0,/device 

    ; Draw meridians every 3 hrs of MLT. 

    FOR i=-5,1 DO radii,0,wedge/2,i*45,0,0,/device 

    ; Draw an extra-long meridian toward noon. 

    radii,0,26./pixres,90,0,0,/device 

    ; Show 5 deg intervals of PGMLat on noon meridian. 

    delta=4. 
    FOR i=2,5 DO arcs,5.*i/pixres,90.-delta/i,90+delta/i,0,0,/device 

    ; Label MLT meridians. 

    mlt_label_offset=[4.,6.,7.,6.,6.,3.,3.,3.] 
    rr=8./pixres 
    FOR i=0,7 DO BEGIN 
        ang=-90.+i*45. 
        dr=0.4*(-SIN(ang*!DTOR)-1.)/pixres 
        xoff=(rr+dr)*COS((ang-mlt_label_offset(i))*!DTOR) 
        yoff=(rr+dr)*SIN((ang-mlt_label_offset(i))*!DTOR) 
        XYOUTS,xoff,yoff,STRTRIM(i*3,2) $ 
              ,ALIGNMENT=0.5,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
        IF i EQ 0 $ 
            THEN XYOUTS,-xoff,yoff,'MLT' $ 
                       ,ALIGNMENT=0.1,CHARSIZE=csz,CHARTHICK=cth 
        ENDFOR 

    ; Label PGMLat parallels. 

    FOR i=1,5 DO BEGIN 
        yoff=(5.*i+0.2)/pixres 
        xoff=-1.5/pixres 
        XYOUTS,xoff,yoff,STRTRIM(90-i*5,2) $ 
              ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
        IF i EQ 1 THEN XYOUTS,0.3/pixres,yoff,'PGMLat' $ 
                             ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
        ENDFOR 

    ENDIF ELSE BEGIN 

    ; How many pixels per 1000 km? 

    PACEpole,emis_ht,glatp,glongp 
    earthrad=relc(glatp) 
    kmperdeg=111.66*(earthrad+emis_ht)/earthrad 
    pixperkm=1./(kmperdeg*pixres) 
    HELP,kmperdeg,pixperkm 

    ; Show 500 km radius and 2500 km distance 

    PLOT,[-600,600],[-600,2400],/DEVICE,/NODATA $ 
        ,POSITION=[wedge/2-600*pixperkm $ 
                  ,wedge/2-600*pixperkm $ 
                  ,wedge/2+600*pixperkm $ 
                  ,wedge/2+2400*pixperkm] $ 
        ,XSTYLE=5,YSTYLE=5,XMARGIN=[0,0],YMARGIN=[0,0] $ 
        ,XRANGE=[-600,600],YRANGE=[-600,2400] 

    ; Draw a circle at 85 deg geographic latitude. 

    arcs,5./pixres,0,360,wedge/2,wedge/2,/device 

    ; Draw an extra-long meridian toward noon. 

    radii,0,26./pixres,90,wedge/2,wedge/2,/device 

    ; Show 5 deg intervals of GGLat on noon meridian. 

    delta=4. 
    FOR i=2,5 DO arcs,5.*i/pixres,90.-delta/i,90+delta/i,wedge/2,wedge/2 $
                     ,/device 

    ; Label PGMLat parallels. 

    FOR i=1,5 DO BEGIN 
        yoff=(5.*i+0.2)/pixres+wedge/2
        xoff=-1.5/pixres+wedge/2
        XYOUTS,xoff,yoff,STRTRIM(90-i*5,2) $ 
              ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
        IF i EQ 1 THEN XYOUTS,0.3/pixres,yoff,'PGMLat' $ 
                             ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
        ENDFOR 

    AXIS,0,0,/DATA,CHARSIZE=csz,CHARTHICK=cth $ 
        ,/XAXIS,XRANGE=[-600,600],/XSTYLE,XTICKS=2,XMINOR=2 
    AXIS,0,0,/DATA,CHARSIZE=csz,CHARTHICK=cth $ 
        ,/YAXIS,YRANGE=[-600,2400],/YSTYLE,YTICKS=5,YMINOR=2 
    !P.THICK=1
    arcs,5.*kmperdeg*pixperkm,0,360,wedge/2,wedge/2,/device
    !P.THICK=cth
 
    ENDELSE 
 
; Write the date, place, emission height, and full scale brightness 
; on the plot somewhere. 
 
dateline=0 
WHILE (isnumber(mtg(dateline)) EQ 0) DO BEGIN 
    dateline=dateline+1 
    IF dateline GT n THEN MESSAGE,'Trouble finding date!' 
    ENDWHILE 
date=STRTRIM(mtg(dateline),2) 
XYOUTS,0,0,date,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
 
MESSAGE,'Got date',/INFORMATIONAL 
 
place='Eureka 88.7!7K!X' 
XYOUTS,0,csz*(1.1*!D.Y_CH_SIZE),place $ 
      ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
 
XYOUTS,wedge-2*csz*!D.X_CH_SIZE,csz*(1.1*!D.Y_CH_SIZE) $ 
      ,STRING(emis_ht,FORMAT='(I3," km")') $ 
      ,ALIGNMENT=1,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
 
fsline=htline 
WHILE STRPOS(mtg(fsline),'Full scale') LT 0 DO fsline=fsline+1 
IF fsline GT 0 THEN BEGIN 
    fullscale=ROUND(FLOAT(getwrd(mtg(fsline),4))) 
    XYOUTS,wedge-2*csz*!D.X_CH_SIZE,1 $ 
          ,STRING(fullscale,FORMAT='("<",I3," R")') $ 
          ,ALIGNMENT=1,/DEVICE,CHARSIZE=csz,CHARTHICK=cth 
    ENDIF 

; Try putting on a colour bar

ghi=STRING(!D.N_COLORS-1,FORMAT='(I3)')
phi=STRING(ROUND(fullscale),FORMAT='(I3)')
bar,layout=['400','300','450','700','0',ghi,'0',phi $
           ,'0',phi,'100.','(I3)','Brightness','2','2' $
           ,ghi,ghi,'2',ghi,'2']
 
MESSAGE,'Grid created',/INFORMATIONAL 
 
; Read back the display as an image and crop it. 
 
ann=TVRD() 
MESSAGE,'Grid TVReaD',/INFORMATIONAL 
 
IF KEYWORD_SET(update) THEN WDELETE,!D.WINDOW 
wann=WHERE(ann) 
armin=MIN(wann/xedge,MAX=armax) 
acmin=MIN(wann MOD xedge,MAX=acmax) 
ann=(TEMPORARY(ann))(acmin:acmax,armin:armax) 

HELP,ann
sann=SIZE(ann)
PRINT,sann
WINDOW,/FREE,XSIZE=sann(2),YSIZE=sann(1)
TV,ROTATE(ann,1)
 
RETURN 
END