; 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

; 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
getsig,year,month,day,mask0,mask1
CASE cam of
 0: mask=mask0
 1: mask=mask1
    ENDCASE
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)
; make list of available files of desired type
files=rfindfile(dir+template,count=nfiles)
CASE nfiles OF
    0: BEGIN
	   PRINT,'No files found; quitting.'
	   GOTO,DONE
	   END
    1: BEGIN
	   PRINT,files(0)+' available; opening file'
	   template=files(0)
	   END
    ELSE: BEGIN
           PRINT,'Files available:'
           FOR i=0,nfiles-1 DO PRINT,i,files(i),FORMAT='(I1,2X,A)'
	   READ,'Select the file to display: ',nd
	   template=files(nd)
	  END
    ENDCASE

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

nimages=rdhbytes(template)

; 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
resp=''
READ,'Make a PostScript plot file? (Y/Z/[N]) ',resp
psplot=(STRUPCASE(resp) EQ 'Y')
zbuff=(STRUPCASE(resp) EQ 'Z')
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 ELSE BEGIN
    IF zbuff THEN BEGIN
        SET_PLOT,'Z'
        DEVICE,Z_BUFFERING=0
        ENDIF
    ENDELSE

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) AND NOT zbuff 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
;   psfile=STRING(root,dd,year MOD 100,month,day,sth,stm,sts,cam,filt,expt $
;         ,FORMAT='(A,A1,3I2.2,"_",3I2.2,"_",2I1,"_",I2.2)')
;
;   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(root,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
    IF zbuff THEN DEVICE,SET_RESOLUTION=[wwidth,wheight] $
             ELSE WINDOW,0,XSIZE=wwidth,YSIZE=wheight
    TVSCL,BINDGEN(256)	; wake up the X window and set !D.N_COLORS correctly
    IF dos AND NOT zbuff THEN WSHOW,0,/ICONIC
    psres=100.		; dummy value
    ENDELSE

; 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 nct EQ '' THEN loadct,ct,/silent ELSE BEGIN
    IF isnumber(nct) THEN BEGIN
        nct=FLOAT(nct)
        ct=FIX(ABS(nct))
        IF nct LT 0 THEN BEGIN
            loadct,ct,/silent
            READ,'Enter the gamma correction factor: ',gamma
            stretch,0,!D.N_COLORS-1,gamma
            ENDIF
        ENDIF ELSE xloadct
    ENDELSE
whitetop

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

; find the zenith pixels (z.a. < 2 deg) in these images

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
zpix=WHERE(za LE 2.0,nzpix)
IF nzpix EQ 0 THEN STOP

; 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
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=ROUND(im/imtrans(0))
            PRINT,utsecs,mean(im(zpix))

; 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

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

; 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=ROUND(colfac*!D.N_COLORS)
; 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=ROUND(colfac*!D.N_COLORS) $
;                  ,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=ROUND(colfac*!D.N_COLORS)
; 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=ROUND(colfac*!D.N_COLORS)
; 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=ROUND(colfac*!D.N_COLORS)
; 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=ROUND(colfac*!D.N_COLORS)
; 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=ROUND(colfac*!D.N_COLORS),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)
            k=''
            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:
IF N_ELEMENTS(imunit) GT 0 THEN FREE_LUN,imunit
IF psplot THEN psclose ELSE IF NOT zbuff THEN WDELETE,0
PRINT,FORMAT='()'
END