forward_function pds_sys_date

;==========================================================================
function pds_sys_date

; Function that calls the system to get the current UTC date and returns it
; in the standard PDS format  (YYYY-MM-DDTHH:MM:SS)

h1= systime(/utc)
cmonth = strmid(h1,4,3)

case cmonth of 
  'Jan': mon = 1
  'Feb': mon = 2
  'Mar': mon = 3
  'Apr': mon = 4
  'May': mon = 5
  'Jun': mon = 6
  'Jul': mon = 7
  'Aug': mon = 8
  'Sep': mon = 9
  'Oct': mon = 10
  'Nov': mon = 11
  'Dec': mon = 12
endcase
dates = strmid(h1,20,4)+"-"+string(mon,format='(i2.2)')+"-"+strmid(h1,8,2)+$
     "T"+strmid(h1,11,8)

return, dates
end
;==========================================================================

function dark_frame_generate, readout_mode, fpa_temp, sim_temp,$
         cover_temp, prim_mirror_temp, sec_mirror_temp,int_seq_gap_time, $
         integ_time, prev_readout_mode, readout_slot, orig_image,$
         dark_header_info, MDARK_FN=mdark_fn
;;+
;
; NAME:
;	dark_frame_generate
;
; PURPOSE:
;	Generate a dark frame for the HRI IR for any set of
;         temperatures and readout slots
;
; EXPLANATION:
;	Produces a synthetic dark frame that is matched to a given set of
;       conditions for an HRI IR image of interest.  The dark frame can then
;       be subtracted from the image.  
;
;       The dark image is produced using a master dark frame, which is
;       scaled to the specific temperatures of the Focal Plane array
;       and the SIM bench (average of the two prism temperatures) and 
;       the integration time.
;
;       If the image is number 1-7 in the sequence, then the dark
;       frame is scaled to account for the temporal dependence of the
;       dark frames.  (From the 8th frame on, it is assumed that the
;       dark level is no longer changing.)  For the temporal scaling,
;       it is assumed that the inter-sequence gap is long enough that
;       the dark has reached the plateau for that mode.  As a check,
;       the value for the real inter-sequence gap is passed to the
;       routine, and it is used to see if the plateau is likely
;       to have been reached.  If not, then a flag is set to indicate
;       that the dark level that was removed is probably incorrect.  
;
;       Compare the dark frame to the original image.  If the dark frame
;       removes too much signal, then adjust the dark level so that it
;       only removes 99% of the signal.  This keeps the routine from removing
;       too much of the signal and having large negative fluxes in the end.
;       NOTE: Nothing is being done to the original image in this routine.
;       It is simply being used for comparison so that the routine does not
;       remove too much of the dark frame.  Actual dark removal is done 
;       outside, in the general pipeline sequence.
;
;       The reference areas at the bottom and sides of each frame are
;       set to zero in the dark frame to preserve the reference areas in
;       the image.
;
;       The master dark files are stored in a separate directory.  
;       The name of that directory must be given in the environment
;       variable HRIIR_MASTER_DARK_DIR which is read by this routine
;       so that the dark files can be found.
;
;
; CALLING SEQUENCE:
;	result = dark_frame_generate (readout_mode, fpa_temp, sim_temp, $
;         cover_temp, prim_mirror_temp, sec_mirror_temp, $
;         int_seq_gap_time, integ_time, prev_readout_mode, readout_slot, $
;         orig_image, dark_header_info, 
;
; INPUTS:
;	readout_mode = (Integer) The mode that was used for reading out the
;	               image. 
;
;       fpa_temp = (Real) Temperature of the focal plane array (K), corrected
;                  for the manganin wire offset. (keyword IRFPAT)
;
;       sim_temp = (Real) Temperature of the SIM Bench (K). (keyword OPTBENT)
;
;       cover_temp = (Real) Temperature of the SIM cover (K) (keyword COVERT - not used)
;
;       prim_mirror_temp = (Real) Primary Mirror Temperature (K) (keyword PRIMIRT - not used)
;
;       sec_mirror_temp = (Real) Primary Mirror Temperature (K) (keyword PRIMIRT - not used)
;
;       int_seq_gap_time = (Real) Inter-sequence gap time.  The time (msec) 
;                          between the last frame in the previous sequence 
;                          and the first frame of the current sequence.  (- not used)
;
;       integ_time = (Real) Exposure integration time (msec) (keyword INTTIME)
;
;       prev_readout_mode = (Integer) The readout mode of the previous sequence.  (- not used)
;
;       readout_slot = (Integer) The number of the image in the current
;                       sequence.  The first frame is in slot 1 
;
;       orig_image = (real) The image for which the dark is being generated.  
;                           Used here only for comparison to make sure that
;                           the routine doesn't remove too much dark.
;
; OPTIONAL INPUTS:
;	MDARK_FN - Filename of the master dark file to use. If this is not set,
;		the dark will be taken from the directory specified by the
;		HRIIR_MASTER_DARK_DIR environemnt variable.
;
; OUTPUTS:
;       result = (Real) Array containing the synthetic dark frame, produced
;                for the given conditions. 
;
;       dark_header_info = (String array) Information to be added to the 
;                          image header describing the dark frame that was 
;                          created and removed from the image.  **** Format
;                          still needs to be defined.
;       
; RESTRICTIONS:
;    I have included no checks for number or format of input parameters, 
;    etc., on the assumption that the pipeline process will be consistent 
;    in what is input into the routine.  
;
;    Modes 5 and 6 will currently return a blank image
;    Mode 7 is not currently supported
;
; ENVIRONMENT VARIABLES:
;     HRIIR_MASTER_DARK_DIR  - directory where the master dark frames
;                              are stored.
;
; MODIFICATION HISTORY:
;	WRITTEN, Tony Farnham  30 November 2004   VERSION  1.00
;	
;	30 Nov 2004 Mark Desnoyer - Added ability to specify the filename of
;		the master dark file to use
;
;       01 June 2005 Tony Farnham - Added keywords to the calling line for 
;                placeholders.  Nothing done internally yet.
;                Added: cover_temp, prim_mirror_temp, sec_mirror_temp,
;                       integ_time
;                Removed: elec_time
;                Changed: int_seq_gap_time (changed from sec to msec)
;
;       13 June 2005 Tony Farnham - Updated coefficients  
;      
;       14 June 2005 Tony Farnham - Added orig_image to the input parameters
;             and added the code to check and make sure that too much dark is
;             not removed.  Removed adjustment for image the image readout
;             slot (all images are now scaled by 1.)  Also changed
;             temperature scaling so only the temp terms are scaled (now take out
;             the constant c0 for the scaling).
;
;       16 June 2005 Tony Farnham - Updated the coefficients for final
;                                   pipeline version
;                                   
;-
On_error, 2

;---- set the environment variable HRIIR_MASTER_DARK_DIR to contain the 
;---- location where the master dark files are stored
master_dir = getenv('HRIIR_MASTER_DARK_DIR')
master_dir = strcompress(master_dir,/remove_all)

dates = pds_sys_date()
rmode = readout_mode
Tfpa = fpa_temp
Tsim = sim_temp
timeisg = int_seq_gap_time/1000.
prmode = prev_readout_mode
rslot = readout_slot
inttime = integ_time

;----- set up the coefficients for each of the different modes:
;  a0, a1, b0, b1, c0 are the coefficients for measuring the
;     temperature dependence.

;  Coefficients computed for all modes combined, in counts per msec
c0 =  0.0583
a0 =  3.006e+10  &   a1 = -3384.40
b0 =  6.563e+07  &   b1 = -2002.00

CASE rmode of
  1: begin master_file = 'master_dark_mode1.fits' &$
     r4x1 = 96   & r4x2 = 160    & r4y1 = 96  & r4y2 = 160  &$
     f1x1 = 0    & f1x2 = 1      & f1y1 = 0   & f1y2 = 255  &$
     f2x1 = 510  & f2x2 = 511    & f2y1 = 0   & f2y2 = 255  &$
     f3x1 = 0    & f3x2 = 511    & f3y1 = 0   & f3y2 = 2 & end  
  2: begin  master_file = 'master_dark_mode2.fits' &$
     r4x1 = 96   & r4x2 = 160    & r4y1 = 32  & r4y2 = 96   &$
     f1x1 = 0    & f1x2 = 1      & f1y1 = 0   & f1y2 = 127  &$
     f2x1 = 510  & f2x2 = 511    & f2y1 = 0   & f2y2 = 127  &$
     f3x1 = 0    & f3x2 = 511    & f3y1 = 0   & f3y2 = 2  & end    
  3: begin  master_file = 'master_dark_mode3.fits' &$
     r4x1 = 96   & r4x2 = 160    & r4y1 = 3   & r4y2 = 61   &$
     f1x1 = 0    & f1x2 = 1      & f1y1 = 0   & f1y2 = 63   &$
     f2x1 = 510  & f2x2 = 511    & f2y1 = 0   & f2y2 = 63   &$
     f3x1 = 0    & f3x2 = 511    & f3y1 = 0   & f3y2 = 2  & end    
  4: begin  master_file = 'master_dark_mode4.fits' &$
     r4x1 = 191  & r4x2 = 319    & r4y1 = 192 & r4y2 = 320  &$
     f1x1 = 0    & f1x2 = 4      & f1y1 = 0   & f1y2 = 507  &$
     f2x1 = 1019 & f2x2 = 1023   & f2y1 = 0   & f2y2 = 507  &$
     f3x1 = 0    & f3x2 = 1023   & f3y1 = 0   & f3y2 = 5  & end    
  5: begin 
     ; -- send back a blank image frame
     dark_frame = fltarr(512,256)*0.0
     dark_header_info = 'COMMENT  MODE 5 - Dark frame is not removed'
     return,dark_frame
     end
  6: begin  
     ; -- send back a blank image frame
     dark_frame = fltarr(1024,512)*0.0
     dark_header_info = 'COMMENT  MODE 6 - Dark frame is not removed'
     return,dark_frame
     end
else: return,0
ENDCASE

; ---- set up the name of the master dark file
if keyword_set(MDARK_FN) then $
     fits_file = MDARK_FN $
else $
     fits_file = master_dir+master_file

;----- read in the master dark frame
fits_read,fits_file,mdark,mdarkh
darkver = fxpar(mdarkh,'VERSION')

;----- scale the master dark frame to the relevant temperatures
rel_dark = a0*exp(a1/Tsim)  + b0*exp(b1/Tfpa)
dark_frame = ((mdark-c0)*rel_dark + c0) * inttime

;  NOTE removed the scaling for the temporal dependence
;;----- scale for the temporal dependence of the first 7 frames
;if rslot ge 1 and rslot le 7 then timescal = sc(rslot)
timescal = 1.000                     ; don't comment out - used in header
;dark_frame = dark_frame * timescal 

;----- Check the dark frame to see if too much will be removed - use region 4
dkavg=moment(dark_frame[r4x1:r4x2,r4y1:r4y2],/NAN)
imavg=moment(orig_image[r4x1:r4x2,r4y1:r4y2],/NAN)
 
; if too much dark is removed, scale to 99% of the image value
if dkavg(0) gt imavg(0) then dark_frame = 0.99*imavg(0)*dark_frame/dkavg(0)

;----- set the reference areas to zero
dark_frame(f1x1:f1x2,f1y1:f1y2) = 0.
dark_frame(f2x1:f2x2,f2y1:f2y2) = 0.
dark_frame(f3x1:f3x2,f3y1:f3y2) = 0.

;=======================================================================
;  GENERATE HEADER INFORMATION
;=======================================================================
header = strarr(500)
blankline='                                                                                '
separator = 'COMMENT ***********************************************************************'
h1=blankline
hcnt = 0
header(hcnt) = blankline
hcnt = hcnt+1

header(hcnt) = separator
hcnt = hcnt+1
header(hcnt) = 'COMMENT **                    DARK LEVEL REMOVAL                             **'
hcnt = hcnt+1
header(hcnt) = separator
hcnt = hcnt+1
hcnt = hcnt+1

;;------ Describe how the dark frame was produced
header(hcnt) = 'COMMENT **                                                                   **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  The dark level was removed using a master dark frame for the     **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  given readout mode, and scaling it for the FPA and prism         **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  temperatures that were measured at the time of the observation.  **' 
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  This is then scaled by the integration time to create the dark.  **'
hcnt = hcnt+1
;header(hcnt) = 'COMMENT **  If the frame is one of the first seven in its particular         **'
;hcnt = hcnt+1
;header(hcnt) = 'COMMENT **  sequence, then a temporal correction is also applied to account  **'
;hcnt = hcnt+1
header(hcnt) = 'COMMENT **  for the overshoot in the first frames of the sequence.           **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  The resulting dark frame was then subtracted from the image to   **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  remove the background dark level.                                **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **                                                                   **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  The temperature scaling is done using the relation:              **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **    TEMPSCAL = C_0 + A_0 * exp(A_1/Tprism) + B_0 * exp(B_1/Tfpa)   **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **  where the coefficients and temperatures are listed below         **'
hcnt = hcnt+1
header(hcnt) = 'COMMENT **                                                                   **'
hcnt = hcnt+1
;header(hcnt) = 'COMMENT **  The time scaling is done from an empirical fit to the data.      **'
;hcnt = hcnt+1
header(hcnt) = 'COMMENT **                                                                   **'
hcnt = hcnt+1
header(hcnt) = separator
hcnt = hcnt+1
;; ----- Input the frame name, version number and readout mode
header(hcnt) = blankline
hcnt = hcnt+1
header(hcnt) = "DARKDATE= '"+dates+"' / Creation date of the dark frame"
hcnt = hcnt+1
header(hcnt) = "MDARKNAM= '"+master_file+"' / Name of the master dark frame used"
hcnt = hcnt+1
header(hcnt) = 'MDARKVER= '+string(darkver,format='(15x,f5.2)')+' / Version number of the master dark frame used'
hcnt = hcnt+1
header(hcnt) = blankline
hcnt = hcnt+1
header(hcnt) = 'DARKMODE= '+string(rmode,format='(19x,i1)')+' / Readout mode used to create the dark frame'
hcnt = hcnt+1
header(hcnt) = "PREVMODE= "+string(prmode,format='(19x,i1)')+$
    " / Previous readout mode"
hcnt = hcnt+1
header(hcnt) = "READSLOT= "+string(rslot,format='(15x,i5)')+$
    " / Readout slot"
hcnt = hcnt+1
header(hcnt) = "TEMPFPA = "+string(Tfpa,format='(12x,f8.3)')+$
    " / Focal Plane Array Temperature"
hcnt = hcnt+1
header(hcnt) = "TEMPSIM = "+string(Tsim,format='(12x,f8.3)')+$
    " / SIM Bench Temperature"
hcnt = hcnt+1
header(hcnt) = "TIMEISG = "+string(timeisg,format='(10x,f10.3)')+$
    " / Inter-sequence Gap time (sec)"
hcnt = hcnt+1
header(hcnt) = blankline
hcnt = hcnt+1

header(hcnt) = 'COMMENT    Dark level temperature coefficients used:'
hcnt = hcnt+1
header(hcnt) = "COEFF_A0= "+string(a0,format='(8x,e12.4)')+$
    " / Dark level temperature coefficient A_0"
hcnt = hcnt+1
header(hcnt) = "COEFF_A1= "+string(a1,format='(12x,f8.2)')+$
    " / Dark level temperature coefficient A_1"
hcnt = hcnt+1
header(hcnt) = "COEFF_B0= "+string(b0,format='(8x,e12.4)')+$
    " / Dark level temperature coefficient B_0"
hcnt = hcnt+1
header(hcnt) = "COEFF_B1= "+string(b1,format='(12x,f8.2)')+$
    " / Dark level temperature coefficient B_1"
hcnt = hcnt+1
header(hcnt) = "COEFF_C0= "+string(c0,format='(11x,f9.4)')+$
    " / Dark level temperature coefficient C_0"
hcnt = hcnt+1
header(hcnt) = blankline
hcnt = hcnt+1
header(hcnt) = 'COMMENT    Scaling factors applied'
hcnt = hcnt+1
header(hcnt) = "TEMPSCAL= "+string(rel_dark,format='(9x,f11.4)')+$
    " / Scaling factor for temperature dependencies"
hcnt = hcnt+1
header(hcnt) = "TIMESCAL= "+string(timescal,format='(9x,f11.4)')+$
    " / Scaling factor for time dependency"
hcnt = hcnt+1

; ----- If the image is first or second in the sequence, put in a flag 
if rslot le 1 then begin
header(hcnt) = blankline
hcnt = hcnt+1
  header(hcnt) = 'COMMENT    ------------------------ CAUTION -----------------------------'
  hcnt = hcnt+1
  header(hcnt) = 'COMMENT    |          First or second frame in the sequence             |'
  hcnt = hcnt+1
  header(hcnt) = 'COMMENT    |           Dark level removed may be incorrect              |'
  hcnt = hcnt+1
  header(hcnt) = 'COMMENT    --------------------------------------------------------------'
  hcnt = hcnt+1
  header(hcnt) = "FLAG-DRK= 'First or second frame in a sequence - Dark level may not be correct'"
  hcnt = hcnt+1
endif 

header(hcnt) = blankline
hcnt = hcnt+1
header(hcnt) = separator
hcnt = hcnt+1
header(hcnt) = 'COMMENT **                  END OF DARK REMOVAL SECTION                      **'
hcnt = hcnt+1
header(hcnt) = separator
hcnt = hcnt+1
hcnt = hcnt+1

dark_header_info = header(0:hcnt)

return, dark_frame
end