C+
C NAME:
C	smei_frm_ped
C PURPOSE:
C	Calculate pedestal for SMEI data frame
C CATEGORY:
C CALLING SEQUENCE:
	integer function smei_frm_ped(cframe,nx,ny,frame,pattern,ped_val,ped_aux)
C INPUTS:
C	cframe		character*(*)	frame name (contains time and camera)
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		"shutter-closed" calibration pattern
C					The pattern is only used if ped_aux(SMEI__PED_RIGHT_OFFSET)
C					is NOT bad, in which case it is subtracted from the
C					pedestal columns prior to updating the left-right offset
C					stored in this element of ped_aux.
C					This is only used for camera 3 when both sides of the frame
C					are used for the pedestal determination.
C	ped_aux(*)	real		six entries:
C
C					ped_aux(SMEI__PED_RIGHT_OFFSET)
C					    Offset between mean left and right pedestal. If not
C					    set to bad than naive means are calculated after taking
C					    out the pattern and this left-to-right offset.
C					    If set to bad, then the pattern array is not accessed.
C
C					ped_aux(SMEI__PED_MEAN_REF    )  mean 'reference' pedestal
C					ped_aux(SMEI__PED_MEAN_EXCESS )  mean excess
C					    frames with naive mean above the threshold
C						ped_aux(SMEI__PED_MEAN_REF)+ped_aux(SMEI__PED_MEAN_EXCESS)
C					    are rejected.
C					    If the right-to-left offset ped_aux(SMEI__PED_RIGHT_OFFSET)
C					    is not bad then the means are calculated after taking out
C					    the pattern and the right-to-left offset.
C
C					ped_aux(SMEI__PED_MEDIAN_DEFICIT)
C					ped_aux(SMEI__PED_MEDIAN_EXCESS )
C					    only pedestal pixels with values between
C					    the naive median pedestal+median_deficit and
C					    naive median pedestal+median_excess.
C					    are used to get the final pedestal.
C					    If set to BadR4() the corresponding
C					    restriction is not applied
C OUTPUTS:
C	smei_frm_ped	integer		# pixels used in final pedestal calculation
C					If only the left side of the frame is
C					used for the pedestal calculation, i.e. if 
C					ped_aux(SMEI__PED_LEFT_SIDE_ONLY) .ne. bad,
C					then the result is multiplied by 2 before
C					returning.
C	ped_val		real		pedestal
C					value will be bad if
C					- no valid pedestal pixels (smei_frm_ped=0)
C					- mean pedestal higher than ped_mean_ref+ped_mean_excess
C	ped_aux(*)	real		six entries:
C					ped_aux(SMEI__PED_RIGHT_OFFSET	)  updated right-to-left offset
C					    Updated only if input was not bad. Update depends on input
C					    value of ped_aux(SMEI__PED_MEAN_EXCESS).
C					ped_aux(SMEI__PED_LEFT_SIDE_ONLY)  use left ped cols only
C					    used only if ped_aux(SMEI__PED_RIGHT_OFFSET) is set to bad
C					    (no pattern used).
C					ped_aux(SMEI__PED_NAIVE_MEAN	)  naive mean pedestal
C					ped_aux(SMEI__PED_NAIVE_MEDIAN	)  naive pedestal median
C					ped_aux(SMEI__PED_MEDIAN_LOW	)  naive median+median_deficit
C					ped_aux(SMEI__PED_MEDIAN_HIGH	)  naive median+median_excess
C					ped_aux(SMEI__PED_VERY_LOW	)  ??
C					ped_aux(SMEI__PED_STDEV		)  standard dev. of ped_val
C CALLS:
C	BadR4, ArrR4Mean, SortR4, smei_frm_pickup, ArrR4MinusArrR4,
C	ArrR4AddConstant, ArrR4Stdev, ArrR4Median, Say
C	smei_frm_ped_mean
C SEE ALSO:
C	smei_frm_clean
C RESTRICTIONS:
C	The median deficit must be negative, the median excess must
C	be positive. If not, all bets are off.
C INCLUDE:
	include		'smei_frm_layout.h'
	include		'smei_frm_basepar.h'
C PROCEDURE:
C	SMEI__* parameters are defined in the two include files.
C
C	Better check the code to find out what 'naive' means
C	these days.
C
C	Only the pedestal values in the frame are accessed (4
C	columns left and right in engineering mode)
C
C	First the 'naive' mean pedestal (average over all pedestal
C	pixels with valid values) is calculated and compared with
C	ped_mean_ref+ped_mean_excess. If the 'naive' mean is higher
C	than ped_mean_ref+ped_mean_excess the pedestal is rejected.
C
C	If this first test is passed then the median pedestal is
C	calculated. All pedestal pixels more than ped_median_excess
C	above the median are rejected. (the median+ped_median_excess
C	is returned in ped_median_high).
C
C	The mean of the remaining pixels is returned as the
C	pedestal ped_val
C MODIFICATION HISTORY:
C	JUN-2004, Paul Hick (UCSD/CASS)
C	FEB-2007, Paul Hick (UCSD/CASS)
C	    Added option to use only left side of frame to calculate
C	    the pedestal when the pattern is not used.
C	NOV-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Added argument cframe.
C	    Mean pedestal is now calculated with smei_frm_ped_mean.
C	    A constant value (hardcoded in smei_base and smei_cal)
C	    was passed in as ped_aux(SMEI__PED_MEAN_REF).
C	    This is not accessed anymore.
C-
	    character	cframe*(*)
	    integer	nx
	    integer	ny
	    real	frame  (*)
	    real	pattern(*)
	    real	ped_val
	    real	ped_aux(*)

	character	cSay*7			/'frm_ped'/
	character	cFlt2Str*14
	integer		smei_frm_pickup
	real		ped_pixs(SMEI__PED_NPIX)		! Scratch array
	real		ped_patt(SMEI__PED_NPIX)		! Scratch array

	ped_right_offset   = ped_aux(SMEI__PED_RIGHT_OFFSET)
	ped_mean_excess	   = ped_aux(SMEI__PED_MEAN_EXCESS )
	!ped_threshold      = ped_aux(SMEI__PED_MEAN_REF    )+ped_mean_excess
	ped_threshold      = smei_frm_ped_mean(cframe)+ped_mean_excess
	ped_median_deficit = ped_aux(SMEI__PED_MEDIAN_DEFICIT)
	ped_median_excess  = ped_aux(SMEI__PED_MEDIAN_EXCESS )
	ped_left_side_only = ped_aux(SMEI__PED_LEFT_SIDE_ONLY)

	bad = BadR4()

	if (ped_right_offset .eq. bad) then			! Don't use pattern

	    if (ped_left_side_only .eq. bad) then
		n = smei_frm_pickup(0+0, nx, ny, frame, ped_pixs)! Left and right side
	    else
		n = smei_frm_pickup(0+4, nx, ny, frame, ped_pixs)! Left side only
	    end if

	else							! Use pattern

	    ! Store left side in bottom half of ped_pix
	    ! Calculate left side mean after taking out pattern

	    n = smei_frm_pickup(0+4, nx, ny, pattern, ped_patt)	! Left side
	    n = smei_frm_pickup(0+4, nx, ny, frame  , ped_pixs)

	    call ArrR4MinusArrR4(-n,ped_pixs,ped_patt,ped_pixs)	! Remove pattern

	    ped_mean_left = ArrR4Mean(-n, ped_pixs , nleft)

	    ! Store right side in top half of ped_pixs

	    nhalf  = n		! Should be SMEI__PED_NPIX/nbin/2
	    iright = n+1

	    ! Calculate right side mean after taking out pattern.
	    ! Also take out ped_right_offset

	    n = smei_frm_pickup(0+8, nx, ny, pattern, ped_patt)	! Right side
	    n = smei_frm_pickup(0+8, nx, ny, frame  , ped_pixs(iright))

	    call ArrR4MinusArrR4 (-n,ped_pixs(iright), ped_patt		,ped_pixs(iright))
	    call ArrR4AddConstant(-n,ped_pixs(iright),-ped_right_offset ,ped_pixs(iright))

	    ped_mean_right = ArrR4Mean(-n, ped_pixs(iright), nright)

	    n = nhalf+n		! Should be SMEI__PED_NPIX/nbin

	end if

	ped_naive_mean = ArrR4Mean(-n, ped_pixs, smei_frm_ped)
	ped_very_low   = 0.0

	! Probably smei_frm_ped is always SMEI__PED_NPIX-4 (=2044) for mode 0.
	! so smei_frm_ped = 0 (all ped_pixs bad) never happens.
	
	if (smei_frm_ped .eq. 0 .or. ped_naive_mean .gt. ped_threshold) then

	    ped_naive_median = bad
	    ped_median_low   = bad
	    ped_median_high  = bad
	    ped_stdev	     = bad

	    ped_val	     = bad

	else

	    if (ped_right_offset .ne. bad) then		! Use pattern

		! Use difference right and left to update ped_right_offset

		ped_mean_diff = ped_mean_right-ped_mean_left

		if (abs(ped_mean_diff) .lt. ped_mean_excess) then
		    ped_right_offset = ped_right_offset+0.1*ped_mean_diff
		else
		    call Say(cSay,'W','L2R offset','not updated; left-right difference of mean is '//cFlt2Str(ped_mean_diff,2))
		    !ped_right_offset = 0.0
		end if

		ped_aux(SMEI__PED_RIGHT_OFFSET) = ped_right_offset 

	    end if

	    ped_naive_median = ArrR4Median(-n, ped_pixs, bad, j, ped_pixs)

	    ! Elements 1..smei_frm_ped in ped_pixs are now sorted.
	    ! Values above smei_frm_ped are bad.

	    if (j .ne. smei_frm_ped) stop 'smei_frm_pedestal: shoot'

	    ilow  = 1
	    ihigh = smei_frm_ped

	    if (ped_median_deficit .ne. bad) then	! Used to exclude 'white measles'

		ped_median_low = ped_naive_median+ped_median_deficit

		do while (ped_pixs(ilow) .le. ped_median_low)
		    ilow = ilow+1			! Find pixel > ped_median_low
		end do

		if (max(ped_pixs(1), ped_pixs(2)) .lt. ped_median_low) then
		    ped_very_low = 1.0
		    !print *, 'BAD PED=',ped_pixs(1), ped_pixs(2), ped_naive_median, ped_median_deficit, ped_aux(SMEI__PED_NAIVE_MEDIAN)
		end if

	    else
		ped_median_low = bad

	    end if

	    if (ped_median_excess .ne. bad) then

		ped_median_high = ped_naive_median+ped_median_excess

		do while (ped_pixs(ihigh) .gt. ped_median_high)
		    ihigh = ihigh-1			! Find pixel <= ped_median_high
		end do

	    else
		ped_median_high = bad

	    end if

	    ! Only include pixels between ped_median_high and ped_median_high
	    ! for final pedestal calculation

	    ihigh = ihigh-ilow+1

	    ped_val   = ArrR4Mean (  -ihigh,ped_pixs(ilow),smei_frm_ped)
	    ped_stdev = ArrR4Stdev(1,-ihigh,ped_pixs(ilow),ped_val,n)

	    if (smei_frm_ped .ne. ihigh) stop 'smei_frm_pedestal: shoot too'

	end if

	ped_aux(SMEI__PED_NAIVE_MEAN	) = ped_naive_mean
	ped_aux(SMEI__PED_VERY_LOW	) = ped_very_low
	ped_aux(SMEI__PED_NAIVE_MEDIAN	) = ped_naive_median
	ped_aux(SMEI__PED_MEDIAN_LOW	) = ped_median_low
	ped_aux(SMEI__PED_MEDIAN_HIGH	) = ped_median_high
	ped_aux(SMEI__PED_STDEV		) = ped_stdev

	if (ped_right_offset .eq. bad) then			! Pattern not used
	    if (ped_left_side_only .ne. bad) smei_frm_ped = 2*smei_frm_ped
	end if

	return
	end