# Memory management
define	T_mask		Memd[t_mask+$1-1]

#---------------------------------------------------------------------------
.help cf_mean Apr94 source
.ih
NAME
cf_mean -- Compute mean with rejection over a specified range.
.endhelp
#---------------------------------------------------------------------------
procedure cf_mean (flux, mask, wave, length, niter, low_reject, high_reject,
                   mean, stddev, n_points, min_wave, max_wave, new_mask)

double	flux[length]		# I:  The flux array.
double	mask[length]		# I:  The mask array.
double	wave[length]		# I:  The wavelength array.
int	length			# I:  Length of the input data.
int	niter			# I:  Number of iterations to perform.
double	low_reject		# I:  Low sigma rejection.
double	high_reject		# I:  High sigma rejection.
double	mean			# O:  The mean.
double	stddev			# O:  The error.
int	n_points		# O:  Number of points in the fit.
double	min_wave, max_wave	# O:  Minimum/maximum wavelengths in mean.
double	new_mask[length]	# O:  The new mask.

# Declarations.
bool	done			# TRUE if mean is found.
double	dx			# Generic.
int	iter			# Current iteration.
double	lcut, hcut		# Flux limits.
int	pix			# Current pixel number.
double	sum, sumsq		# Sum accumulators.
pointer	t_mask			# Temporary mask array.

begin
	# Copy mask into temporary mask.
	call malloc (t_mask, length, TY_DOUBLE)
	call amovd (mask, T_mask(1), length)

	# No do the looping.
	done = false
	iter = 0
	repeat {

            # Get starting minimum/maximum wavelengths.
            call alimd (wave, length, max_wave, min_wave)
            
            # Accumulate sums.
            sum = 0.d0
            sumsq = 0.d0
            n_points = 0
	    do pix = 1, length
		if (T_mask(pix) <= 0.d0) {
                    sum = sum + flux(pix)
                    sumsq = sumsq + flux(pix)**2
		    n_points = n_points + 1
                    if (wave(pix) < min_wave)
                        min_wave = wave(pix)
                    else if (wave(pix) > max_wave)
                        max_wave = wave(pix)
		}

            # Calculates mean and standard deviation.
            switch (n_points) {
	    case 0:
                mean = INDEFD
                stddev = INDEFD
                min_wave = INDEFD
                max_wave = INDEFD
                call error (1, "cf_mean: no points found to average")

	    case 1:
		mean = sum
		stddev = INDEFD
                done = true
		break

	    default:
		mean = sum / n_points
		dx = (sumsq - mean * sum) / (n_points - 1)
                if (dx < 0.d0) {
                    stddev = 0.d0
                    done = true
                    break
                } else
                    stddev = sqrt (dx)
            }

            # Now, remove those points outside the rejection range.
            if (IS_INDEFD(low_reject) || IS_INDEFD(high_reject) ||
		niter <= 0 || IS_INDEFI(niter)) {
                done = true
                break
            }
            lcut = mean - (low_reject * stddev)
            hcut = mean + (high_reject * stddev)
            done = true
            do pix = 1, length
                if (T_mask(pix) <= 0.d0)
                    if (flux(pix) < lcut || flux(pix) > hcut) {
                        T_mask(pix) = 1.d0
                        done = false
                    }
            if (done)
                break

	    # Increment iteration count.
	    iter = iter + 1
	} until (iter > niter)

	# If a new mask array is provided, fill it with the calculated
	# mask.
	if (!IS_INDEFD(new_mask[1]))
	    call amovd (T_mask(1), new_mask, length)
        
	# That's all folks.
	call mfree (t_mask, TY_DOUBLE)

        # If all iterations are performed without getting an answer,
        # error out.
        if (!done)
            call error (2, "cf_mean: mean does not converge")
end
#---------------------------------------------------------------------------
# end of cf_mean
#---------------------------------------------------------------------------