;+
; NAME:
;       BUILDMTG
; PURPOSE:
;       Creating a montage of Polar Camera images that show F-layer
;       patches convecting over Eureka.
; CATEGORY:
;
; CALLING SEQUENCE:
;       BUILDMTG, MTGNAME [, /UPDATE ] [, /TEST ] [, COLOR = color ] $
;              [, FULLSCALE = fullscale ] [, /BREAKUP ] $
;              [, /REVERSE_CT ] [, /MULTIPAGE ]
; INPUTS:
;       MTGNAME: The name of a .mtg file (located in ~cnsr3/pgm)
;               that contains the information needed to build the
;               montage.  This includes the numbers of columns and
;               rows to shift between adjacent images.  This file
;               is produced by the IDL routine MONTAGE.PRO.
; KEYWORD PARAMETERS:
;       /UPDATE: If set, specifies that the MONTAGE array is to be
;               displayed after each addition to it, to keep the
;               user apprised of the procedure's progress.  This
;               option slows the procedure down by a factor of at
;               least 3.
;       /TEST:  If set, the montage array is built but not printed
;               to a PostScript file.  This allows debugging to
;               proceed efficiently.
;       COLOR:  If defined, gives the color table to be used in
;               producing a color PostScript plot; if not defined,
;               the plot is produced in black and white.
;       FULLSCALE: Optionally sets the full-scale brightness in
;               the montage.  If not set, 300 R is assumed.
;       /BREAKUP: If set, stash the montage array into an ASSOCiated
;               file every 250 rows to avoid overflowing virtual
;               memory.
;       /REVERSE_CT: If set, reverses the colour table to plot the
;               montage as a 'negative', i.e., with the brightest
;               regions appearing darkest in the plot.
;        /MULTIPAGE: If set, writes the output to as many PostScript
;                files as there are pages to be printed, with one page
;                per file.
; OUTPUTS:
;       None.
; COMMON BLOCKS:
;       None.
; SIDE EFFECTS:
;       A colour PostScript file is produced in ~cnsr3/pgm, containing
;       the montage array separated into 8" sections with 0.8" of
;       overlap between sections.
; RESTRICTIONS:
;       An image file having the same name as the .mtg file, but
;       with extension .pgm instead, must exist in ~cnsr3/pgm.
;       This file must contain the Polar Camera images, transformed
;       to PACE geomagnetic coordinates, that will be assembled
;       into a montage.
; PROCEDURE:
;       Later, 'gator.
; MODIFICATION HISTORY:
;       Written 94/06/20 by D P Steele.
;       Documented 94/06/21 by D P Steele.
;-

PRO buildmtg,name $
            ,update=update,test=test,color=color,fullscale=fullscale $
            ,breakup=breakup,reverse_ct=reverse_ct,multipage=multipage

IF N_PARAMS() EQ 0 THEN BEGIN
    doc_library,'buildmtg'
    RETURN
    ENDIF

; 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

; Set the background pixel value appropriate for the colour table to be used.

background=(1B-KEYWORD_SET(reverse_ct))*(!D.N_COLORS-1)

; 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)'

; Read source camera and filter numbers from file.

cfline=0
WHILE STRMID(mtg(cfline),0,12) NE 'Images from ' DO cfline=cfline+1
cam=STRMID(mtg(cfline),19,1)
filt=STRMID(mtg(cfline),29,1)

; 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=times/1800.
close2hh=(ABS(hh-ROUND(hh)) LT (1.5/30.))	; i.e., within 1.5 min

; 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

; Decide how big montage has to be, given net shifts found

trsh=FIX(TOTAL(rshift))
tcsh=FIX(TOTAL(cshift))
width=tcsh+edge
PRINT,trsh,tcsh,FORMAT= $
      '("Total of ",I4," rows and ",I4," columns of shift")'

IF KEYWORD_SET(breakup) THEN BEGIN      ; setup to save montage(*,0:249)
    montage=BYTARR(width,500)           ; to disk every time it's filled
    OPENW,mtun,'montage.tmp',/GET_LUN
    mtgsave=ASSOC(mtun,BYTARR(width,250))
    nsaved=0
    pixmax=0
    cmin=width-1
    cmax=0
    ENDIF ELSE montage=BYTARR(width,trsh+edge)
montage(*,*)=background                        ; white background

; Open .pgm file that contains the transformed images

ifile=path
CASE STRLOWCASE(STRMID(ifile,STRLEN(ifile)-1,1)) OF
'g': ext='geo'
'p': ext='pgm'
ELSE: MESSAGE,'unknown montage type'
    ENDCASE
STRPUT,ifile,ext,STRLEN(ifile)-3
OPENR,data,ifile,/GET_LUN
image=ASSOC(data,BYTARR(edge,edge))

; Set up to begin filling montage array.

midrow=edge/2
uprow=midrow+rshift(0)/2        ; use actual shift for this image
c0=MAX([0,-tcsh])               ; position for L edge of first image
chunk=image(firstn)
chunk(outside(chunk,dil=2,ero=1))=background        ; make image background white
montage(c0,0)=chunk(*,0:uprow)
c0=c0+cshift(0)                 ; L edge moves by cshift(0)
lastrow=uprow                   ; bottom unfilled row of montage

IF KEYWORD_SET(update) THEN BEGIN       ; display progress
    IF KEYWORD_SET(breakup) THEN ylen=250 ELSE ylen=1000
    WINDOW,0,XSIZE=500,YSIZE=ylen
    mwin=!D.WINDOW
    ENDIF ELSE BEGIN
    old_dev=!D.NAME
    SET_PLOT,'Z'
    IF dos THEN BEGIN
        psfile=ifile
        pslen=STRLEN(psfile)

        ; Inject the plot type description into the existing suffix

        IF N_ELEMENTS(color) GT 0 THEN suffix='cs' ELSE suffix='ps'
        STRPUT,psfile,suffix,pslen-2

        ; Compress the cam/filt information into a digit 0..9 and
        ; use the other digit (later) to index the plot file.

        cfnum=STRTRIM(STRING(5*FIX(cam)+FIX(filt)),2)
        STRPUT,psfile,cfnum,pslen-6
        ENDIF ELSE BEGIN
        IF N_ELEMENTS(color) GT 0 THEN psfile=ifile+'.cps' $
                                  ELSE psfile=ifile+'.ps'
        ENDELSE
    HELP,psfile
    ENDELSE

IF N_ELEMENTS(color) GT 0 THEN loadct,color,/silent ELSE loadct,0,/silent
!P.COLOR=(!D.N_COLORS-1)-background                ; black text

; Add each image to the montage

csz=2                           ; big characters for time labels
cth=csz                         ; big, THICK characters!
FOR i=1,npix2-1 DO BEGIN
    chunk=image(firstn+i)
    ; The following line is a temporary debugging measure (DPS, 960221)
;   flag=(48 LE i) AND (i LE 55)
;   chunk(outside(chunk,dil=2,ero=1,view=flag))=background        ; white
    chunk(outside(chunk,dil=2,ero=1))=background        ; white
    chunk=chunk(*,uprow+1-rshift(i-1):uprow)
    ; only paste rshift rows from current image
;   montage(c0:c0+edge-1,lastrow+1:lastrow+rshift(i-1))= $
;       montage(c0:c0+edge-1,lastrow+1:lastrow+rshift(i-1))+chunk
    montage(c0,lastrow)=chunk   ; MUCH neater, if it works
    PRINT,i,MAX(chunk),FORMAT='(I3,I6)'
    IF close2hh(i) THEN BEGIN                           ; add time label

; Turn time string into an image to paste into the montage

        ut=strsec(times(i)+.5)
        xs=csz*STRLEN(ut)*!D.X_CH_SIZE
        ys=csz*!D.Y_CH_SIZE
        IF KEYWORD_SET(update) THEN BEGIN        ; implies (X) Windows device
            WINDOW,/FREE,XSIZE=xs,YSIZE=ys,/PIXMAP
            pwin=!D.WINDOW
            ENDIF ELSE DEVICE,SET_RESOLUTION=[xs,ys],Z_BUFFERING=0
        ERASE,background
        XYOUTS,0,0,ut,/DEVICE,CHARSIZE=csz,CHARTHICK=cth
        utim=TVRD()
        IF KEYWORD_SET(update) THEN WDELETE,pwin
        x0=MAX([0,c0-xs])
        y0=lastrow+1
;       montage(x0:x0+xs-1,y0:y0+ys-1)= $
;           montage(x0:x0+xs-1,y0:y0+ys-1)+BYTSCL(utim,TOP=MAX(image(i)))
;        montage(x0:x0+xs-1,y0:y0+ys-1)= $
;            montage(x0:x0+xs-1,y0:y0+ys-1)+BYTSCL(utim)
        montage(x0,y0)=BYTSCL(utim)
        ENDIF

; Update L edge and bottom unfilled row pointers

    c0=c0+cshift(i)
    uprow=uprow-rshift(i-1)/2+rshift(i)/2
    lastrow=lastrow+rshift(i-1)

    IF KEYWORD_SET(breakup) THEN BEGIN
        IF lastrow GE 250 THEN BEGIN    ; save 250 rows of montage
            mtgsave(nsaved)=montage(*,0:249)    ; copy 250 rows to disk
            pixmax=MAX([pixmax,MAX(montage)])   ; highest pixel value
            wmtg=WHERE(montage NE background)        ; save L and R limits
            smtg=SIZE(montage) & x=wmtg MOD smtg(1) ; of montage as
            cmin=MIN([cmin,MIN(x,MAX=xmax)])    ; cmin and
            cmax=MAX([cmax,xmax])               ; cmax
            nsaved=nsaved+1                     ; increment record pointer
            PRINT,nsaved
            montage(0,0)=montage(*,250:499)     ; move remaining rows down
            montage(0,250)=REPLICATE(background,width,250) ; and zero out old position
            lastrow=lastrow-250
            ENDIF
        ENDIF

    IF KEYWORD_SET(update) THEN BEGIN
        WSET,mwin
        TVSCL,extrac(montage,c0+382-500,lastrow-ylen,500,ylen)
        ENDIF
    ENDFOR

; No further need for .pgm file

FREE_LUN,data

IF KEYWORD_SET(breakup) THEN BEGIN      ; save last part of montage array
    mtgsave(nsaved)=montage(*,0:249)
    nsaved=nsaved+1
    PRINT,nsaved
    mtgsave(nsaved)=montage(*,250:499)
    nsaved=nsaved+1
    PRINT,nsaved
    wmtg=WHERE(montage NE background)   ; save L and R limits
    smtg=SIZE(montage) & x=wmtg MOD smtg(1) ; of montage as
    cmin=MIN([cmin,MIN(x,MAX=xmax)])    ; cmin and
    cmax=MAX([cmax,xmax])               ; cmax
    pixmax=MAX([pixmax,MAX(montage)])
    FREE_LUN,mtun                       ; and close the ASSOC file
    ENDIF

; 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
!P.COLOR=(!D.N_COLORS-1)-background
!P.BACKGROUND=background
ERASE,background

; 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,edge/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,[-500,500],[-500,2000],/DEVICE,/NODATA $
        ,POSITION=[wedge/2-500*pixperkm $
                  ,wedge/2-500*pixperkm $
                  ,wedge/2+500*pixperkm $
                  ,wedge/2+2000*pixperkm] $
        ,XSTYLE=5,YSTYLE=5,XMARGIN=[0,0],YMARGIN=[0,0] $
        ,XRANGE=[-500,500],YRANGE=[-500,2000]

    AXIS,0,0,/DATA,CHARSIZE=csz,CHARTHICK=cth $
        ,/XAXIS,XRANGE=[-500,500],/XSTYLE,XTICKS=2,XMINOR=2
    AXIS,0,0,/DATA,CHARSIZE=csz,CHARTHICK=cth $
        ,/YAXIS,YRANGE=[-500,2000],/YSTYLE,YTICKS=5,YMINOR=2
;   arcs,5.*kmperdeg*pixperkm,0,360,0,0,/data

    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)')
tclr=STRING((!D.N_COLORS-1)-background,FORMAT='(I0)')
bar,layout=['400','300','450','700','0',ghi,'0',phi $
           ,'0',phi,'100.','(I3)','Rayleighs','2','2' $
           ,tclr,tclr,'2',tclr,'2']

MESSAGE,'Grid created',/INFORMATIONAL

; Read back the display as an image and crop it.

ann=TVRD()
IF KEYWORD_SET(update) THEN WDELETE,!D.WINDOW
wann=WHERE(ann NE background)
armin=MIN(wann/xedge,MAX=armax)
acmin=MIN(wann MOD xedge,MAX=acmax)
ann=(TEMPORARY(ann))(acmin:acmax,armin:armax)

MESSAGE,'Grid TVReaD',/INFORMATIONAL

; Shrink montage array to minimum dimensions

IF NOT KEYWORD_SET(breakup) THEN BEGIN
    wm=WHERE(montage NE background)
    rmin=MIN(wm/(width),MAX=rmax)
    cmin=MIN(wm MOD (width),MAX=cmax)
    montage=(TEMPORARY(montage))(cmin:cmax,rmin:rmax)
    pixmax=MAX(montage)
    MESSAGE,'Montage shrunk',/INFORMATIONAL
    ENDIF

; Insert coordinate grid where it will fit best, i.e., obstruct as
; few filled pixels as possible.  Try both bottom corners.

asize=SIZE(ann)

IF KEYWORD_SET(breakup) THEN BEGIN      ; recover needed part of
    nrecs=(asize(2)/250)+1              ; montage array to insert
    montage=BYTARR(width,nrecs*250) ; annotation
    OPENU,mtun,'montage.tmp',/GET_LUN
    mtgsave=ASSOC(mtun,BYTARR(width,250))
    FOR i=0,nrecs-1 DO montage(0,250*i)=mtgsave(i)
    ENDIF

IF KEYWORD_SET(breakup) THEN BEGIN
    msize=[2L,width,(nsaved-1)*250L+lastrow-1,1L]
    msize=[msize,msize(1)*msize(2)]
    ENDIF ELSE msize=SIZE(montage)

popll=TOTAL(montage(0:asize(1)-1,0:asize(2)-1) NE background)
poplr=TOTAL(montage(msize(1)-asize(1):msize(1)-1,0:asize(2)-1) NE background)
x0=(poplr LE popll)*(msize(1)-asize(1))
REPEAT BEGIN
    dx=(TOTAL(montage(x0-1,0:asize(2)-1) EQ background) EQ asize(2))
    x0=x0-dx
    ENDREP UNTIL dx EQ 0
PRINT,'Coordinate grid inserted at column '+STRTRIM(x0,2)
montage(x0,0)=BYTSCL(ann)
ann=[0B]
wann=ann


IF KEYWORD_SET(breakup) THEN BEGIN      ; save modified part
    FOR i=0,nrecs-1 DO mtgsave(i)=montage(*,250*i:250*(i+1)-1)
    FREE_LUN,mtun
    ENDIF

; Set up for PostScript output

; psres=125.
psres=250.      ; try this value; see how the plot looks
maxwid=FLOOR(8.1*psres)
pglen=8*psres
step=0.9*pglen
npages=CEIL((msize(2)-pglen)/step)+1
PRINT,'The montage will appear on ',npages,' pages',FORMAT='(A,I0,A)'

IF NOT KEYWORD_SET(test) THEN BEGIN
    IF NOT KEYWORD_SET(multipage) THEN BEGIN
        thispsfile=psfile
        psopen,thispsfile,color=(N_ELEMENTS(color) GT 0),/long
        ENDIF
    ENDIF

IF N_ELEMENTS(color) GT 0 THEN loadct,color,/silent ELSE loadct,0,/silent
IF KEYWORD_SET(reverse_ct) THEN BEGIN
    TVLCT,r,g,b,/GET
    r=reverse(r)
    g=reverse(g)
    b=reverse(b)
    TVLCT,r,g,b
    PRINT,'Color table reversed'
    ENDIF

IF KEYWORD_SET(breakup) THEN BEGIN
    OPENR,mtun,'montage.tmp',/GET_LUN
    mtgsave=ASSOC(mtun,BYTARR(width,250))
    ENDIF

FOR i=0,npages-1 DO BEGIN
    ; IF i NE npages-1 THEN GOTO,endline
    IF (NOT KEYWORD_SET(test)) AND KEYWORD_SET(multipage) THEN BEGIN
        thispsfile=psfile
        IF dos $
            THEN STRPUT,thispsfile,STRING(i,FORMAT='(I0)'),pslen-5 $
            ELSE thispsfile=thispsfile+STRING(i,FORMAT='(".",I0)')
        HELP,thispsfile
        psopen,thispsfile,color=(N_ELEMENTS(color) GT 0),/long
        ENDIF

    ; Determine upper and lower bounds of this page of plot

    firstline=i*step
    lastline=MIN([i*step+pglen-1,msize(2)-1])
    IF KEYWORD_SET(breakup) THEN BEGIN
        ; read in records from montage.tmp to restore page of montage
        firstrec=LONG(firstline)/250
        lastrec=(LONG(lastline)/250)<(nsaved-1)
        tmp=BYTARR(width,(lastrec-firstrec+1)*250)
        FOR ip=firstrec,lastrec DO tmp(0,250*(ip-firstrec))=mtgsave(ip)
        firstline=firstline-firstrec*250
        lastline=lastline-firstrec*250
        ENDIF ELSE BEGIN
        tmp=montage(*,firstline:lastline)
        lastline=lastline-firstline
        firstline=0
        ENDELSE

    ; Determine left and right bounds of this page of plot

    wtmp=WHERE(tmp NE background)
;   one2two,wtmp,tmp,x,y
    stmp=SIZE(tmp)
    x=wtmp MOD stmp(1)
    ledge=MIN(x,MAX=redge)
    ; ledge=MIN(wtmp MOD msize(1),MAX=redge)
    mid=(redge+ledge)/2
    ledge=MAX([mid-maxwid/2,0])
    redge=MIN([mid+maxwid/2,msize(1)-1])

    IF KEYWORD_SET(test) THEN BEGIN
        xsiz=redge-ledge+1
        ysiz=lastline-firstline+1
        fac=CEIL(xsiz/1000.)>CEIL(ysiz/1000.)
        IF (xsiz MOD fac) GT 0 THEN redge=redge-(xsiz MOD fac)
        IF (ysiz MOD fac) GT 0 THEN lastline=lastline-(ysiz MOD fac)
        IF i EQ 0 THEN WINDOW,2,XSIZE=1000,YSIZE=1000
        TVSCL,REBIN(tmp(ledge:redge,firstline:lastline) $
                   ,xsiz/fac,ysiz/fac)
        PRINT,spc(80)+'Page '+STRTRIM(i,2)
        ENDIF ELSE BEGIN
        IF i NE 0 THEN ERASE
        mtv,tmp(ledge:redge,firstline:lastline) $
           ,0,0,res=psres,high=pixmax,/temp
        ENDELSE
; endline:
    IF KEYWORD_SET(multipage) THEN psclose,/normps3
    ENDFOR

IF KEYWORD_SET(breakup) THEN FREE_LUN,mtun

montage=0B
IF NOT KEYWORD_SET(test) THEN BEGIN
    IF NOT KEYWORD_SET(multipage) THEN psclose,/normps3
    CD,pgmroot
    slash=rstrpos(thispsfile,'/')
    psfile=STRMID(thispsfile,slash+1,STRLEN(psfile)-slash-1)
    ENDIF ELSE PRINT,spc(80)+'Done'

IF N_ELEMENTS(old_dev) EQ 1 THEN SET_PLOT,old_dev

RETURN
END