; CNA_NLINCORR -- Correct NICMOS science image values for detector nonlinearity.
; The SCI and ERR arrays are updated, and the DQ values are propogated.
; The SAMP and TIME arrays are not modified.

; Revision history:
; 10-Jan-2000, I. Barg - Modified structure names (and tag names) to match 
;             "/iraf/stsdas/lib/hstio.h".
; 08-Jan-2001, I. Barg - updated to reflect new NICMOS linearity FILES
;             released by STScI on 01-01-2001.
; 10-Jan-2001, I. Barg - Changed all array indexes from () to [].

function cna_nlincorr, nic, nlin, input, zsig

@calnica_common

  ; Arguments:
  ;	nic	 i: NICMOS info structure
  ;	nlin	 i: nonlinearity reference data
  ;	input	io: input image to be corrected
  ;	zsig	io: MULTIACCUM zero-read signal image

  ; Add in the MULTIACCUM zero-read signal only if
  ; ZSIG step is turned on, and only for groups other than the
  ; zeroth-read itself (Vsn 3.2)
  if (nic.ZSIG.corr eq PERFORM and nic.group ne nic.ngroups) then begin
    input.sci.data.data = input.sci.data.data + $
      zsig.sci.data.data
    mask = where((zsig.dq.data.data and ZEROSIG) gt 0,nzsig)
    if nzsig gt 0 then begin
      tempdq = input.dq.data.data
      tempdq[mask] = tempdq[mask] or ZEROSIG
      input.dq.data.data = tempdq
    endif
  endif

  ; Find where the data pixels fall into the various regimes.
  test1 = input.sci.data.data lt nlin.nodes[0].data.data
  test2 = (input.sci.data.data ge nlin.nodes[0].data.data) and $
	  (input.sci.data.data lt nlin.nodes[1].data.data)
  test3 = input.sci.data.data ge nlin.nodes[1].data.data

  ; See how many saturated pixels, report it.
  tmp = where(test3 eq 1,nsp)
  print, ' NLINCORR detected ',nsp,' saturated pixels'

  ; test1, test2, test3, and test4 will be arrays containing ones and zeros.
  ; These can be used as masks.

  r_1_data = input.sci.data.data * float(test1)
  r_2_data = input.sci.data.data * float(test2)
  r_3_data = input.sci.data.data * float(test3)

  r_1_err = input.err.data.data * float(test1)
  r_2_err = input.err.data.data * float(test2)
  r_3_err = input.err.data.data * float(test3)
  r_2_err_save = r_2_err

  r_1_dq = (nlin.dqual.data.data) * test1
  r_2_dq = (nlin.dqual.data.data) * test2 
  r_3_dq = (nlin.dqual.data.data or SATURATED) * test3

  ; Update data quality.
  input.dq.data.data = input.dq.data.data or $
			 (r_1_dq + r_2_dq + r_3_dq)

  ; Update the errors.
  r_2_data2 = r_2_data * r_2_data
  
  e11 = nlin.error[0].data.data
  e22 = nlin.error[1].data.data
  e12 = nlin.error[2].data.data
  e33 = nlin.error[3].data.data
  e13 = nlin.error[4].data.data
  e23 = nlin.error[5].data.data
  
  r_2_err = e11 + r_2_data2*e22 + r_2_data2*r_2_data2*e33 + $
            2*(abs(r_2_data)*e12 + r_2_data2*e13 + abs(r_2_data)*r_2_data2*e23)
  
  r_2_err = r_2_err * float(test2)

  ; Now we will apply the quadratic. Assume the coefficients are in arrays.
  ; The saturated pixels and pixels in regime '1' just pass through.
  r_2_data = (nlin.coeff[0].data.data + nlin.coeff[1].data.data*r_2_data + $
              nlin.coeff[2].data.data*r_2_data2) * r_2_data
  
  ; Next, recombine the data and errors.
  input.sci.data.data = r_1_data + r_2_data + r_3_data

  ; Remove the MULTIACCUM zero-read signal only if
  ; ZSIG step is turned on, and only for groups other than the
  ; zeroth-read itself (Vsn 3.2)
  if (nic.ZSIG.corr eq PERFORM and nic.group ne nic.ngroups) then begin
    input.sci.data.data = input.sci.data.data - $
      zsig.sci.data.data
  endif

  ;ee =  r_1_err  + r_2_err  + r_3_err
  ;eee = input.err.data.data * input.err.data.data + ee
  r_2_err = r_2_err_save * r_2_err_save + r_2_err
  pos = where (r_2_err ge 0.0,poscount)
  neg = where (r_2_err lt 0.0,negcount)
  if (poscount gt 0) then begin
    r_2_err[pos] = sqrt(r_2_err[pos])
  endif
  if (negcount gt 0) then begin
    r_2_err[neg] = 0.0
  endif

  input.err.data.data = r_1_err  + r_2_err  + r_3_err

  temp = where((input.dq.data.data and SATURATED) gt 0, numsat)
  print,' NLINCORR detected ',numsat,' saturated pixels;'

  return, 0 ; successful return

end