; CNA_LINFIT: Compute a linear fit of the form y = a + bx to a set of
; data points x, y with individual standard deviations sig. Returned
; are a, b and their respective probable uncertainties siga and sigb.
; Taken from Numerical Recipes (with modifications to handle masked data).

pro cna_linfit, x, y, sig, mask, ndata, a, b, siga, sigb

  ; Local variables.
  ; int i
  ; int ngood
  ; int goodpt
  ; float ss, sx, sy, st2, t, wt, sxoss

  ; Initialize accumulators and results.
  ss = 0.0
  sx = 0.0
  sy = 0.0
  st2 = 0.0
  a = 0.0
  b = 0.0
  siga = 0.0
  sigb = 0.0

  ; Count number of non-masked points.
  ngood = 0
  for i = 0, ndata-1 do begin
    if (mask(i) eq 0) then begin
      goodpt = i
      ngood = ngood + 1
    endif
  endfor

  ; Check for trivial cases.
  if (ngood eq 0) then begin
    return
  endif else if (ngood eq 1 and x(goodpt) ne 0) then begin
    b  = y(goodpt) / x(goodpt)
    sigb = sig(goodpt) / x(goodpt)
    return
  endif

  ; Accumulate sums.
  for i = 0, ndata-1 do begin
    if (sig(i) ne 0 and mask(i) eq 0) then begin
      wt = 1.0 / (sig(i)*sig(i))
      ss = ss + wt
      sx = sx + x(i)*wt
      sy = sy + y(i)*wt
    endif
  endfor

  if (ss ne 0) then begin
    sxoss = sx / ss
  endif else begin
    sxoss = 0
  endelse

  for i = 0, ndata-1 do begin
    if (sig(i) ne 0 and mask(i) eq 0) then begin
      t = (x(i) - sxoss) / sig(i)
      st2 = st2 + t*t
      b = b + t*y(i)/sig(i)
    endif
  endfor

  ; Solve for b and sigb.
  if (st2 ne 0) then begin
    b = b / st2
    sigb = sqrt (1.0 / st2)
  endif else begin
    b = 0
    sigb = 0
  endelse

  ; Solve for a and siga.
  if (ss ne 0 and st2 ne 0) then begin
    a = (sy - sx*b) / ss
    siga = sqrt ((1.0+sx*sx/(ss*st2)) / ss)
  endif else begin
    a = 0
    siga = 0
  endelse

end


; CNA_FITSAMPS: Fit accumulating counts vs. time to compute mean countrate,
; iteratively rejecting CR hits and refitting until no new samples are
; rejected. The rejection test is based on the input error value for
; each sample.

pro cna_fitsamps, nsamp, sci, err, dq, time, thresh, out_sci, out_err, $
		  out_samp, out_time, i, j

@calnica_common

  ; Local variables.
  ; int k      ; loop index.
  ; int nrej    ; number of rejected samples.
  ; float *tsci    ; list of cleaned fluxes.
  ; float *terr    ; list of cleaned errors.
  ; float *ttime    ; list of cleaned times.
  ; float bad_sci    ; sci  value of bad sample(s).
  ; float bad_time    ; time value of bad sample(s).
  ; int   bad_samp    ; number of bad samples.
  ; float a, siga    ; linear fit zeropoint and error.
  ; float b, sigb    ; linear fit slope and error.
  ; float *diff    ; sample differences relative to fit.
  ; float prev_diff  ; temporary difference value.

  ; Function definitions.
  ; void cna_linfit (float , float , float , short , short, float ,
  ;      float , float , float )

  ; Initialize the output values.
  out_sci  = 0.0
  out_err  = 0.0
  out_samp = 0.0
  out_time = 0.0

  ; Allocate memory for local arrays.
  tsci   = fltarr(nsamp)
  terr   = fltarr(nsamp)
  ttime  = fltarr(nsamp)
  diff   = fltarr(nsamp)

  ; Fit and test samples for this pixel until no new samples are rejected.
  nrej = 1
  while (nrej gt 0) do begin

    ; Copy non-rejected sample values to tmp arrays.
    bad_samp = 0
    bad_time = 0
    bad_sci  = 0
    for k = 0, nsamp-1 do begin
      ; Temporarily flag zeroth read so it's not used in fit.
      ; but use the zeroth read data if it's OK and the first
      ; read is saturated (Vsn 3.2)
      if (k eq 0) then begin
	tsci(k) = sci(k)
	terr(k) = err(k)
	ttime(k) = time(k)
        if (not(dq[0] eq 0 and (dq[1] and SATURATED) gt 0)) then begin
	  dq(k) = dq(k) or BADPIX
	endif
      endif else if (dq(k) eq 0) then begin
	; Copy unflagged samples to tmp arrays if a previous
	; sample was flagged, subtract its values from
	; accumulated counts and time, but use original errs..
	if (bad_samp eq 0) then begin
	  tsci(k)  = sci(k)
	  ttime(k) = time(k)
	endif else begin
	  tsci(k)  = sci(k)  - bad_sci
	  ttime(k) = time(k) - bad_time
	endelse
	terr(k) = err(k)
      endif else if ((dq(k) and  BADPIX) eq 0) then begin
	; Accumulate bad sci and time values for flagged samples,
	; but only if they're not flagged as BADPIX. BADPIX is
	; used to flag electronic noise which does not affect
	; later samples.
	bad_sci  = bad_sci  + (sci(k)  - sci(k-1))
	bad_time = bad_time + (time(k) - time(k-1))
	bad_samp = bad_samp + 1
      endif


    endfor

    ; Compute mean countrate using linear fit.
    cna_linfit, ttime, tsci, terr, dq, nsamp, a, b, siga, sigb

    ; Reject outliers.
    nrej    = 0
    out_samp = 0
    out_time = 0

    ; First, compute difference of each unflagged sample
    ; relative to the fitted counts.
    for k = 1, nsamp-1 do begin
      if (terr(k) ne 0 and dq(k) eq 0) then begin
        diff(k) = (tsci(k)-(a+b*ttime(k))) / terr(k)
      endif else begin
        diff(k) = 0
      endelse
    endfor

    ; Next, do a special check for electronic noise spikes.
    ; These show up as large deviations (positive or negative)
    ; for a single sample relative to the samples on either
    ; side of it.
    for k = 2, nsamp-2 do begin
      if (dq(k-1) eq 0 and dq(k) eq 0 and dq(k+1) eq 0) then begin

        ; Check for positive spikes.
        if (diff(k)-diff(k-1) gt thresh  and  $
          diff(k)-diff(k+1) gt thresh) then begin
          if ((tsci(k)-tsci(k-1)) gt 3*terr(k-1) and $
            (tsci(k)-tsci(k+1)) gt 3*terr(k+1)) then begin
            dq(k) = dq(k) or BADPIX
            nrej = nrej + 1
	  endif
	endif

	; Check for negative spikes.
	if ((diff(k-1)-diff(k)) gt thresh and $
	  (diff(k+1)-diff(k)) gt thresh) then begin
	  if ((tsci(k-1)-tsci(k)) gt 3*terr(k-1) and $
	    (tsci(k+1)-tsci(k)) gt 3*terr(k+1)) then begin
	    dq(k) = dq(k) or BADPIX
	    nrej = nrej + 1
	  endif
	endif
      endif
    endfor

    ; Finally, check for CR hits by searching for large
    ; negative-to-positive going changes in diff.
    prev_diff = 0
    for k = 2, nsamp-1 do begin
      if (dq(k-1) eq 0) then prev_diff = diff(k-1)
      if (dq(k) eq 0) then begin
        if ((diff(k)-prev_diff) gt thresh) then begin 
          dq(k) = dq(k) or CR_HIT
          prev_diff = diff(k)
          nrej = nrej + 1
        endif
      endif
    endfor

  endwhile ; iterate the fit.

  ; Set output values.
  out_sci = b
  out_err = sigb

  ; Compute number of good samples and total integration time.
  for k = 0, nsamp-1 do begin
    if (dq(k) eq 0) then begin
      out_samp = out_samp + 1
      out_time = ttime(k)
    endif
  endfor

  ; Remove BADPIX flag from zeroth read.
  if (dq(0) and BADPIX) gt 0 then dq(0) = dq(0) - BADPIX

end


; CNA_SCRID: Identify cosmic ray hits in single NICMOS images.
; Detected hits are only flagged in the DQ array SCI image pixel
; values are NOT changed..

function cna_scrid, nic, input

; Arguments:
;  nic   i: NICMOS info structure
;  input   i: input image


@calnica_common

  ; Local variables.

  ; Function definitions.

  ; Routine not available yet.
  print, ' WARNING, CRIDCALC not yet implemented will be skipped'

  nic.CRID.corr = OMIT

  ; Successful return.
  return, 0
end


; CNA_MCRID: Identify cosmic ray hits in a stack of MultiACCUM samples.
; Detected hits are flagged in the DQ arrays of the individual samples,
; but the SCI images in the individual samples are NOT changed.
; A final single-group image is produced using only unflagged data
; from the individual samples. The accumlated counts vs. time are fitted
; with a first-order (linear) polynomial using standard linear regression.
; Outliers are identified by their distance from the fit and rejected..

function cna_mcrid, nic, input, crimage

; Arguments:
;  nic   i: NICMOS info structure
;  input  io: input image stack
;  crimage   o: CR reject output image (single-group)

@calnica_common

  THRESH = 4.0

  ; Local variables.
  ; short i, j, k        ; pixel and loop indexes.
  ; short nsamp          ; number of samples for pixel.
  ; float current_sci    ; current sci value.
  ; float current_err    ; current err value.
  ; short current_dq     ; current dq  value.
  ; float *sci           ; list of sci  values for pixel.
  ; float *err           ; list of err  values for pixel.
  ; short *dq            ; list of dq   values for pixel.
  ; float *time          ; list of time values for pixel.
  ; float out_sci        ; output sci  value.
  ; float out_err        ; output err  value.
  ; short out_dq         ; output dq   value.
  ; short out_samp       ; output samp value.
  ; float out_time       ; output time value.

  ; Function definitions.
  ; cna_copyGroup (SingleNicmosGroup , SingleNicmosGroup )
  ; void cna_fitsamps (short, float , float , short , float , float,
  ; float , float , short , float , short, short)

  ; Allocate memory for local arrays.
  sci  = fltarr(nic.ngroups)
  err  = fltarr(nic.ngroups)
  dq   = intarr(nic.ngroups)
  time = fltarr(nic.ngroups)

  ; Initialize the output crimage by copying the first input
  ; group to the crimage. 
  ; 06Mar2000, I. Barg - 'cna_copyGroup' is done by CNA_CRIDCALC.PRO.
  ; cna_copyGroup, input(0), crimage

  ; Set CR reject threshold use default if no user input.
  if (nic.crthresh eq 0.0) then nic.crthresh = THRESH

  ; Print info to processing log.
  print, ' CRIDCALC using ',nic.crthresh,' sigma rejection threshold'

  PRINT,'Start mcrid routine:'
  print,"crimage.sci_dat.sci_array.data[0:5,0:5]:"
  print,FORMAT = '(6(F12.6))',crimage.sci_dat.sci_array.data[0:5,0:5]

  ; Loop over image array, computing mean countrate at each pixel.
  for j = 0, input(0).sci_dat.sci_array.tot_ny-1 do begin
    for i = 0, input(0).sci_dat.sci_array.tot_nx-1 do begin

      ; Initialize the output image DQ value.
      out_dq = 0

      ; Create list of samples for this pixel.
      nsamp = 0
      for k = nic.ngroups-1, 0, -1 do begin

	; Retrieve SCI, ERR, and DQ values for this sample.
	current_sci = input(k).sci_dat.sci_array.data(i,j)
	current_err = input(k).err_dat.err_array.data(i,j)
	current_dq  = input(k).dq_dat.dq_array.data(i,j)

	; Temporarily convert countrates back to counts.
	if (nic.bunit(0) eq COUNTRATE) then begin

	  ; Use sampzero for exposure time of zeroth read
	  ; if ZSIGCORR is on, otherwise use regularly
	  ; assigned value of exptime (Version 3.2)

	  if k eq nic.ngroups-1 and nic.ZSIG.corr eq PERFORM then begin
	    current_sci = current_sci * nic.sampzero
	    current_err = current_err * nic.sampzero
	  endif else begin
	    current_sci = current_sci * nic.exptime(k)
	    current_err = current_err * nic.exptime(k)
	  endelse
	endif

	; Propagate DQ values to output DQ.
	out_dq = out_dq or current_dq

        ; Remove ZEROSIG and GROT bits from DQ, if present,
        ; because those samples are OK to use here (Vsn 3.2)
        if ((current_dq and ZEROSIG) gt 0) then begin
           current_dq = current_dq - ZEROSIG
	endif
        if ((current_dq and GROT) gt 0) then $
           current_dq = current_dq - GROT

	sci(nsamp)  = current_sci
	err(nsamp)  = current_err
	dq(nsamp)   = current_dq
	time(nsamp) = nic.exptime(k)
        
        ; Use sampzero for zeroth read exptime (Vsn 3.2)
        if (k eq nic.ngroups-1 and nic.ZSIG.corr eq PERFORM) then $
           time(nsamp) = nic.sampzero

        nsamp = nsamp + 1
      endfor

      ; Do iterative rejection and computation of slope.
      cna_fitsamps, nsamp, sci, err, dq, time, nic.crthresh, $
          out_sci, out_err, out_samp, out_time, i, j

      ; If at least one good sample was found, reset output DQ
      ; value to zero. But if the original data contained a
      ; ZEROSIG or GROT value, leave those values set in output DQ
      ; (Version 3.2) 
      if (out_samp gt 0) then begin
        current_dq = 0
	if (out_dq and ZEROSIG) gt 0 then begin
	  current_dq = current_dq + ZEROSIG
        endif
	if (out_dq and GROT) gt 0 then current_dq = current_dq + GROT
	out_dq = current_dq
      endif

      ; Convert results back to counts if necessary.
      if (nic.bunit(0) eq COUNTS) then begin
        out_sci = out_sci * out_time
        out_err = out_err * out_time
      endif

      ; Store final values in output crimage.
      crimage.sci_dat.sci_array.data(i,j) = out_sci
      crimage.err_dat.err_array.data(i,j) = out_err
      crimage.dq_dat.dq_array.data(i,j) = out_dq
      crimage.smpl_dat.smpl_array.data(i,j) = out_samp
      crimage.intg_dat.intg_array.data(i,j) = out_time

      ; Set CR_HIT DQ for all samples following a hit.
      for k = 0, nic.ngroups-2 do begin
        if ((dq(k) and CR_HIT) gt 0) then dq(k+1) = dq(k+1) or CR_HIT
      endfor

      ; Update input DQ values for detected outliers.
      for k = nic.ngroups-1, 0, -1 do begin
        if ((dq(nic.ngroups-1-k) and CR_HIT) gt 0) then begin
	  input(k).dq_dat.dq_array.data(i,j) =  $
	    input(k).dq_dat.dq_array.data(i,j) or CR_HIT
	endif

        if ((dq(nic.ngroups-1-k) and BADPIX) gt 0) then begin
	  input(k).dq_dat.dq_array.data(i,j) =  $
	    input(k).dq_dat.dq_array.data(i,j) or BADPIX
	endif
      endfor

    endfor ; end of loop over nx.
  endfor ; end of loop over ny.

  PRINT,'After mcrid routine:'
  print,"crimage.sci_dat.sci_array.data[0:5,0:5]:"
  print,FORMAT = '(6(F12.6))',crimage.sci_dat.sci_array.data[0:5,0:5]

  ; Successful return.
  return, 0
end


; IDL_CRIDCALC: Identify cosmic ray hits in NICMOS images. This routine
; calls separate subroutines for working with single (ACCUM,RAMP,BRTOBJ)
; images or multiple (MULTIACCUM) images.
;
; Revision history:
; H.Bushouse  Feb. 1996  Written for Build 2 (Version 2.0)
; H.Bushouse  05-Jun-1997  New cna_mcrid routine that performs a linear
;        fit to flux vs time and uses NOISFILE to
;        compute errors. (Version 2.3)
; H.Bushouse  23-Jul-1997  Changed cna_linfit to return zero when given 0
;        data points and return the input value when
;        given 1 data point. (Version 3.0)
; H.Bushouse  28-Jul-1997  Changed NicInfo.bunit from scalar to vector
;        (Version 3.0)
; H.Bushouse  11-Sep-1997  Added use of nic.crthresh to allow user input
;        (Version 3.0)
; H.Bushouse  29-Sep-1997  Complete rewrite of cna_mcrid - eliminated use of
;        asigclip,aerrclip,merrclip,mean,wtmean routines
;        to compute mean differences and reject samples.
;        Rejection and fit to samples are now both done
;        together using new cna_fitsamps routine. No need
;        for NOISFILE anymore. (Version 3.0)
; D.Lytle 12-Nov-1997 IDL translation
; H.Bushouse   15-Feb-1998     Modified n_mcrid and n_fitsamps for special
;         handling of MultiAccum zeroth reads when first read is saturated,
;         including special handling of new ZEROSIG and GROT DQ values
;         (Version 3.2)
; H.Bushouse   05-May-1998     Modified n_linfit to fix division by zero bug
;          when ngood=1 and x[goodpt]=0 (Version 3.2)
; D.Lytle 30-May-1998 IDL update
; I. Barg 06Mar2000 - CNA_CRIDCALC.PRO renamed to IDL_CRIDCALC.PRO.  
;          CNA_CRIDCALC.PRO now either calls an external "C" cridcalc routine
;          or this IDL_CRIDCALC procedure.


function idl_cridcalc, nic, input, crimage

; Arguments:
;  nic     i: NICMOS info structure
;  input  io: input image(s)
;  crimage o: CR reject output image (for MultiACCUM only)

@calnica_common

  ; Function definitions.
  ; int cna_mcrid (NicInfo , MultiNicmosGroup , SingleNicmosGroup )
  ; int cna_scrid (NicInfo , MultiNicmosGroup )

  ; Call multi-image routine for MultiACCUM datasets.
  if (nic.obsmode eq MULTIACCUM) then begin
    if (cna_mcrid (nic, input, crimage) ne 0) then return, 8086
  endif else begin
    ; Call single-image routine for all others.
    if (cna_scrid (nic, input) ne 0) then return, 8086
  endelse

  ; Successful return.
  return, 0
end