;+
; NAME: TRAY2BGS
;	
; PURPOSE:
;	To subtract images from a given interval from the minimum
;	image produced over that interval, or from some other
;	minimum image (i.e., another camera-filter combination)
;	over that same interval and store it in a file.
;
; CATEGORY:
;	
; CALLING SEQUENCE:
;	TRAY2BGS,YEAR,MONTH,DAY,HH,MM,SS,SPAN,CN,FN $
;	        [,EXP=EXP] [,CAM=CAM] [,FILT=FILT] [,VRHIJN=VRHIJN]	
;
; INPUTS:
;	YEAR, MONTH, DAY:
;		The year, month, day that the images of interest
;		were created.
;	HH, MM, SS:    
;		The hour, minute and second in UT, that marks the
;		beginning of the interval.
;	SPAN:   The time, in hours, that the interval extends from
;		the start time.
;	CN:	The camera number.
;	FN:	The filter number.
;
; KEYWORD PARAMETERS:
;	EXP:    The exposure time of the images to be used. If it
;		is not set then default will be 60 seconds.
;	CAM:    The camera number of the minimum image to be used.
;		If it is not set, then default will be the number
;		specified by CN.
;	FILT:   The filter number of the minimum image to be used.
;		If it is not set, then default will be the number
;		specified by FN.
;	VRHIJN: A van Rhijn layer will be used. If not set, the
;		default will be the normal minimum image.
;
; OUTPUTS:
;	None.
;
; COMMON BLOCKS:
;	HEADER,FNAME,HBYTES,NIMAGES,UTSEC,EXPTIMES,KEEP	
;	FNAME:  Stores the name of the .ray file that has been
;		last used. It will only change when a new .ray file
;		has been opened.
;	NIMAGES: The number of images in the .ray file.
;	HBYTES: Stored the 512 header bytes of all the images in
;		the .ray file.
;	UTSEC:  Stores the startime in UT seconds of all the images
;		in the .ray file.
;	EXPTIMES: Stores all the exposure time of the images in
;		the .ray file.
;	KEEP:   Holds the indexes into the HBYTES, EXPTIMES, and
;		UTSEC, arrays of all the images with exposure times
;		specified by the user.
;
; SIDE EFFECTS:
;	A file is made in the ~cnsr3/save directory. The file has
;	the form yearmonthdaycamfilt.bgs.
;
; RESTRICTIONS:
;	None known.
;
; PROCEDURE:
;	.... 
;
; MODIFICATION HISTORY:
;	Written August 1994, by T A Oliynyk	
;-

PRO tray2bgs,year,month,day,hh,mm,ss,span,cn,fn $
   ,exp=exp,cam=cam,filt=filt,vrhijn=vrhijn

IF N_PARAMS() LT 9 THEN BEGIN
    doc_library,'tray2bgs'
    RETURN
    ENDIF

COMMON header,fname,hbytes,nimages,utsec,exptimes,keep

site=PoCasite(year,month,day)

sy=STRING(year MOD 100,FORMAT='(I2.2)')
sm=STRING(month,FORMAT='(I2.2)')
sd=STRING(day,FORMAT='(I2.2)')
date=sy+'/'+sm+'/'+sd
scn=STRING(cn,FORMAT='(I1)')
sfn=STRING(fn,FORMAT='(I1)')

path='/jasper/cnsr3_data1/save/'
rayfile=path+sy+sm+sd+scn+sfn+'.ray'

f=FINDFILE(rayfile,COUNT=existr)

IF existr NE 0 THEN BEGIN 

    IF NOT KEYWORD_SET(exp) THEN exp=60. ELSE exp=FLOAT(exp)
    IF (N_ELEMENTS(fname) EQ 0)  THEN fname=''
    IF (N_ELEMENTS(utsec) EQ 0) OR (fname NE rayfile) THEN BEGIN

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

        fname=rayfile
        OPENR,hbunit,rayfile,/GET_LUN
        hbfile=ASSOC(hbunit,BYTARR(512))
        hbstat=FSTAT(hbunit)
        hbbytes=hbstat.SIZE
        nimages=hbbytes/66048L
        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 GT 0) AND (i MOD 20 EQ 0) THEN PRINT,i,FORMAT='(I4,$)'
            IF (i GT 0) AND (i MOD 400 EQ 0) THEN PRINT,FORMAT='()'
            ENDFOR	; (i=0L,nimages-1)
        PRINT,FORMAT='()'
        FREE_LUN,hbunit
        ENDIF		; (N_ELEMENTS(utsec) EQ 0)
  
    keep=WHERE(exptimes eq exp) ; keep only the images with exposure
    				; times less than exp.

    IF keep(0) NE -1 THEN BEGIN
        bpath='/jasper/cnsr3_data1/bkgrd/'

; if keywords cam and filt are set, then use them to identify the
; proper .bkg file.

        IF (N_ELEMENTS(cam) GT 0) OR (N_ELEMENTS(filt) GT 0) THEN BEGIN
            IF N_ELEMENTS(cam) GT 0 THEN cm=cam ELSE cm=cn
            IF N_ELEMENTS(filt) GT 0 THEN ft=filt ELSE ft=fn
            factor=tpbfactor(cn,fn)*(exp/60.)	; find the filter passband 
            ENDIF ELSE BEGIN
            cm=cn
            ft=fn
            factor=1. 
            ENDELSE	; (N_ELEMENTS(cam) GT 0) OR ...

; Define 'background' cam/filt strings.
        bcm=STRING(cm,FORMAT='(I1)')
        bft=STRING(ft,FORMAT='(I1)')
    
        bkgfile=bpath+sy+sm+sd+bcm+bft+'.bkg'
        f=FINDFILE(bkgfile,COUNT=existbkg)  ; find appropriate .bkg file.
     
        IF existbkg NE 0 THEN BEGIN
	    get_bkghds,year,month,day,cm,ft,hd ; get header information from .bkg file.

; find minimum images produced from the interval specified by the user.
      
            mspan=LONG(ROUND(span*3600))
            ind=WHERE(ABS(hd.utsec.int-mspan) LT 1 )
      
            IF ind(0) NE -1 THEN BEGIN
         
; find images from the .ray file that were produced between the times
; specified by the user.

	        butsec=3600L*hh+60L*mm+1L*ss            ; start time
	        eutsec=ROUND(butsec+FLOAT(span)*3600L)  ; stop time
	        rimgs=WHERE((utsec GE butsec) AND (utsec LT eutsec))
		; index of images between the start and stop times.

                IF rimgs(0) NE -1 THEN BEGIN
         
; find the minimum image from the .bkg file that was produced from
; the same interval as the identifed images from the .ray file.

	            kbkg=WHERE(ABS(hd(1:*).utsec.sut-butsec) LT mspan/2)+1      
                    IF kbkg(0) NE 0 THEN BEGIN 
	                com=where_array(ind,kbkg)
	                kbkg=kbkg([com(uniq(com))])

		        kbkg=FIX(MEDIAN(kbkg))   ; index of the
		        ; minimum image to be used for subtraction

; take the background image indexed by kbkg and subtract it from
; every image in the .ray file indexed by rimgs.

		        OPENR,runit,rayfile,/GET_LUN
		        OPENR,bunit,bkgfile,/GET_LUN

		        rayfiles=ASSOC(runit,BYTARR(66048)) ; .ray
		       				; images and headers
		        bkgfiles=ASSOC(bunit,BYTARR(262168)); .bkg
		       				; images and headers

		        bkgim=FLOAT(bkgfiles(kbkg),24,256,256)
		        ; minimum image to be used for subtraction

		        IF KEYWORD_SET(vrhijn) THEN BEGIN
			    vr=FLOAT(bkgfiles(0),24,256,256) ; van
			    	; Rhijn layer stored in .bkg file
			    bkgim=bkgim/vr	; van Rhijn corrected
			    			; .bkg minimum image

; find the intensity of the van Rhijn layer (for the minimum image)

	                    pcent=pcentile(bkgim,0.999)
	                    hmax=CEIL(MAX(bkgim<pcent))
	                    hmin=FLOOR(MIN(bkgim<pcent))
	                    hist=HISTOGRAM(bkgim,MAX=hmax,MIN=hmin)
	                    nbins=N_ELEMENTS(hist)
	                    intensity=INDGEN(nbins)+hmin
	                    kpoint=WHERE(intensity GT 0)
	                    hist=hist(kpoint)
	                    intensity=intensity(kpoint)
	                    ghist=GAUSSFIT(intensity,hist)
			    ; fit a gaussian to the histogram
	                    z=MAX(WHERE(ghist EQ MAX(ghist)))
	     
	                    bright=BYTE(intensity(z))
			    ; brightness of the layer
			    bkgim=bright*vr     ; van Rhijn layer
			    ; corresponding to the minimum image.

			    ENDIF ELSE bright=255B ; (KWS(vrhijn))
			    ; set layer brightness to 255, to
			    ; identify that a van Rhijn layer was
			    ; not used.

; If background and feature images are from different cameras, we
; have to sample the background image at the right look directions.

			IF cn NE cm THEN BEGIN
			    colf=REBIN(INDGEN(256),256,256)
			    rowf=TRANSPOSE(colf)
			    cr2az,year,month,day,cn,fn,colf,rowf $
				 ,azf,zaf,/degrees,bin=2,site=site
                            az2cr,year,month,day,cm,ft,azf,zaf $
				 ,colb,rowb,/degrees,bin=2,site=site
			    colb=colb>0<255
			    rowb=rowb>0<255
			    ; WINDOW,/FREE,XSIZE=512,YSIZE=256
			    ; twin=!D.WINDOW
			    ; TVSCL,bkgim<100,0
                            bkgim=bkgim(colb,rowb)
			    ; TVSCL,bkgim<100,1
			    ; xxx=GET_KBRD(1)
			    ; WDELETE,twin
                            ENDIF	; (cn NE cm)

; open appropriate .bgs file if it exist. If it does not exist then
; create it.

	                bgsfile=path+sy+sm+sd+scn+sfn+'.bgs'
	                f=FINDFILE(bgsfile,COUNT=existbgs)
	                IF existbgs EQ 0 THEN BEGIN
                            OPENW,sunit,bgsfile,/GET_LUN
                            ENDIF ELSE BEGIN
; get header bytes from .bgs file (write procedure to do this)
	                    OPENU,sunit,bgsfile,/APPEND,/GET_LUN
                            ENDELSE	; (existbgs EQ 0)
 
	                mspan=BYTE(ROUND(60L*span))

	                FOR i=0, N_ELEMENTS(rimgs)-1 DO BEGIN 
                            file=rayfiles(rimgs(i))
	                    rayim=BYTE(file,512,256,256)

; If an image from the .ray file was renormalized to fit 8 bits,
; reverse this.

	                    IF file(290) GT 0 THEN BEGIN
	                        IF file(291) EQ 1 THEN $
				    rayim=LONG(rayim)*(2L^file(290))
	                        IF file(291) EQ 0 THEN $
				    rayim=LONG(rayim)*file(290)
	                        ENDIF

; subract minimum background image from the images from the .ray file. 

	                    bgsim=rayim-factor*bkgim
                            bgsim=bgsim<pcentile(bgsim,0.995)>0

	                    bgsmax=MAX(bgsim)
	                    WRITEU,sunit,file(0:289)

; normalize the images to be stored in the .bgs file to 8 bits

	                    IF bgsmax LE 255L*255L THEN BEGIN
	                        scla=255L*INDGEN(255)
                                scale=BYTE(MAX(WHERE(scla LE bgsmax)))+1B
	                        bgsim=BYTE(bgsim/FLOAT(scale)) 
	                        WRITEU,sunit,scale,0B 
	                        ENDIF ELSE BEGIN
                                scla=255L*2^LINDGEN(20)
	                        scale=BYTE(MAX(WHERE(scla LE bgsmax)))+1B
	                        bgsim=BYTE(bgsim/FLOAT(2L^scale))
	                        WRITEU,sunit,scale,1B
                                ENDELSE	; (bgsmax LE 255L*255L)

; write images and header bytes to the .bgs file

                            WRITEU,sunit,BYTE(cn),BYTE(fn) $
				  ,BYTE(hh),BYTE(mm),BYTE(ss) $
		                  ,mspan,bright,file(299:511),bgsim  
	                    ENDFOR	; (i=0, N_ELEMENTS(rimgs)-1)

                        FREE_LUN,runit
	                FREE_LUN,bunit
	                FREE_LUN,sunit

                        ENDIF ELSE BEGIN
	                PRINT,bcm,bft,date $
	    		     ,FORMAT='("The .bkg file for camera ",A,' $
				    +'",filter ",A,", on date ",A,' $
				    +'", does not contain a background image")'
	                PRINT,hh,mm,ss,span $
	                     ,FORMAT='("corresponding to the interval ' $
				    +'starting at",1X,I2.2,":",I2.2,":",' $
				    +'I2.2," and spanning ",F5.2," hours.")'
                        ENDELSE
	            ENDIF ELSE BEGIN
	            PRINT,scn,sfn,date $
                         ,FORMAT='("The .ray file for camera ",A,",' $
				+'filter ",A,", on date ",A,' $
				+'", does not contain any images")'
	            PRINT,hh,mm,ss,span $
	                 ,FORMAT='("in the interval starting at",' $
				+'1X,I2.2,":",I2.2,":",I2.2," and ' $
				+'spanning ",F5.2," hours.")'
                    ENDELSE 
                ENDIF ELSE BEGIN
                PRINT,bcm,bft,date $
	             ,FORMAT='("The .bkg file for camera ",A,",' $
			    +'filter ",A,", on date ",A,", does ' $
			    +'not contain")'
                PRINT,span $
		     ,FORMAT='("any background images produced ' $
			    +'from ",I2," hour intervals")'
                ENDELSE 
            ENDIF ELSE BEGIN
            PRINT,bcm,bft,date $
                 ,FORMAT='("The .bkg file for camera ",A,",filter ",' $
			+'A,", on date ",A,", does not exist.")'
            ENDELSE
        ENDIF ELSE BEGIN
        PRINT,scn,sfn,date $
             ,FORMAT='("The .ray file for camera ",A,",filter ",A,' $
		    +'", on date ",A,", does contain any")'
        PRINT,exp $
	     ,FORMAT='(1X,I2,1X,"second exposure images.")'
        ENDELSE
    ENDIF ELSE BEGIN
    PRINT,scn,sfn,date $
         ,FORMAT='("The .ray file for camera ",A,",filter ",A,", on ' $
		+'date ",A,", does not exist.")'
    ENDELSE

END