include	"apercen.h"

# cdends -- obtain the mean values of the highest and lowest points of an array

procedure cdends (data, mask, npts, np_ceiling, np_floor, ceiling, floor,
			ceiling_err, floor_err)

real	data[npts]	# i: input mask array
real	mask[npts]	# i: input mask array
int	npts		# i: number of points
int	np_ceiling	# i: number of points used to determine ceiling level
int	np_floor	# i: number of points used to determine floor level
real	ceiling		# o: mean ceiling level
real	floor		# o: mean floor level
real	ceiling_err	# o: mean ceiling level's mean error
real	floor_err	# o: mean floor level's mean error

int	i, j
pointer	sp, ptr, ptmsk
real	sum, sumsq
#-------------------------------------------------------------------------------
begin
	call smark (sp)
	call salloc (ptr, npts, TY_REAL)
	call salloc (ptmsk, npts, TY_REAL)

	call amovr (data, Memr[ptr], npts)
	call amovr (mask, Memr[ptmsk], npts)

	call xt_sort2 (Memr[ptr], Memr[ptmsk], npts)

	# calculate the mean ceiling value by taking the average of the 
	# top VALID np_ceiling points
	if (np_ceiling <= 0)
	    np_ceiling = 1
	j = 0
	sum = 0.
	sumsq = 0.
	i = npts
	while (i >= 1 && j < np_ceiling) {
	    if (abs(Memr[ptmsk+i-1]-OKVAL) < IOTA) {
		sum = sum + Memr[ptr+i-1]
		sumsq = sumsq + Memr[ptr+i-1]**2
	    	j = j + 1
	    }
	    i = i - 1
	}
	if (j <= 0)
	    call error (1, "No valid ceiling points")
	else {
	    ceiling = sum / real(j)
	    if (j == 1)
		ceiling_err = sqrt(ceiling)
	    else
		ceiling_err = sqrt((sumsq/real(j) - ceiling**2)/real(j-1))
	    if (j < np_ceiling) {
		call printf ("There is only %d valid ceiling points\n")
		    call pargi (j)
	    }
	}

	# ditto for the floor level
	if (np_floor <= 0) {
	    floor = 0.
	    floor_err = 0.
	} else {
	    j = 0
	    sum = 0.
	    sumsq = 0.
	    i = 1
	    while (i <= npts && j < np_floor) {
	    	if (abs(Memr[ptmsk+i-1]-OKVAL) < IOTA) {
		    sum = sum + Memr[ptr+i-1]
		    sumsq = sumsq + Memr[ptr+i-1]**2
	    	    j = j + 1
	    	}
	    	i = i + 1
	    }
	    if (j <= 0) {
	 	floor = 0.
		floor_err = 0.
	    	call printf ("No valid floor points\n")
	    } else {
	        floor = sum / real(j)
	        if (j == 1)
		    floor_err = sqrt(floor)
	        else
		    floor_err = sqrt((sumsq/real(j) - floor**2)/real(j-1))
	        if (j < np_floor) {
		    call printf ("There is only %d valid floor points\n")
		    	call pargi (j)
	        }
	    }
	}
	call sfree (sp)
end