; Revision history:
; I. Barg 10-Jan-2000 Modified structure names (and tag names) to match 
;             "/iraf/stsdas/lib/hstio.h".

; CNA_NFLAGS -- Count the number of pixels assigned each DQ flag value..

pro cna_nflags, in, nflags

; Arguments:
;  in      i: input image
;  nflags  o: flag counter array (13)

  ; List of flag values to count.
  fval = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]

  ; Initialize the flag counter array.
  for k=0,13-1 do nflags(k) = 0

  ; Loop through the DQ image, counting the pixels with flags.
  tmp = where ((in.dq.data.data eq 0),cnt)
  nflags(0) = cnt
  for k=0,12-1 do begin
    tmp = where(((in.dq.data.data and fval(k)) ne 0),cnt)
    nflags(k+1) = cnt
  endfor

end


; CNA_STATS -- Compute mean, stdv, min, max of unflagged pixels
; in a NICMOS SCI image.

pro cna_stats, in, x1, x2, y1, y2, mean, medn, stdv, min, max

; Arguments:
;  in      i: input image
;  x1, x2  i: column limits for statistics calculation
;  y1, y2  i: row limits for statistics calculation
;  mean    o: mean value of unflagged pixels
;  stdv    o: standard deviation of unflagged pixels
;  min     o: minimum unflagged pixel value
;  max     o: maximum unflagged pixel value


  ; Initialize the counters and results.
  ngood = 0
  sumx  = 0
  sumx2 = 0
  mean = 0.0d0
  stdv = 0.0d0
  min  = 0.0d0
  max  = 0.0d0

  ; Compute the mean, standard deviation, min and max
  ; for the image section specified.

  ; Create an array containing only those data values corresponding to
  ; data quality values of zero.
  dqsubarray = in.dq.data.data[x1:x2,y1:y2]
  test = where((dqsubarray eq 0),count)
  if (count ne 0) then begin
    subarray = in.sci.data.data[x1:x2,y1:y2]
    testarray = subarray[test]

    ; Calculate mean,median,sd,min, and max for this masked array.
    min = min(testarray)
    max = max(testarray)
    medn = median(testarray)
    if (min eq max) then begin
      stdv = 0.0
      mean = min
    endif else begin
      ma = moment(testarray,sdev=stdv)
      mean = ma(0)
    endelse

  endif else begin

    ; All the pixels in this array are bad.
    stdv = 0.0
    mean = 0.0
    medn = 0.0
    min =  0.0
    max =  0.0

  endelse

end


; CNA_STATCALC -- Compute statistics for NICMOS science images. The min, max,
; mean and standard deviation of 'good' (unflagged) pixels are computed
; for the entire image, as well as each 128x128 quadrant. A tally of the
; number of pixels with a each DQ flag is also computed. All results are
; written to keywords in the science image extension header. This
; routine does NOT modify the input image data in any way.


function cna_statcalc, input

@calnica_common

; Arguments:
;  input  io: input image


  ; Local variables.
  nflags = lonarr(13)     ; flag counter array.
  mean = 0.0d0
  stdv = 0.0d0
  min  = 0.0d0
  max  = 0.0d0

  ; Function definitions.
  ; void cna_nflags (SingleNicmosGroup , int )
  ; void cna_stats (SingleNicmosGroup , int, int, int, int, double , 
  ;       double , double , double )
  ; int putKeyI (Hdr , char , int,   char )
  ; int putKeyF (Hdr , char , float, char )

  ; Count the number of pixels with each DQ flag value.
  cna_nflags, input, nflags

  tmpscihdr = input.sci.hdr

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'NQUAL00',nflags(0)
  sxaddpar,tmpscihdr,'NQUAL01',nflags(1)
  sxaddpar,tmpscihdr,'NQUAL02',nflags(2)
  sxaddpar,tmpscihdr,'NQUAL03',nflags(3)
  sxaddpar,tmpscihdr,'NQUAL04',nflags(4)
  sxaddpar,tmpscihdr,'NQUAL05',nflags(5)
  sxaddpar,tmpscihdr,'NQUAL06',nflags(6)
  sxaddpar,tmpscihdr,'NQUAL07',nflags(7)
  sxaddpar,tmpscihdr,'NQUAL08',nflags(8)
  sxaddpar,tmpscihdr,'NQUAL09',nflags(9)
  sxaddpar,tmpscihdr,'NQUAL10',nflags(10)
  sxaddpar,tmpscihdr,'NQUAL11',nflags(11)
  sxaddpar,tmpscihdr,'NQUAL12',nflags(12)

  ; Calculate statistics for the entire SCI image.
  cna_stats, input, 0, 255, 0, 255, mean, medn, stdv, min, max

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'GOODMEAN',mean
  sxaddpar,tmpscihdr,'GOODMEDN',medn
  sxaddpar,tmpscihdr,'GOODSTDV',stdv
  sxaddpar,tmpscihdr,'GOODMIN', min
  sxaddpar,tmpscihdr,'GOODMAX', max

  ; Calculate statistics for SCI image quadrant A.
  cna_stats, input, 0, 127, 0, 127, mean, medn, stdv, min, max

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'QAMEAN',mean
  sxaddpar,tmpscihdr,'QAMEDN',medn
  sxaddpar,tmpscihdr,'QASTDV',stdv
  sxaddpar,tmpscihdr,'QAMIN', min
  sxaddpar,tmpscihdr,'QAMAX', max

  ; Calculate statistics for SCI image quadrant B.
  cna_stats, input, 0, 127, 128, 255, mean, medn, stdv, min, max

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'QBMEAN',mean
  sxaddpar,tmpscihdr,'QBMEDN',medn
  sxaddpar,tmpscihdr,'QBSTDV',stdv
  sxaddpar,tmpscihdr,'QBMIN', min
  sxaddpar,tmpscihdr,'QBMAX', max

  ; Calculate statistics for SCI image quadrant C.
  cna_stats, input, 128, 255, 128, 255, mean, medn, stdv, min, max

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'QCMEAN',mean
  sxaddpar,tmpscihdr,'QCMEDN',medn
  sxaddpar,tmpscihdr,'QCSTDV',stdv
  sxaddpar,tmpscihdr,'QCMIN', min
  sxaddpar,tmpscihdr,'QCMAX', max

  ; Calculate statistics for SCI image quadrant D.
  cna_stats, input, 128, 255, 0, 127, mean, medn, stdv, min, max

  ; Store results in header keywords.
  sxaddpar,tmpscihdr,'QDMEAN',mean
  sxaddpar,tmpscihdr,'QDMEDN',medn
  sxaddpar,tmpscihdr,'QDSTDV',stdv
  sxaddpar,tmpscihdr,'QDMIN', min
  sxaddpar,tmpscihdr,'QDMAX', max

print,size(input.sci.hdr)
print,size(tmpscihdr)
  input.sci.hdr = tmpscihdr[0:599]

  return, 0

end

