; SHOWRAY.PRO
;
; Originally written to show images from Polar Camera initial field test
; at Ryning Farm Observatory
; Modified 1993/1/29 by DPS to handle images from image tree structure
; on /jasper/cnsr3_data1
; Modified 1994/2/7 by DPS to read images from associated files,
; where the latter are formatted as 512-byte headers followed by
; 256x256 bytes.  The images are scaled to 8 bits for storage, and
; are reconstituted by multiplying by 2^p, where the value of p is
; stored in byte 290 of the header (counting from 0).
;

COMMON date,year,month,day
COMMON headbytes,oldfname,hbytes,utsec,exptimes,fwtemps,bad_image

; create byte-valued input parameters
cam=0 & filt=cam & day=cam & month=cam
@isitdos
s=''
psplot=0B
take2=0B
!ORDER=1
btype=['','LOG']
IF !D.WINDOW EQ 0 THEN WDELETE,0

; tell user what imaging modes and filters are available
PRINT,"Num   Cam0 Filters            Cam1 Filters"
wavelengths=STRARR(5)
FOR i=0,4 DO BEGIN
    wv=STRTRIM(filters(i),2)
    tmp=STRMID(wv,0,STRPOS(wv,'nm')+2)
    wavelengths(i)=tmp+spc(24-STRLEN(tmp))
ENDFOR
col2=wavelengths & maxlen2=MAX(STRLEN(col2))
FOR i=0,4 DO BEGIN
    wv=STRTRIM(filters(i+5),2)
    tmp=STRMID(wv,0,STRPOS(wv,'nm')+2)
    wavelengths(i)=tmp+spc(24-STRLEN(tmp))
ENDFOR
col3=wavelengths & maxlen3=MAX(STRLEN(col3))
FOR i=0,4 DO BEGIN
    line=STRING(i,col2(i),col3(i),FORMAT='(I1,5X,2A24)')
    PRINT,line
    ENDFOR

; prompt user and get input
REPEAT BEGIN
    READ,'Enter the year, month, day, camera, and filter desired: ' $
        ,year,month,day,cam,filt
    ENDREP UNTIL ((cam EQ 0) OR (cam EQ 1)) AND $
                 ((0 LE filt) AND (filt LE 4))

hour=0
minute=0
second=0
sthms=''
READ,'Enter the hour, minute, and second to start (default 00:00:00 UT): ' $
    ,sthms
IF sthms NE '' THEN BEGIN
    nw=nwrds(sthms)
    IF nw GE 1 THEN hour=FIX(getwrd(sthms,0))
    IF nw GE 2 THEN minute=FIX(getwrd(sthms,1))
    IF nw GE 3 THEN second=FIX(getwrd(sthms,2))
    ENDIF ELSE BEGIN
    hour=0
    minute=0
    second=0
    startut=0L
    ENDELSE

; ensure a 4-digit year for display; make reasonable assumptions
IF year LT 1900 THEN BEGIN
    IF year LT 50 THEN year=year+2000
    IF year GT 50 THEN year=year+1900
    ENDIF

; create title for image plot
wtitle=STRING(filters(cam*5+filt),FORMAT='("Polar Camera",2X,A)')
wtitle=wtitle+'  '+ymd2date(year,month,day,format='d$ n$ y$')

; Display only those pixels viewing the sky
resp=''
Read,Prompt='Was the moonshade in place? ([Y]/N) ',resp
moonshade=(StrUpCase(resp) NE 'N')
IF moonshade THEN BEGIN
    getsig,year,month,day,mask0,mask1
    CASE cam of
     0: mask=mask0
     1: mask=mask1
        ENDCASE
    ENDIF ELSE BEGIN
    mk,mask
    mask(*)=1
    ENDELSE
wmask=WHERE(mask)
nwmask=TOTAL(mask)
PRINT,nwmask,FORMAT='("Mask has ",F6.0," elements.")'

; define path to desired images
template=filocate(year,month,day,cam,filt)

IF template EQ '' THEN BEGIN
    MESSAGE,'No files found; quitting',/INFORMATIONAL
    textout=0B                ; give it a value so quitting won't cause an error
    GOTO,DONE
    ENDIF

; Open the file as a series of small (header-size) units and get
; the header information

nimages=rdhbytes(template)
IF N_ELEMENTS(bad_image) GT 1 THEN BEGIN
    PRINT,'More than 1 bad image; override and display? ([Y]/N) ',resp
    IF STRUPCASE(resp) NE 'N' THEN bad_image=[-1]
    ENDIF
; STOP

; Here we just load (or create) the data to correct for the reduced
; transmission.

reltrans,oldfname,filtemp,filtrans

; Calculate solar azimuths for all images now, to get ahead of the game
uth=utsec/3600
utm=(utsec-3600*uth)/60
uts=utsec-3600*uth-60*utm
sunazel,year,month,day,uth,utm,uts,sunaz

; identify image sequences and prompt user for sequence(s) to display
dt=ABS(utsec-SHIFT(utsec,1))        ; rectify big negative number in front
seqstart=WHERE(dt GT 600)        ; semi-arbitrary large number (10 min)

RESTART:
; Set plotting device
READ,'Make a PostScript plot file? (Y/[N]) ',resp
psplot=(STRUPCASE(resp) EQ 'Y')
colfac=0.8+0.195*psplot                ; 0.8*!D.N_COLORS for X, !D.N_COLORS-1 for PS
IF psplot THEN BEGIN
    READ,'B/W or Colour? ([B]/C) ',resp
    colour=(STRUPCASE(resp) EQ 'C')
    ENDIF

; Offer option to create a text file tabulating "zenith" brightnesses.
text=''
READ,'Save zenith brightnesses to a text file? (Y/[N]) ',text
textout=(STRUPCASE(text) EQ 'Y')
IF textout THEN READ,'Enter the maximum zenith angle to include: ',maxza

IF seqstart(0) GE 0 THEN BEGIN
    PRINT,'Sequences begin with the following image numbers(s):'
    PRINT,seqstart
    ENDIF

mnex=MIN(exptimes,MAX=mxex)
IF mnex NE mxex THEN BEGIN
    PRINT,mnex,mxex $
    ,FORMAT='("Exposure times range from ",F4.1," to ",F4.1," s")'
    texp=''
    READ,'Enter the exposure time of interest (<Return> for all): ',texp
    selexp=isnumber(texp,expt)
    ENDIF ELSE BEGIN
    selexp=0
    expt=FIX(ROUND(mnex))
    ENDELSE

; Would a van Rhijn correction help?
READ,'Apply a van Rhijn correction? (Y/[N]) ',resp
vR=(STRUPCASE(resp) EQ 'Y')
IF vR THEN BEGIN
    READ,'Enter the layer height (km): ',vRh
    READ,'Enter the layer zenith brightness (R): ',vRb
    ; Set up azimuth and zenith angle images for the specified camera
    ; and site
    PCsite=PoCasite(year,month,day)
    col=BINDGEN(256)#REPLICATE(1B,1,256)
    row=TRANSPOSE(col)
    cr2az,year,month,day,cam,filt,col,row,az,za,/deg,site=PCsite
    v=vanRhijn(za,vRh)*vRb

    ; Do a test to see how well we're estimating the brightness

    OPENW,vrun,'showray.vrb',/GET_LUN
    PRINTF,vrun,year,month,day,cam,filt,hour,minute,second,vRh,vRb $
               ,FORMAT='(I4,2I3,2I2,3I3,2F6.1)'
    FREE_LUN,vrun
    WINDOW,3,XSIZE=256,YSIZE=256

    ENDIF ELSE v=FLTARR(256,256)        ; subtract nothing

IF NOT take2 THEN BEGIN
    IF sthms EQ '' THEN BEGIN
        instring=''
        READ,'Enter the range of image numbers to display (<Return> for all images): ',instring
        f=isnumber(instring,inf)
        IF f NE 0 THEN BEGIN                ; user entered at least one number
            inf=FIX(inf)
            IF f GT 0 $                        ; user entered only one number
                THEN inl=nimages-1 $
                ELSE BEGIN
                    f=isnumber(getwrd(instring,1),inl)
                    IF f NE 0 THEN inl=FIX(inl) ELSE inl=nimages-1
                ENDELSE
            ENDIF ELSE BEGIN
            inf=0
            inl=nimages-1
            ENDELSE
        ENDIF ELSE BEGIN
        startut=(hour*60L+minute)*60L+second
        inf=MIN(WHERE(utsec GE startut))
        inl=nimages-1
        ENDELSE
    ENDIF

; Set up image display
X_charheight=10                ; try to keep scalings that apply in X windows
X_charwidth=8
barwidth=16
ninrow=4
nincol=3
ninwin=ninrow*nincol
wwidth=ninrow*256+barwidth
wheight=nincol*256+4*X_charheight

; Create PostScript plot file name and open the file
IF psplot THEN BEGIN
    IF !D.WINDOW EQ 0 THEN WDELETE,0
    IF startut EQ 0 THEN BEGIN
        stheader=gethd(hbytes(*,inf(0)))
        sth=FIX(stheader.misc.tm.hour)
        stm=FIX(stheader.misc.tm.minute)
        sts=FIX(stheader.misc.tm.second)
        ENDIF ELSE BEGIN
        sth=hour
        stm=minute
        sts=second
        ENDELSE
;
;   A new algorithm for naming PS files is needed because DOS doesn't
;   allow enough characters for an unambiguous file name.  The solution
;   is to store YYMMDD, 5*C+F, and a serial number that distinguishes
;   different PS files for the same date and CAM/FILT combination.  The
;   serial number is derived from the number of existing files that start
;   with the same 7 characters.  The start and exposure time information
;   are stored in the PS file as a comment using the DEVICE,OPTION='...'
;   command.
;   DPS, 950817

    psfile=STRING(psroot,dd,year MOD 100,month,day,5*cam+filt $
                 ,FORMAT='(A,A1,3I2.2,I1)')
    already=rfindfile(psfile+'*.*',count=ser)
    psfile=psfile+STRING(ser,FORMAT='(I1)')

    IF KEYWORD_SET(colour) THEN suf='.cps' ELSE suf='.ps'
    psfile=psfile+suf
    psopen,psfile,/landscape,color=colour

;   The following command stores the start and exposure times in the PS
;   file as a comment string.

    DEVICE,OUTPUT=STRING(sth,stm,sts,expt,FORMAT= $
                        '("%%! Start ",3I2.2," UT, ",I2," s exposure")')

    psres=MAX([wwidth/8.9,wheight/(!D.Y_SIZE/2540.)])
    ENDIF ELSE BEGIN
    WINDOW,0,XSIZE=wwidth,YSIZE=wheight
    TVSCL,BINDGEN(256)        ; wake up the X window and set !D.N_COLORS correctly
    IF dos THEN WSHOW,0,/ICONIC
    psres=100.                ; dummy value
    ENDELSE

; Open file to hold the zenith brightnesses.
IF textout THEN BEGIN
    zname=STRING(year MOD 100,month,day,5*cam+filt $
                ,FORMAT='(3I2.2,I1)')
    zfiles=rfindfile(zbroot+dd+zname+'*',count=nzf)
    zname=zname+STRING(nzf,FORMAT='(I1,".TXT")')
    OPENW,zlun,zbroot+dd+zname,/GET_LUN
    PRINTF,zlun,year,month,day,cam,filt,maxza $
               ,FORMAT='(";",I5,2I3,2I2,F5.1," deg half-angle"/)'

    ; Set up the arrays needed to sample the zenith angle cone.

    IF N_ELEMENTS(az) EQ 0 THEN BEGIN
        PCsite=PoCasite(year,month,day)
        col=BINDGEN(256)#REPLICATE(1B,1,256)
        row=TRANSPOSE(col)
        cr2az,year,month,day,cam,filt,col,row,az,za,/deg,site=PCsite
        inview=WHERE(za LT maxza,ninview)
        PRINT,maxza,ninview $
             ,FORMAT='(F4.1," deg half-angle field of view includes "' $
                    +',I0," pixels.")'
        ENDIF
    ENDIF

; Now load the colour table
cf=5*cam+filt
CASE cf OF
    2: ct=3
    7: ct=1
    8: ct=8
    ELSE: ct=0
    ENDCASE
nct=''
READ,'Enter desired colour table (default='+STRTRIM(ct,2)+'): ',nct
IF isnumber(nct) THEN BEGIN
    nct=FLOAT(nct)
    ct=FIX(ABS(nct))
    loadct,ct,/silent
    IF nct LT 0 THEN BEGIN
        READ,'Enter the gamma correction factor: ',gamma
        stretch,0,!D.N_COLORS-1,gamma
        ENDIF
    ENDIF ELSE xloadct
whitetop

; Offer option of plotting "negative" images
IF ct EQ 0 THEN BEGIN
    Read,Prompt='Do you want to plot images as negatives? (Y/[N]) ',resp
    neg=(StrUpCase(resp) EQ 'Y')
    ENDIF ELSE neg=0

; Prompt user for fixed saturation brightness
tops=''
PRINT,'If you want to assign the top of the color bar to a fixed number of counts,'
PRINT,'enter that number now, otherwise enter anything else for auto-scaling.'
READ,tops
flag=isnumber(tops,topnum)        ; topnum=0 if tops doesn't represent a number
IF topnum(0) LT 0 THEN topnum=LONG(topnum(0))+(2L^16)
topnum=topnum(0)

; Prompt user for brightness scale type (log/lin; default=lin)
baxtyp=''
READ,'Plot brightnesses logarithmically? (Y/[N]) ',baxtyp
baxis=(STRUPCASE(baxtyp) EQ 'Y')        ; 0 => linear, 1 => logarithmic
IF (baxis EQ 1) THEN BEGIN
    logfac=30.
    READ,'Enter the multiplicative full-scale range: ',logfac
    ENDIF

; Open the single big file and associate it to variable.
OPENR,imunit,template,/GET_LUN
files=ASSOC(imunit,BYTARR(66048))

; read in each image, decode the header, rebin to
; 256x256 if necessary, scale to 8 bits and display
IF dos AND (NOT psplot) THEN WSHOW,0,ICONIC=0
ip=-1
mid=-1
nch=0
; STOP
FOR i=inf(0),inl(0) DO BEGIN
    whbad=WHERE(bad_image EQ i,nwhbad)
    IF nwhbad EQ 0 THEN BEGIN
        file=files(i)
        im=BYTE(file,512,256,256)
    ; If file was renormalized to fit 8 bits, reverse this.
        KIh=gethd(hbytes(*,i))
        IF KIh.proc.scale_factor GT 0 THEN BEGIN
            multfac=KIh.proc.scale_factor
            IF KIh.proc.base2 THEN multfac=2^multfac
            im=LONG(im)*multfac
            ENDIF

; extract exposure time if necessary
        IF selexp THEN BEGIN
            thisexp=KIh.misc.exp_scale*KIh.exp.exposure/1000.
            go_on=(ABS(thisexp-expt) LT 1.0)
            ENDIF ELSE go_on=1

        IF go_on THEN BEGIN
; extract start time of image from header
            itime=BYTARR(4)
            FOR itptr=0,3 DO itime(itptr)=KIh.misc.tm.(itptr)
            utsecs=TOTAL(itime*[3600,60,1,.01])

; extract filter temperature, determine reduced filter transmission,
; and divide image by reduced transmission to recover actual emission
; brightnesses
            fwtemp=fwtemps(i)
            imtrans=interpol(filtrans,filtemp,fwtemp)
            im=LONG(ROUND(im/imtrans(0)))

; Save the zenith brightness to the file, if specified.
        IF textout THEN PRINTF,zlun,utsecs,mean(im(inview)) $
                              ,FORMAT='(F7.1,F8.1)'

; If autoscaling is desired, find lmaxi, otherwise use topnum
            CASE 1 OF
topnum(0) EQ 0: BEGIN        ; show all except top and bottom 0.1%
                    lmaxi=pcentile(im(wmask),0.999)
                    lmaxbh=pcentile(im(wmask),0.001)
                    END
topnum(0) LT 1: BEGIN        ; adjust percentile limits for autoscaling
                    lmaxi=pcentile(im(wmask),topnum(0))
                    lmaxbh=pcentile(im(wmask),1-topnum(0))
                    END
          ELSE: BEGIN        ; fix maximum brightness to display
                    lmaxi=topnum(0)
                    lmaxbh=0
                    END
                ENDCASE
            lmaxi=lmaxi(0)
            lmaxbh=lmaxbh(0)
            im(255,98)=lmaxi        ; force image to scale properly

; If logarithmic image display selected, warn user of effects on lmaxbh
            olmaxbh=lmaxbh
            IF (baxis EQ 1) THEN lmaxbh=(lmaxi/logfac)>olmaxbh
            IF lmaxbh NE olmaxbh THEN $
                PRINT,olmaxbh,lmaxbh $
                     ,FORMAT='("Lower limit was ",F5.0," R; now "' $
                            +',F5.0," R")'

            ip=ip+1                ; index number for screen position of image

; show image orientation, once per screen or plot page
            IF (psplot AND (ip MOD 12 EQ 0)) OR $
               ((NOT psplot) AND (ip EQ 0)) THEN BEGIN

; label image display window
; keep 'X' char size
                IF ip EQ 0 THEN csz=((1.*X_charheight)/wheight) $
                                   /((1.*!D.Y_CH_SIZE)/!D.Y_VSIZE)
                XYOUTS,scl(res=psres,(wwidth-barwidth)/2) $
                      ,scl(res=psres,wheight-3*X_charheight) $
                      ,wtitle $
                      ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=2.0*csz,WIDTH=twidthn
                IF ip EQ 0 THEN BEGIN
                    twidthd=CONVERT_COORD([twidthn,0.],/NORMAL,/TO_DEVICE)
                    twidthd=twidthd(0)        ; we only need the x-component
                    ENDIF
; get user's name and display it
                IF ip EQ 0 THEN uname=username()
                XYOUTS,scl(res=psres,wwidth-X_charwidth) $
                      ,scl(res=psres,wheight-3.9*X_charheight) $
                      ,'User '+STRUPCASE(uname) $
                      ,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
; get date/time and display them
                datetime=SYSTIME()
                XYOUTS,scl(res=psres,wwidth-X_charwidth) $
                      ,scl(res=psres,wheight-2.5*X_charheight) $
                      ,STRMID(datetime,11,8)+' UT' $
                      ,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
                XYOUTS,scl(res=psres,wwidth-X_charwidth) $
                      ,scl(res=psres,wheight-1.1*X_charheight) $
                      ,STRMID(datetime,4,7)+STRMID(datetime,20,4) $
                      ,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
; create and display color bar
                IF ip EQ 0 THEN BEGIN
                    bar=BINDGEN(256,barwidth)
                    bar=REBIN(BYTSCL(bar,TOP=!D.N_COLORS-1) $
                             ,nincol*256,barwidth)
                    bar=ROTATE(bar,3)
                    ENDIF
                mtv,bar,ninrow*256,0,res=psres,neg=neg

                nflag=isnumber(KIh.misc.az,azimuth)
                IF nflag NE 0 THEN BEGIN ; if the camera azimuth is known,
                    azrad=!DTOR*azimuth        ; plot arrows showing the azimuth
                    leftend=scl(res=psres,wwidth/2.)-twidthd/2.
                    maxlen=(leftend-scl(res=psres,1.5*X_charwidth)) $
                           /scl(res=psres,(COS(azrad)+SIN(azrad))*X_charheight)
                    arlen=MIN([FIX(maxlen),4])*X_charheight
                    tailx=REPLICATE(scl(res=psres,0.5*X_charwidth+SIN(azrad)*arlen),2)
                    taily=REPLICATE(scl(res=psres,nincol*256),2)
                    headx=[tailx(0)-scl(res=psres,arlen*SIN(azrad)), $
                           tailx(1)+scl(res=psres,arlen*COS(azrad))]
                    heady=[scl(res=psres,nincol*256+arlen*COS(azrad)), $
                           scl(res=psres,nincol*256+arlen*SIN(azrad))]
                    arrow,tailx,taily,headx,heady,/device,hsize=-0.2
                    XYOUTS,headx(0)+scl(res=psres,X_charwidth) $
                          ,MIN([heady(0),scl(res=psres,wheight-X_charheight)]) $
                          ,'N',ALIGNMENT=0,/DEVICE,CHARSIZE=csz
                    XYOUTS,headx(1)+scl(res=psres,X_charwidth),heady(1),'E' $
                          ,ALIGNMENT=0.5,/DEVICE,CHARSIZE=csz
                    ENDIF ELSE BEGIN                ; otherwise show cardinal directions
                    XYOUTS,scl(3*X_charwidth,res=psres) $
                          ,scl(nincol*256+3*X_charheight,res=psres) $
                          ,'N',/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                    XYOUTS,scl(res=psres,X_charwidth) $
                          ,scl(res=psres,nincol*256+1.5*X_charheight),'W' $
                          ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                    XYOUTS,scl(res=psres,5*X_charwidth) $
                          ,scl(res=psres,nincol*256+1.5*X_charheight),'E' $
                          ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                    XYOUTS,scl(res=psres,3*X_charwidth) $
                          ,scl(res=psres,nincol*256),'S',/DEVICE $
                          ,ALIGNMENT=0.5,CHARSIZE=csz
                    ENDELSE
                ENDIF

; Display the image in the form requested
            CASE baxis OF
             0: mtv,ROTATE(((im-v)*mask)>0,5),low=lmaxbh,high=lmaxi,res=psres $
                   ,(ip MOD ninrow)*256,(nincol-1-(ip MOD ninwin)/ninrow)*256 $
                   ,neg=neg
             1: mtv,ROTATE((ALOG((im-v)>1)*mask),5),low=ALOG(lmaxbh>1) $
                   ,high=ALOG(lmaxi),res=psres $
                   ,(ip MOD ninrow)*256,(nincol-1-(ip MOD ninwin)/ninrow)*256 $
                   ,neg=neg
                ENDCASE

            IF vR THEN BEGIN
                h=HISTOGRAM(im(wmask)-v(wmask),MIN=-50,MAX=100,BIN=10)
                OPENU,vrun,'showray.vrb',/GET_LUN,/APPEND
                PRINTF,vrun,itime(0:2),h(0:14)/nwmask,TOTAL(h(0:14))/nwmask $
                           ,FORMAT='(3(I2.2,":"),16F5.2)'
                FREE_LUN,vrun

                WSET,3
                tmp=ROTATE((im-v)*mask,5)
                tmp=(tmp+(tmp LT 0)*(lmaxi-tmp))
                ; tmp(WHERE(tmp LT 0))=lmaxi
                TVSCL,tmp<lmaxi
                WSET,0

                ENDIF

; Set up to display image information
            textcol=ROUND(colfac*!D.N_COLORS)
            IF Keyword_Set(neg) THEN textcol=(!D.N_Colors-1)-textcol

; display image sequence number
            ext=STRTRIM(KIh.misc.sequence,2)
            XYOUTS,scl(res=psres,(ip MOD ninrow)*256+1) $
                  ,scl(res=psres $
                      ,(nincol-(ip MOD ninwin)/ninrow)*256-1-X_charheight) $
                  ,ext,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz $
                  ,COLOR=textcol
; and image time
            XYOUTS,scl(res=psres,(ip MOD ninrow)*256+1) $
                  ,scl(res=psres $
                      ,(nincol-1-(ip MOD ninwin)/ninrow)*256+3+X_charheight) $
                  ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz,COLOR=textcol $
                  ,STRING(itime(0:2),FORMAT='(I2.2,":",I2.2,":",I2.2)')
; and exposure time
            XYOUTS,scl(res=psres,(ip MOD ninrow)*256+1) $
                  ,scl(res=psres,(nincol-1-(ip MOD ninwin)/ninrow)*256+1) $
                  ,STRING(KIh.misc.exp_scale*KIh.exp.exposure/1000. $
                  ,FORMAT='(F5.1,1X,"s")'),/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz $
                  ,COLOR=textcol
; and maximum displayed Rayleighs
            XYOUTS,scl(res=psres,((ip MOD ninrow)+1)*256-1) $
                  ,scl(res=psres $
                      ,(nincol-(ip MOD ninwin)/ninrow)*256-1-X_charheight) $
                  ,STRING(LONG(lmaxi),FORMAT='("< ",I0.0," R")') $
                  ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz $
                  ,COLOR=textcol
; and minimum displayed Rayleighs
            XYOUTS,scl(res=psres,((ip MOD ninrow)+1)*256-1) $
                  ,scl(res=psres,(nincol-1-(ip MOD ninwin)/ninrow)*256+1) $
                  ,STRING(LONG(lmaxbh),FORMAT='("> ",I0.0," R")'),/DEVICE $
                  ,ALIGNMENT=1.0,CHARSIZE=csz $
                  ,COLOR=textcol
; and image format
            XYOUTS,scl(res=psres,((ip MOD ninrow)+1)*256-1) $
                  ,scl(res=psres $
                      ,(nincol-(ip MOD ninwin)/ninrow)*256-1-2.5*X_charheight) $
                  ,btype(baxis) $
                  ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz $
                  ,COLOR=textcol
; and solar azimuth
            solaz=sunaz(i)
            IF solaz EQ 0. THEN sunptr='' ELSE BEGIN
                IF solaz LT 180.  THEN BEGIN
                    sunptr='Sun>'
                    ENDIF ELSE BEGIN
                    sunptr='<Sun'
                    solaz=solaz-180.
                    ENDELSE
                ENDELSE
            XYOUTS,scl(res=psres,((ip MOD ninrow)+1)*256-2.5*X_charwidth-1) $
                  ,scl(res=psres $
                      ,(nincol-1-(ip MOD ninwin)/ninrow)*256+3*X_charheight+1) $
                  ,sunptr $
                  ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz $
                  ,COLOR=textcol,ORIENTATION=90-solaz

; Interpret user's latest keystroke
            PRINT,'.',FORMAT='(A,$)'
            nch=nch+1
            IF (ip+1) MOD ninwin EQ 0 THEN BEGIN
                PRINT,'|',FORMAT='(A,$)'
                IF psplot THEN ERASE
                ENDIF
            nch=nch+1

; Wrap ...|... display at the end of the line.

            IF nch GE 80 THEN BEGIN
                PRINT,FORMAT='()'
                nch=0
                ENDIF
            k=GET_KBRD(1)
            CASE STRUPCASE(k) OF
           'S': ip=ip-1
           'Q': BEGIN
                    IF dos AND (NOT psplot) THEN WSHOW,0,/ICONIC
                    PRINT,'Aborting this plot; enter hr,min,sec' $
                         +' OR image number to restart'
                    READ,resp
                    CASE nwrds(resp) OF
                     0: GOTO,DONE
                     1: BEGIN
                            sthms=resp
                            inf=FIX(getwrd(resp,0))
                            inl=nimages-1
                            PRINT,inf(0),inl(0)
                            IF psplot THEN psclose
                            take2=1B
                            GOTO,RESTART
                            END
                  ELSE: BEGIN
                            sthms=resp
                            hour=FIX(getwrd(resp,0))
                            minute=FIX(getwrd(resp,1))
                            IF nwrds(resp) GT 2  $
                                THEN second=FIX(getwrd(resp,2)) $
                                ELSE second=0
                            startut=(hour*60L+minute)*60L+second
                            IF startut GT MAX(utsec) THEN BEGIN
                                lastut=STRING(uth(nimages-1) $
                                             ,utm(nimages-1) $
                                             ,uts(nimages-1) $
                                             ,FORMAT='(3(I2.2,:,":"))')
                                PRINT,'No images after '+lastut+'!'
                                inf=0
                                inl=nimages-1
                                ENDIF ELSE BEGIN
                                inf=MIN(WHERE(utsec GE startut))
                                inl=nimages-1
                                ENDELSE
                            IF psplot THEN psclose
                            take2=1B
                            GOTO,RESTART
                            END
                        ENDCASE
                    END
           'X': GOTO,DONE
          ELSE:
                ENDCASE
            ENDIF
            EMPTY
        ENDIF
    ENDFOR

DONE:
PRINT,FORMAT='()'
IF textout THEN BEGIN
    FREE_LUN,zlun
    PRINT,'Zenith brightnesses are in '+zbroot+dd+zname
    ENDIF
IF N_ELEMENTS(imunit) GT 0 THEN FREE_LUN,imunit
IF psplot THEN psclose ELSE IF !D.WINDOW EQ 0 THEN WDELETE,0
PRINT,FORMAT='()'
END