; NEWSHOWRAY.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).
;
; create byte-valued input parameters
cam=0B & filt=cam & day=cam & month=cam
; load various string-valued arrays needed
initialize,filters,months,days,modes,mpaths,firstchars
s=''
psplot=0B
!ORDER=1
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
PRINT,'Enter the camera, filter, year, month, and day desired'
READ,cam,filt,year,month,day
hour=0
minute=0
second=0
PRINT,'Enter the hour, minute, and second to start (default 00:00:00 UT)'
sthms=''
READ,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),year,months(month-1), $
	      day,FORMAT='("Polar Camera",2X,A,2X,I4,1X,A3,1X,I0.2)')

; Display only those pixels viewing the sky
sigfile='/jasper/cnsr3_data1/sig'+STRTRIM(fix(cam),2)+'.EU'
OPENR,sigunit,sigfile,/GET_LUN
nsigs=0L
READU,sigunit,nsigs
sigpix=LONARR(nsigs)
READU,sigunit,sigpix
CLOSE,sigunit
FREE_LUN,sigunit

mk,mask
mask(sigpix)=1

; define path to desired images
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)'
	   PRINT,'Select the file to display'
	   READ,nd
	   template=files(nd)
	  END
    ENDCASE

; Open the file as a series of small (header-size) units and get
; the header information
OPENR,hbunit,template,/GET_LUN
hbfile=ASSOC(hbunit,BYTARR(512))
hbstat=FSTAT(hbunit)
hbbytes=hbstat.SIZE
nimages=hbbytes/131584L
PRINT,'There are '+STRTRIM(nimages,2)+' images'
hbytes=BYTARR(512,nimages)
utsec=LONARR(nimages)
exptimes=FLTARR(nimages)
FOR i=0L,nimages-1 DO BEGIN
    hbytes(0,i)=hbfile(129*i)
    header=gethd(hbytes(*,i))
    utsec(i)=3600L*header.misc.tm.hour $
	     + 60L*header.misc.tm.minute $
	     +  1L*header.misc.tm.second
    exptimes(i)=header.misc.exp_scale*header.exp.exposure/1000.
    IF (i MOD 20 EQ 0) THEN PRINT,i,FORMAT='(I4,$)'
    IF (i GT 0) AND (i MOD 400 EQ 0) THEN PRINT,FORMAT='(/)'
    ENDFOR
PRINT,''
CLOSE,hbunit
FREE_LUN,hbunit

; 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=''
PRINT,'Make a PostScript plot file? (Y/[N])'
READ,resp
psplot=(STRUPCASE(resp) EQ 'Y')

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")'
    PRINT,'Enter the exposure time of interest (<Return> for all)'
    texp=''
    READ,texp
    selexp=isnumber(texp,expt) 
    ENDIF ELSE BEGIN
    selexp=0
    expt=FIX(ROUND(mnex))
    ENDELSE

IF sthms EQ '' THEN BEGIN
    PRINT,'Enter the range of image numbers to display (<Return> for all images)'
    instring='' & READ,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

; 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(year MOD 100,month,day,sth,stm,sts,cam,filt,expt $
          ,FORMAT='("/jasper/cnsr3_data1/",3I2.2,".",3I2.2,".",2I1,".",I2.2,".ps")')
    psopen,psfile,/landscape,/col
    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
    0: ct=0
    2: ct=3
    5: ct=0
    7: ct=1
    8: ct=8
    ELSE: ct=0
    ENDCASE
loadct,ct,/silent

; Prompt user for fixed saturation brightness
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.'
tops=''
READ,tops
flag=isnumber(tops,topnum)	; topnum=0 if tops doesn't represent a number

; label image display window
csz=((1.*X_charheight)/wheight)/((1.*!D.Y_CH_SIZE)/!D.Y_VSIZE)	; keep 'X' char size
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
twidthd=CONVERT_COORD([twidthn,0.],/NORMAL,/TO_DEVICE)
twidthd=twidthd(0)			; we only need the x-component
; get user's name and display it
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
bar=REBIN(BINDGEN(256,barwidth),nincol*256,barwidth)
bar=ROTATE(bar,3)
mtv,bar,ninrow*256,0,res=psres

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

; read in each image, decode the header, subtract dark frame, rebin to 
; 256x256 if necessary, scale to 8 bits and display
ip=-1
mid=-1
FOR i=inf(0),inl(0) DO BEGIN
    file=files(i)
    im=BYTE(file,512,256,256)
    ; If file was renormalized to fit 16 bits, reverse this.
    IF file(290) GT 0 THEN im=LONG(im)*(2^file(290))
    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])

; If autoscaling is desired, find lmaxi, otherwise use topnum
        IF topnum EQ 0 $
; Is there any useful information left in the image?
; identify useful min and max values for byte-scaling
	    THEN lmaxi=percentile(im,0.999) $
            ELSE lmaxi=topnum
        lmaxbh=percentile(im,0.001)

        ip=ip+1		; index number for screen position of image
; show image orientation, one time only
        IF ip EQ 0 THEN BEGIN
            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
        mtv,ROTATE(im*mask,5),low=lmaxbh,high=lmaxi,res=psres $
           ,(ip MOD ninrow)*256,(nincol-1-(ip MOD ninwin)/ninrow)*256

; 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=0.8*!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=0.8*!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=0.8*!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(FIX(lmaxi),FORMAT='("< ",I0.0," R")') $
              ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz,COLOR=0.8*!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(FIX(lmaxbh),FORMAT='("> ",I0.0," R")'),/DEVICE $
              ,ALIGNMENT=1.0,CHARSIZE=csz,COLOR=0.8*!D.N_COLORS
        PRINT,'.',FORMAT='(A,$)'
        IF (ip+1) MOD ninwin EQ 0 THEN PRINT,'|',FORMAT='(A,$)'
        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
			        inf=FIX(getwrd(resp,0))
			        inl=nimages-1
				IF psplot THEN psclose
			        GOTO,RESTART
			        END
		         ELSE: BEGIN
			           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
			           inf=MIN(WHERE(utsec GE startut))
			           inl=nimages-1
				   IF psplot THEN psclose
			           GOTO,RESTART
			           END
		         ENDCASE
		     END
            'X': GOTO,DONE
	    ELSE:
	    ENDCASE
        ENDIF
    ENDFOR

DONE:
IF psplot THEN psclose
PRINT,''
END