; 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
dir=saveroot+dd
;dir='/jasper/cnsr3_data1/save/'
; define file name for desired file
template=STRING(year MOD 100,month,day,cam,filt,FORMAT='(3I2.2,2I1,".*")')
; make list of available files of desired type
files=FINDFILE(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/[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

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
    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)')
    IF KEYWORD_SET(colour) THEN suf='.cps' ELSE suf='.ps'
    psfile=psfile+suf
    psopen,psfile,/landscape,color=colour
    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
    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 NE '' THEN ct=FIX(nct)
loadct,ct,/silent
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 LT 0 THEN topnum=LONG(topnum)+(2L^16)

; 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
ip=-1
mid=-1
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.
        IF file(290) GT 0 THEN BEGIN
            IF file(291) EQ 1 THEN im=LONG(im)*(2^file(290))
	    IF file(291) EQ 0 THEN im=LONG(im)*file(290)
	    ENDIF
        KIh=gethd(hbytes(*,i))

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

; If autoscaling is desired, find lmaxi, otherwise use topnum
	    CASE 1 OF
   topnum 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 LT 1: BEGIN	; adjust percentile limits for autoscaling
		    lmaxi=pcentile(im(wmask),topnum)
		    lmaxbh=pcentile(im(wmask),1-topnum)
		    END
          ELSE: BEGIN	; fix maximum brightness to display
		    lmaxi=topnum
		    lmaxbh=0
		    END
		ENDCASE
	    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 SPAWN,'whoami',username
                XYOUTS,scl(res=psres,wwidth-X_charwidth) $
		      ,scl(res=psres,wheight-3.9*X_charheight) $
                      ,'User '+STRUPCASE(username) $
		      ,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
; get date/time and display them
                SPAWN,'date -u',datetime
                XYOUTS,scl(res=psres,wwidth-X_charwidth) $
		      ,scl(res=psres,wheight-2.5*X_charheight) $
                      ,STRMID(datetime,11,9)+' 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,24,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,$)'
            IF (ip+1) MOD ninwin EQ 0 THEN BEGIN
		PRINT,'|',FORMAT='(A,$)'
		IF psplot THEN ERASE
		ENDIF
            k=GET_KBRD(1)
            CASE STRUPCASE(k) OF
           'S': ip=ip-1
           'Q': BEGIN
                    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 WDELETE,0
PRINT,FORMAT='()'
END