PRO buildmtg,name,update=update,test=test,color=color,fullscale=fullscale
;+
; 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 ] ]
; 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 set, produce a color PostScript plot; if not
;		set, produce the plot in black and white.
;       FULLSCALE: Optionally sets the full-scale brightness in
;		the montage.  If not set, 300 R is assumed.
; 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(color) EQ 0 THEN color=0B

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.))
; STOP,'PRINT,strsec(times(WHERE(close2hh)))'

; 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))
PRINT,trsh,tcsh,FORMAT= $
      '("Total of ",I4," rows and ",I4," columns of shift")'
montage=BYTARR(tcsh+edge,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))
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
    WINDOW,0,XSIZE=500,YSIZE=1000
    mwin=!D.WINDOW
    ENDIF
loadct,3,/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
    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)
    lastrow=lastrow+rshift(i-1)
    IF KEYWORD_SET(update) THEN BEGIN
	WSET,mwin
	TVSCL,extrac(montage,c0+382-500,lastrow-1000,500,1000)
	ENDIF
    ENDFOR

; No further need for .pgm file

FREE_LUN,data

; 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)
rmin=MIN(wann/wedge,MAX=rmax)
cmin=MIN(wann MOD wedge,MAX=cmax)
ann=(TEMPORARY(ann))(cmin:cmax,rmin:rmax)

MESSAGE,'Grid TVReaD',/INFORMATIONAL

; Shrink montage array to minimum dimensions

wm=WHERE(montage)
rmin=MIN(wm/(tcsh+edge),MAX=rmax)
cmin=MIN(wm MOD (tcsh+edge),MAX=cmax)
montage=(TEMPORARY(montage))(cmin:cmax,rmin:rmax)
pixmax=MAX(montage)

MESSAGE,'Montage shrunk',/INFORMATIONAL

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

asize=SIZE(ann)
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)

; Set up for PostScript output

; psres=125.
psres=250.	; try this value; see how the plot looks
pglen=8*psres
step=0.9*pglen
npages=CEIL(msize(2)/step)
PRINT,'The montage will appear on ',npages,' pages',FORMAT='(A,I0,A)'
IF KEYWORD_SET(color) THEN psfile=ifile+'.cps' ELSE psfile=ifile+'.ps'
IF NOT KEYWORD_SET(test) THEN psopen,psfile,color=color,/long
; loadct,3*KEYWORD_SET(color),/silent	; this should be parameter-dependent
loadct,0,/silent			; this is not parameter-dependent
FOR i=0,npages-1 DO BEGIN
    ; IF i NE npages-1 THEN GOTO,endline
    IF i NE 0 THEN ERASE
    tmp=montage(*,i*step:MIN([i*step+pglen-1,msize(2)-1]))
    wtmp=WHERE(tmp)
    ledge=MIN(wtmp MOD msize(1),MAX=redge)
    maxwid=FLOOR(8.1*psres)
    mid=(redge+ledge)/2
;   tmp=tmp(MAX([mid-maxwid/2,0]):MIN([mid+maxwid/2,msize(1)-1]),*)
    IF NOT KEYWORD_SET(test) THEN $
	mtv,tmp(MAX([mid-maxwid/2,0]):MIN([mid+maxwid/2,msize(1)-1]),*) $
	   ,0,0,res=psres,high=pixmax,/temp
; endline:
    ENDFOR

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)
    SPAWN,'rmps3 '+psfile+' n'
    IF KEYWORD_SET(color) THEN queue='qms' ELSE queue='ps'
    SPAWN,'lpr -l -P'+queue+' -s '+psfile
    SPAWN,'lpq -P'+queue
    ENDIF

RETURN
END