C+
C NAME:
C	smei_frm_dark
C PURPOSE:
C	Calculates dark current for SMEI data frame
C CATEGORY:
C	ucsd/camera/for/lib
C CALLING SEQUENCE:
	integer function smei_frm_dark(bMedian, nx, ny, frame, pattern, dark_val, dark_pixs, dark_aux)
C INPUTS:
C	bMedian		logical		.TRUE. : return median
C					.FALSE.: return mean
C	nx		integer		horizontal frame size
C	ny		integer		vertical frame size
C	frame(nx,ny)	real		frame data (after massaging by href=smei_frm_clean=)
C	pattern(nx,ny)	real		dark current pattern
C					    (not accessed if dark_aux(SMEI__DRK_RATIO) is bad)
C	dark_aux(*)	real		dark_aux(SMEI__DRK_PEDESTAL)
C					    pedestal (return value of href=smei_frm_ped=)
C					dark_aux(SMEI__DRK_HEADROOM)
C					    used in median calculation of dark current
C
C					dark_aux(SMEI__DRK_RATIO)
C					    ratio of dark current in frame and dark current
C					    in pattern. If set bad, then the pattern array is
C					    not used. The caller sets the ratio by using some
C					     estimate for the dark current in the frame.
C					dark_aux(SMEI__DRK_CUT)
C					    threshold above pedestal used to reject dark
C					    current pedestal.
C					    (not accessed if dark_aux(SMEI__DRK_RATIO) is bad)
C					dark_aux(SMEI__DRK_POWER)
C					    not used (yet?)
C OUTPUTS:
C	smei_frm_dark	real		# pixels contributing to dark current
C	dark_val	real		dark current value (with pedestal subtracted)
C					    if dark_val is bad then smei_frm_dark=0 and v.v.
C	dark_pixs(*)	real		all dark current pixel values (with pedestal subtracted)
C					    used to calculate dark_val. If the pattern was
C					    used than 'cosmic ray' pixels were flagged as bad.
C	dark_aux(*)	real		three entries
C					dark_aux(SMEI__DRK_MEAN  ) mean-based dark current
C					dark_aux(SMEI__DRK_MEDIAN) median-based dark current
C					dark_aux(SMEI__DRK_STDEV ) standard deviation of ped_val
C INCLUDE:
	include		'smei_frm_layout.h'
	include		'smei_frm_basepar.h'
C CALLS:
C	BadR4, ArrR4AddConstant, ArrR4Median, ArrR4Mean, ArrR4Stdev, smei_frm_pickup
C RESTRICTIONS:
C	The values in the pattern (when used) should all be valid,
C	i.e. no BadR4() values at least for those pixels in the
C	input frame which are expected to contain valid data.
C SEE ALSO:
C	smei_frm_clean, smei_frm_ped
C PROCEDURE:
C	Only the dark current values in the frame are accessed (4 columns left and right
C	in engineering mode).
C
C	If dark_aux(SMEI__DRK_RATIO) is bad then the return value 'dark_val' is the median
C	or the mean (depending on setting of bMedian) of the dark current pixels, with the
C	pedestal value in dark_aux(SMEI__DRK_PEDESTAL) subtracted.
C
C	If the pattern and dark ratio are available (dark_aux(SMEI__DRK_RATIO) is NOT bad)
C	then the best estimate of the expected dark current is
C	   dark_aux(SMEI__DRK_PEDESTAL)+dark_aux(SMEI__DRK_RATIO)*pattern
C	
C	The threshold dark_aux(SMEI__DRK_CUT) is used to detect spikes ('cosmic rays'):
C	by rejected dark current pixels with values above
C	   dark_aux(SMEI__DRK_PEDESTAL)+dark_aux(SMEI__DRK_RATIO)*pattern+dark_aux(SMEI__DRK_CUT).
C	The median (or mean) is then calculated after the spikes have been flagged as bad.
C
C	Note that the pattern is multiplied with the dark ratio, and NOT with some power of the
C	dark ratio as in ???.
C MODIFICATION HISTORY:
C	JUN-2004, Paul Hick (UCSD/CASS)
C	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Both mean and median are now calculated and returned in
C	    dark_aux(SMEI__DRK_MEAN) and dark_aux(SMEI__DRK_MEDIAN).
C	    The return value dark_val is one of these two, depending
C	    on setting of bMedian.
C-
	    logical	bMedian
	    integer	nx
	    integer	ny
	    real	frame  (nx,ny)
	    real	pattern(nx,ny)
	    real	dark_val
	    real	dark_pixs(*)
	    real	dark_aux (*)

	character	cSay*8			/'frm_dark'/
	integer		smei_frm_pickup
	real		dark_tmp(SMEI__DRK_NPIX)	! Scratch array

	bad = BadR4()

	! Pick up dark current pixels left and right (may contain bad values)

	n = smei_frm_pickup(1, nx, ny, frame, dark_pixs)

	ped_val    = dark_aux(SMEI__DRK_PEDESTAL)
	dark_ratio = dark_aux(SMEI__DRK_RATIO   )

	if (dark_ratio .ne. bad) then

	    ! This section is not accessed when the pattern is made, but only
	    ! for the measles/cosmic ray removal when the pattern is known already

	    ! Pick up pattern in dark current areas left and right
	    ! These pattern values should all be valid (except for four
	    ! pixels in top row at right of CCD).

	    n = smei_frm_pickup(1, nx, ny, pattern, dark_tmp)

	    ! the CR removal differs for first 10; OK here, no CRs found!
	    ! CR cut for the pattern, = 200 e-

	    dark_threshold = ped_val+dark_aux(SMEI__DRK_CUT)
	    !power = dark_aux(SMEI__DRK_POWER)
	    !if (power .ne. 1.0) dark_ratio = dark_ratio**power

	    do i=1,n

		val = dark_pixs(i)

		if (val .ne. bad) then
		    if (val .gt. dark_ratio*dark_tmp(i)+dark_threshold) val = bad
		end if

		dark_pixs(i) = val

	    end do

	end if
						! Subtract pedestal
	call ArrR4AddConstant(-n, dark_pixs, -ped_val, dark_pixs)

	! Calculate mean and median of remaining bins

	wbin = 1.0/dark_aux(SMEI__DRK_HEADROOM)
	dark_aux(SMEI__DRK_MEDIAN) = ArrR4Median(-n,dark_pixs,wbin,smei_frm_dark,dark_tmp)
	dark_aux(SMEI__DRK_MEAN  ) = ArrR4Mean  (-n,dark_pixs,     smei_frm_dark)

	if (bMedian) then			! Return median
	    dark_val = dark_aux(SMEI__DRK_MEDIAN)
	else					! Return mean
	    dark_val = dark_aux(SMEI__DRK_MEAN  )
	end if

	dark_aux(SMEI__DRK_STDEV) = ArrR4Stdev(1,-n,dark_pixs,dark_val,i)

	return
	end