# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<imhdr.h>
include	"../icombine.h"

$for (sird)
# IC_AVERAGE -- Compute the average (or summed) image line.
# Options include a weighted average/sum.

procedure ic_average$t (d, m, n, wts, npts, doblank, doaverage, average)

pointer	d[ARB]			# Data pointers
pointer	m[ARB]			# Image ID pointers
int	n[npts]			# Number of points
real	wts[ARB]		# Weights
int	npts			# Number of output points per line
int	doblank			# Set blank values?
int	doaverage		# Do average?
$if (datatype == sil)
real	average[npts]		# Average (returned)
$else
PIXEL	average[npts]		# Average (returned)
$endif

int	i, j, k
real	sumwt, wt
$if (datatype == sil)
real	sum
$else
PIXEL	sum
$endif

include	"../icombine.com"

begin
	# If no data has been excluded do the average/sum without checking
	# the number of points and using the fact that the weights are
	# normalized.  If all the data has been excluded set the average/sum
	# to the blank value if requested.

	if (dflag == D_ALL) {
	    if (dowts) {
		do i = 1, npts {
		    k = i - 1
		    wt = wts[Memi[m[1]+k]]
		    sum = Mem$t[d[1]+k] * wt
		    do j = 2, n[i] {
			wt = wts[Memi[m[j]+k]]
			sum = sum + Mem$t[d[j]+k] * wt
		    }
		    average[i] = sum
		}
	    } else {
		do i = 1, npts {
		    k = i - 1
		    sum = Mem$t[d[1]+k]
		    do j = 2, n[i]
			sum = sum + Mem$t[d[j]+k]
		    if (doaverage == YES)
			average[i] = sum / n[i]
		    else
			average[i] = sum
		}
	    }
	} else if (dflag == D_NONE) {
	    if (doblank == YES) {
		do i = 1, npts
		    average[i] = blank
	    }
	} else {
	    if (dowts) {
		do i = 1, npts {
		    if (n[i] > 0) {
			k = i - 1
			wt = wts[Memi[m[1]+k]]
			sum = Mem$t[d[1]+k] * wt
			sumwt = wt
			do j = 2, n[i] {
			    wt = wts[Memi[m[j]+k]]
			    sum = sum + Mem$t[d[j]+k] * wt
			    sumwt = sumwt + wt
			}
			if (doaverage == YES)
			    average[i] = sum / sumwt
			else
			    average[i] = sum
		    } else if (doblank == YES)
			average[i] = blank
		}
	    } else {
		do i = 1, npts {
		    if (n[i] > 0) {
			k = i - 1
			sum = Mem$t[d[1]+k]
			do j = 2, n[i]
			    sum = sum + Mem$t[d[j]+k]
			if (doaverage == YES)
			    average[i] = sum / n[i]
			else
			    average[i] = sum
		    } else if (doblank == YES)
			average[i] = blank
		}
	    }
	}
end
$endfor