PRO buildmtg,name $
	    ,update=update,test=test,color=color,fullscale=fullscale $
	    ,breakup=breakup,reverse_ct=reverse_ct
;+
; 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 ] $
;              [, FULLSCALE = fullscale ] [, /BREAKUP ] $
;              [, /REVERSE_CT ]
; 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.
; 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.  The file is also queued for
;	printing on the QMS printer ( -Pqms ).
; 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.
;-

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

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='/jasper/cnsr3_data1/pgm/'+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

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

; 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+(trsh/(2*npix2))	; use average rshift
uprow=midrow+rshift(0)/2	; use actual shift for this image
c0=MAX([0,-tcsh])		; position for L edge of first image
montage(c0,0)=(image(firstn))(*,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
IF N_ELEMENTS(color) GT 0 THEN loadct,color,/silent ELSE loadct,0,/silent

; 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))(*,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
	WINDOW,/FREE,XSIZE=xs,YSIZE=ys,/PIXMAP
	pwin=!D.WINDOW
	XYOUTS,0,0,ut,/DEVICE,CHARSIZE=csz,CHARTHICK=cth
	utim=TVRD()
	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)
	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)			; save L and R limits
	    one2two,wmtg,montage,x,y		; 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)=BYTARR(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)			; save L and R limits
    one2two,wmtg,montage,x,y		; 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
WINDOW,/FREE,XSIZE=wedge,YSIZE=5*wedge/2,/PIXMAP
!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,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)

    ; 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

    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

MESSAGE,'Grid created',/INFORMATIONAL

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

ann=TVRD()
WDELETE,!D.WINDOW
wann=WHERE(ann)
armin=MIN(wann/wedge,MAX=armax)
acmin=MIN(wann MOD wedge,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)
    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) GT 0)
poplr=TOTAL(montage(msize(1)-asize(1):msize(1)-1,0:asize(2)-1) GT 0)
x0=(poplr LT popll)*(msize(1)-asize(1))
REPEAT BEGIN
    dx=(TOTAL(montage(x0-1,0:asize(2)-1)) EQ 0)
    x0=x0-dx
    ENDREP UNTIL dx EQ 0
PRINT,'Coordinate grid inserted at column '+STRTRIM(x0,2)
montage(x0,0)=montage(x0,0)+BYTSCL(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 N_ELEMENTS(color) GT 0 THEN psfile=ifile+'.cps' ELSE psfile=ifile+'.ps'
IF NOT KEYWORD_SET(test) THEN BEGIN
    psopen,psfile,/color,/long
    ; IF N_ELEMENTS(color) EQ 0 THEN DEVICE,bits_per_pixel=4
    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

    ; 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)
    one2two,wtmp,tmp,x,y
    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:
    ENDFOR

IF KEYWORD_SET(breakup) THEN FREE_LUN,mtun

montage=0B
IF NOT KEYWORD_SET(test) THEN BEGIN
    psclose,/normps3
    CD,'/jasper/cnsr3_data1/pgm'
    slash=rstrpos(psfile,'/')
    psfile=STRMID(psfile,slash+1,STRLEN(psfile)-slash-1)
    ; STOP
    ; SPAWN,'rmps3 '+psfile+' n'
    ; IF N_ELEMENTS(color) GT 0 THEN queue='qms' ELSE queue='ps'
    ; SPAWN,'lpr -l -P'+queue+' -s '+psfile
    ; SPAWN,'lpq -P'+queue
    ENDIF ELSE PRINT,spc(80)+'Done'

RETURN
END