; CNA_ZSIGCALC: Calculate the amount of signal from sources in a MultiAccum
; zeroth-read image. This is done by computing the amount of signal that
; arrived between the zeroth and first reads, and then scaling that to the
; exposure time of the zeroth read in order to estimate how much signal
; probably came in between the reset and the zeroth read. Note that this
; does not work well if the source flux is so high that it's beginning to
; go non-linear already in either the zeroth or first read. The raw data
; values in the zeroth and first readouts are also differenced with the
; super zero read reference image and checked for saturation in those
; readouts.
;
; Revision history:
; H.Bushouse  25-Jul-1997  Created (Version 3.0)
; H.Bushouse  18-Aug-1997  Modified to use super zero read image in
;        NlinData structure (Version 3.0)
; H.Bushouse  11-Sep-1997  Added use of nic.zsthresh to allow user input
;        of threshold value (Version 3.0)
; H.Bushouse  28-Oct-1997  Added checks for dummy ref data before
;        doing computations (Version 3.0)
; D.Lytle  10-Nov-1997 Corrupted into IDL!!
; I. Barg 10-Jan-2000 Modified structure names (and tag names) to match 
;             "/iraf/stsdas/lib/hstio.h".


function cna_zsigcalc, nic, input, mask, nois, dark, nlin, zsig

@calnica_common

ZTHRESH =  5.0  ; default signal clipping threshold (sigma).
ZTIME =  0.203  ; zeroth read exposure time (seconds).

; Arguments:
;  nic   i: NICMOS info structure
;  input i: input images
;  mask  i: bad pixel mask image
;  nois  i: read noise image
;  dark  i: dark current image
;  nlin  i: non-linearity correction data
;  zsig  o: zeroth read signal image

  ; Local variables.
  ; int i, j            ; pixel indexes.
  ; int nsat0, nsat1    ; number of saturated pixels.
  ; float zscale        ; exposure time scale factor.

  ; Function definitions.
  ; int cna_copyGroup (SingleNicmosGroup , SingleNicmosGroup )
  ; int cna_asub (SingleNicmosGroup , SingleNicmosGroup )
  ; int cna_aor  (SingleNicmosGroup , SingleNicmosGroup )
  ; int cna_noiscalc (NicInfo , SingleNicmosGroup , SingleNicmosGroup )

  ; Initialize counters.
  nsat0 = 0
  nsat1 = 0

  ; Set proper threshold value use default if no user input.
  if (nic.zsthresh eq 0.0) then begin
    nic.zsthresh = ZTHRESH
  endif

  ; Copy first read to zsig image.
  tempgroup = input.group[nic.ngroups-2]
  cna_copyGroup, tempgroup, zsig

  ; If any of the ref data is dummy, zero out the zsig image
  ; and abort the routine.
  if (nic.MASK.ref.dummy eq 1 or nic.NOIS.ref.dummy eq 1  or $
    nic.DARK.ref.dummy eq 1 or nic.NLIN.ref.dummy eq 1) then begin

    if (cna_asub (zsig, zsig) ne 0) then return, 8482

    print, 'Warning!, ZSIGCALC skipped because ref data is dummy'

    return, 0
  endif

  ; Subtract the zeroth read from the first read.
  if (cna_asub (zsig, (input.group[nic.ngroups-1])) ne 0) then return, 8482

  ; Mask bad pixels in the zsig image.
  if (cna_aor (zsig, mask) ne 0) then return, 8482

  ; Compute noise in the zsig image.
  if (cna_noiscalc (nic, nois, zsig) ne 0) then return, 8482

  ; Subtract dark current from the zsig image.
  if (cna_asub (zsig, dark) ne 0) then return, 8482

  ; Compute the scale factor to match the zeroth-read exptime.
  zscale = ZTIME / nic.exptime(nic.ngroups-2)

  ; Loop over the zsig image.
  for j = 0, zsig.sci.data.tot_ny-1 do begin
  for i = 0, zsig.sci.data.tot_nx-1 do begin

    ; Only look at unflagged pixels.
    if (zsig.dq.data.data[i,j] eq 0) then begin

      ; Scale pixels with greater than zsthresh*err signal to
      ; the exposure time of the zeroth read.
      if (zsig.sci.data.data[i,j] ge nic.zsthresh * $
	zsig.err.data.data[i,j]) then begin
	zsig.sci.data.data[i,j] = zsig.sci.data.data[i,j] * $
	zscale

      ; Mask out low signal pixels by setting them to zero.
      endif else zsig.sci.data.data[i,j] = 0.0

      ; Flag pixels that are saturated in the zeroth read.  
      if ((input.group[nic.ngroups-1].sci.data.data[i,j] - $
	nlin.zsci.data.data[i,j]) gt $
	nlin.nodes(1).data.data[i,j]) then begin
	tempgroup = input.group[nic.ngroups-1]
	zsig.dq.data.data[i,j] = zsig.dq.data.data[i,j] $
	or SATURATED
	tempgroup.dq.data.data[i,j] = $
	  tempgroup.dq.data.data[i,j] or SATURATED
	nsat0 = nsat0 + 1
	input.group[nic.ngroups-1] = tempgroup
      endif

      ; Flag pixels that are saturated in the first read.
      if (input.group[nic.ngroups-2].sci.data.data[i,j] - $
	nlin.zsci.data.data[i,j] gt $
	nlin.nodes(1).data.data[i,j]) then begin
	tempgroup = input.group[nic.ngroups-2]
	zsig.dq.data.data[i,j] = zsig.dq.data.data[i,j] $
	or SATURATED
	tempgroup.dq.data.data[i,j] = $
	  tempgroup.dq.data.data[i,j] or SATURATED
	nsat1 = nsat1 + 1
	input.group[nic.ngroups-2] = tempgroup
      endif

      ; Set flagged pixels to zero in the zsig image.
    endif else zsig.sci.data.data[i,j] = 0.0
  endfor
  endfor

  ; Report the number of saturated pixels detected.
  print, ' ZSIGCALC detected ',nsat0,' saturated pixels in 0th read'
  print, ' ZSIGCALC detected ',nsat1,' saturated pixels in 1st read'

  ; Successful return.
  return, 0

end