; MODSEQ.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 1993/5/11 from SHOW.PRO to display images from both cameras
; in the order they were acquired.
;
; create byte-valued input parameters
cam=0B & filt=cam
day=0 & month=day
mask=BYTARR(256,256)
; load various string-valued arrays needed
initialize,filters,modes,mpaths,firstchars
s=''
!ORDER=1
; set colour table 
loadct,0,/silent
; uncover screen where text is to be printed
IF !D.NAME EQ 'X' THEN WDELETE,0 

; Warn user that individual data files are needed.
PRINT,'This program reads individual raw data files, rather than .ray or .bgs files.'
PRINT,'If you don''t have these files on disk, stop now and load them.'
STOP
Print,'NOTE: This modified version of SEQ.PRO displays only images from'
Print,'Camera 0.  It can be used to create a PostScript plot if the'
Print,'PostScript device has been selected and initialized before'
Print,'running this program.'

; prompt user and get input
PRINT,'Enter the starting year, month, day, UT hour, minute, and second desired'
READ,year,month,day,hr,min,sec
resp=''
Read,'Subtract dark frame before displaying images? ([Y]/N) ',resp
nodk=StrUpCase(resp) EQ 'N'
IF !D.NAME NE 'PS' THEN BEGIN
    Read,'Make .GIF file of each filled image frame? (Y/[N]) ',resp
    makegif=StrUpCase(resp) EQ 'Y'
    IF makegif THEN BEGIN
        gif_ctr=0
        gif_root='C:\DATA\'
        ENDIF
    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(year,months(month-1),day $
;	     ,FORMAT='("Polar Camera",2X,I4,1X,A3,1X,I0.2)')
wtitle="Polar Camera  "+ymd2date(year,month,day,format='y$ n$ d$')
PRINT,wtitle

; define mask for images, with fixed horizon at 150 superpixels 
prad=150
d256=SHIFT(dist(256),128,128) & base=WHERE(d256 GT 150)
wprad=WHERE(d256 LE prad) & mask(wprad)=1

; make lists of available files of desired type
rdilist,0,year,month,day,hr,min,sec,n0,t0,files0
; rdilist,1,year,month,day,hr,min,sec,n1,t1,files1

; identify and remove setup images (25 ms exposure time)
wncal=WHERE(STRPOS(files0,'cal') EQ -1,n0)
;wn25ms=WHERE(STRMID(files0,STRLEN(files0)-7,3) NE '000',n0)
t=t0(wncal)
files=files0(wncal)
; wncal=WHERE(STRPOS(files1,'cal') EQ -1,n1)
;wn25ms=WHERE(STRMID(files1,STRLEN(files1)-7,3) NE '000',n1)
; t1=t1(wncal)
; files1=files1(wncal)

; zip lists together
; IF (n0 EQ n1) THEN BEGIN	; equal numbers in both lists
;     IF (MIN(ABS(t0-t1)) LT .5) THEN BEGIN
;         t=TRANSPOSE([TRANSPOSE(t0),TRANSPOSE(t1)])
;         files=TRANSPOSE([TRANSPOSE(files0),TRANSPOSE(files1)])
;     ENDIF
; ENDIF ELSE BEGIN		; lists have different lengths
;     t=FLTARR(MAX([n0,n1]),2)
;     files=STRARR(MAX([n0,n1]),2)
;     IF n0 GT n1 THEN BEGIN	; more from Camera 0
;         t(*,0)=t0
;         files(*,0)=files0
;         FOR i=0,n1-1 DO BEGIN	; find best match for each Cam 1 image
;             d=MIN(ABS(t0-t1(i)),p)
;             IF d LT 0.5 THEN BEGIN
;                 t(p,1)=t1(i)
;                 files(p,1)=files1(i)
;             ENDIF
;         ENDFOR
;     ENDIF ELSE BEGIN		; more from Cam 1
;         t(*,1)=t1
;         files(*,1)=files1
;         FOR i=0,n0-1 DO BEGIN   ; find best match for each Cam 0 image
;             d=MIN(ABS(t1-t0(i)),p)
;             IF d LT 0.5 THEN BEGIN
;                 t(p,0)=t0(i)
;                 files(p,0)=files0(i)
;             ENDIF
;         ENDFOR
;     ENDELSE
; ENDELSE
ext_list=extension(files)

; read in or generate required dark frames
dkt=0
dkim=INTARR(1,256)
dkpair=INTARR(256,256)
FOR i=0,n0-1 DO BEGIN
    IF STRPOS(files(i),'_No') EQ -1 THEN BEGIN
        exp0=FIX(STRMID(files(i),STRLEN(files(i))-7,3))
        we=WHERE(exp0 EQ dkt,nwe)
        IF nwe EQ 0 THEN BEGIN	; haven't got this exposure time yet
            IF NOT nodk THEN BEGIN
                rdmeandk,0,exp0,dk0
;                 rdmeandk,1,exp0,dk1
                ENDIF ELSE dk0=INTARR(256,256)
            dkpair(*,0:255)=dk0
;             dkpair(*,256:511)=dk1
            dkt=[dkt,exp0]
            dkim=[dkim,dkpair]
        ENDIF
    ENDIF
ENDFOR
; FOR i=0,n1-1 DO BEGIN
;     IF STRPOS(files1(i),'_No') EQ -1 THEN BEGIN
;         exp1=FIX(STRMID(files1(i),STRLEN(files1(i))-7,3))
;         we=WHERE(exp1 EQ dkt,nwe)
;         IF nwe EQ 0 THEN BEGIN	; haven't got this exposure time yet
;             rdmeandk,0,exp1,dk0
;             rdmeandk,1,exp1,dk1
;             dkpair(*,0:255)=dk0
;             dkpair(*,256:511)=dk1
;             dkt=[dkt,exp1]
;             dkim=[dkim,dkpair]
;         ENDIF
;     ENDIF
; ENDFOR
dkt=dkt(1:*)		; drop dummy entry at start of array
dkim=dkim(1:*,*)	; drop dummy column ...

; create and label image display window
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
IF (!D.NAME EQ 'X') OR (!D.NAME EQ 'WIN') THEN $
    WINDOW,0,XSIZE=wwidth,YSIZE=wheight $
    ELSE DEVICE,/COLOR,BITS_PER_PIXEL=8	; PostScript
csz=((1.*X_charheight)/wheight)/((1.*!D.Y_CH_SIZE)/!D.Y_VSIZE)	; keep 'X' char size
XYOUTS,scl(res=120.,(wwidth-barwidth)/2) $
      ,scl(res=120.,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
IF (STRLOWCASE(STRMID(!VERSION.OS,0,3)) NE 'win') $
    THEN SPAWN,'whoami',username $
    ELSE username='CNSR3'
XYOUTS,scl(res=120.,wwidth-X_charwidth),scl(res=120.,wheight-3.9*X_charheight) $
      ,'User '+STRUPCASE(username),ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
; get date/time and display them
dt_tm_brk,systime(),date,time
XYOUTS,scl(res=120.,wwidth-X_charwidth),scl(res=120.,wheight-2.5*X_charheight) $
      ,time+' UT',ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
XYOUTS,scl(res=120.,wwidth-X_charwidth),scl(res=120.,wheight-1.1*X_charheight) $
      ,date,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz
; create and display color bar
bar=REBIN(BINDGEN(256,barwidth),nincol*256,barwidth)
bar=ROTATE(bar,3)
mtv,BYTSCL(bar,TOP=!D.N_COLORS-1),ninrow*256,0,res=120.

; 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
itime=BYTARR(4)
FOR i=0,n0-1 DO BEGIN
    ti=t(i)
;     p=WHERE(ABS(t(*,1)-ti),nm)
;     IF nm GT -1 THEN p=[i,p(0)] ELSE p=[i,i]
;     FOR j=0,1 DO BEGIN
        fnlen=STRLEN(files(i))
        IF fnlen GT 0 THEN BEGIN
            IF StrPos(files(i),'__No') EQ -1 THEN BEGIN
                rdkimg,files(i),h,im
                KIh=gethd(h)
                FOR k=0,3 DO itime(k)=KIh.misc.tm.(k)
                ENDIF ELSE BEGIN
                mk,im           ; make a 256x256 array of 0B
                ENDELSE
; if image isn't 256x256, rebin it
            im=REBIN(im,256,256,/SAMPLE)
; ensure counts are actually summed if a 512x512 image is compressed
            cfac=512./(256*KIh.exp.sbin)
            rfac=512./(256*KIh.exp.pbin)
            im=im*cfac*rfac
; subtract appropriate dark frame
            expt=FIX(KIh.misc.exp_scale*KIh.exp.exposure/1000.)
            dp=WHERE(dkt EQ expt)
; if this is a dark frame, store it in the dkim array
	    IF (STRPOS(files(i),'dark') NE -1) AND (NOT nodk) THEN $
		dkim(256*dp(0),0)=MEDIAN(im,3)
            im=FIX((im-dkim(256*dp(0):256*dp(0)+255,0:255))>0)
; identify useful min and max values for byte-scaling
;            ih=HISTOGRAM(im(wprad),MIN=0,BIN=1)	; histogram of exposed image
;            cih=ih	; next line computes cumulative pixel value pdf
;            FOR jj=1,N_ELEMENTS(cih)-1 DO cih(jj)=cih(jj)+cih(jj-1)
;            cih=cih/FLOAT(MAX(cih))		; normalize to unity
;            lmaxi=MAX(WHERE(cih LE .999))	; top 0.1% of pixels saturate
            lmaxi=pcentile(im(wprad),0.999)
; if viewing sky image, take dark 'skirt' as min value
            exptype=STRMID(files(i),fnlen-12,1)
            shut=(exptype EQ 'b') OR (exptype EQ 'd')	; bias or dark frame
            IF NOT shut THEN BEGIN
                bh=HISTOGRAM(im(base),MIN=0,BIN=1)
                maxbh=MAX(bh,lmaxbh)		; bottom = most prob. skirt value
            ENDIF ELSE BEGIN
; otherwise take lowest level in dark frame
                lmini=pcentile(im(wprad),0.001)
                lmaxbh=lmini	; rename this parameter for compatibility below
            ENDELSE
        ENDIF ELSE BEGIN
; if no file name, generate an identifiably fake image
            im=BYTARR(256,256)
            im(RANDOMU(s,100))=1
            itime=BYTARR(4)
            expt=-1
            lmaxi=1
            lmaxbh=0
        ENDELSE

        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=120.,wwidth/2.)-twidthd/2.
                maxlen=(leftend-scl(res=120.,1.5*X_charwidth)) $
                       /scl(res=120.,(COS(azrad)+SIN(azrad))*X_charheight)
                arlen=MIN([FIX(maxlen),4])*X_charheight
                tailx=REPLICATE(scl(res=120.,0.5*X_charwidth+SIN(azrad)*arlen),2)
                taily=REPLICATE(scl(res=120.,nincol*256),2)
                headx=[tailx(0)-scl(res=120.,arlen*SIN(azrad)), $
                       tailx(1)+scl(res=120.,arlen*COS(azrad))]
                heady=[scl(res=120.,nincol*256+arlen*COS(azrad)), $
                       scl(res=120.,nincol*256+arlen*SIN(azrad))]
                arrow,tailx,taily,headx,heady,/device,hsize=-0.2
                XYOUTS,tailx(0),heady(0)-scl(res=120.,X_charheight),'N' $
                       ,ALIGNMENT=0,/DEVICE,CHARSIZE=csz
                XYOUTS,headx(1)+scl(res=120.,X_charwidth),heady(1),'E' $
                       ,ALIGNMENT=0.5,/DEVICE,CHARSIZE=csz
            ENDIF ELSE BEGIN		; otherwise show cardinal directions
                XYOUTS,scl(3*X_charwidth,res=120.) $
                      ,scl(nincol*256+3*X_charheight,res=120.) $
                      ,'N',/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                XYOUTS,scl(res=120.,X_charwidth) $
                      ,scl(res=120.,nincol*256+1.5*X_charheight),'W' $
                      ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                XYOUTS,scl(res=120.,5*X_charwidth) $
                      ,scl(res=120.,nincol*256+1.5*X_charheight),'E' $
                      ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz
                XYOUTS,scl(res=120.,3*X_charwidth) $
                      ,scl(res=120.,nincol*256),'S',/DEVICE $
                      ,ALIGNMENT=0.5,CHARSIZE=csz
            ENDELSE
        ENDIF
; display the image
        mtv,ROTATE(im*mask,5),low=lmaxbh,high=lmaxi,res=120. $
           ,(ip MOD ninrow)*256,(nincol-1-(ip MOD ninwin)/ninrow)*256
; display image camera, filter, and sequence numbers
        IF fnlen GT 0 $
            THEN camfilt=STRMID(files(i),fnlen-9,2) $
            ELSE camfilt=STRMID(files(i),fnlen-9,2)
        XYOUTS,scl(res=120.,(ip MOD ninrow)*256+1) $
              ,scl(res=120.,(nincol-(ip MOD ninwin)/ninrow)*256-1-X_charheight) $
              ,camfilt+':'+ext_list(i),/DEVICE,ALIGNMENT=0.0,CHARSIZE=2*csz $
              ,COLOR=!D.N_COLORS
; and image time
        XYOUTS,scl(res=120.,(ip MOD ninrow)*256+1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*256+3+X_charheight) $
              ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=2*csz,COLOR=!D.N_COLORS $
              ,STRING(itime(0:2),FORMAT='(I2.2,":",I2.2,":",I2.2)')
; and exposure time
        XYOUTS,scl(res=120.,(ip MOD ninrow)*256+1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*256+1) $
              ,STRING(expt,FORMAT='(F5.1,1X,"s")'),/DEVICE,ALIGNMENT=0.0 $
              ,CHARSIZE=2*csz,COLOR=!D.N_COLORS
; and maximum displayed DN
        XYOUTS,scl(res=120.,((ip MOD ninrow)+1)*256-1) $
              ,scl(res=120.,(nincol-(ip MOD ninwin)/ninrow)*256-1-X_charheight) $
              ,STRING(FIX(lmaxi),FORMAT='("< ",I0.0," DN")') $
              ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=2*csz,COLOR=!D.N_COLORS
; and minimum displayed DN
        XYOUTS,scl(res=120.,((ip MOD ninrow)+1)*256-1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*256+1) $
              ,STRING(FIX(lmaxbh),FORMAT='("> ",I0.0," DN")'),/DEVICE $
              ,ALIGNMENT=1.0,CHARSIZE=2*csz,COLOR=!D.N_COLORS
        PRINT,'.',FORMAT='(A,$)'
        IF (ip+1) MOD ninwin EQ 0 THEN BEGIN
            PRINT,'|',FORMAT='(A,$)'
            IF makegif THEN BEGIN
                !Order=1-!Order
                gif=TVRD()
                !Order=1-!Order
                Help,gif
                gif_name=String(gif_root,year MOD 100,month,day,gif_ctr $
                               ,Format='(A,4I2.2,".GIF")')
                Help,gif_name
                Write_GIF,gif_name,gif
                gif_ctr=gif_ctr+1
                ENDIF
            ENDIF
;    ENDFOR
    IF GET_KBRD(1) EQ 'q' THEN GOTO,DONE
ENDFOR
DONE:
WDelete,!D.Window
PRINT,''
END