define	OKVAL	0

#  PHASE_COADD -- Sort data into bins and accumulate the sums in each bin.
#
#  Description:
#  ------------
#  Accumulate the weighted data pixel values and their squares in each phase 
#  bin.  If there is an input mask file, the "bad" pixels will not be included 
#  in the accumulation.
#
#  The scheme of binning is illustrated by the following example:  if there are
#  5 bins between phases 0 and 1, the points in [0, 0.1) and in [0.9, 1) go to
#  bin #1, points in [0.1, 0.3) go to bin #2, points in [0.3, 0.5) go to bin #3,
#  etc.
#  
#  Date		Author			Description
#  ----		------			-----------
#  20-Feb-1990  J.-C. Hsu		design and code
#------------------------------------------------------------------------------

procedure phase_coadd (arrphase, arrdata, arrmask, maskflag, npix, nbins,
			npts, sum, sumwt, sumsq)

double	arrphase[npix]	# input: phase array
real	arrdata[npix]	# input: data array
real	arrmask[npix]	# input: mask array
bool	maskflag	# input: is there input mask file?
int	npix		# input: size of input array
int	nbins		# input: number of (output) bins
int	npts[nbins]	# input/output: number of data points in each bin
double	sum[nbins]	# input/output: accumulated data values in each bin
double	sumwt[nbins]	# input/output: accumulated weight in each bin
double	sumsq[nbins]	# input/output: accumulated data squares in each bin

double	weight		# weight of the data point
int	indx		# bin index
int	i
#==============================================================================
begin

	# put each data point into the apropriate bin
	do i = 1, npix {
	    if (maskflag && arrmask[i] != OKVAL) 
		;
	    else {
		indx = int(arrphase[i] * nbins + 0.5d0) + 1
		if (indx == nbins+1)
		    indx = 1
		if (indx < 1 || indx > nbins) {
		    call printf ("point no. %s out of range: %g\n")
		    	call pargi(i)
			call pargd(arrphase[i])
		} else {
		    npts[indx] = npts[indx] + 1

		    # use equal weights now since there is no input variance
		    # and input count(rate) might be zero
		    weight = 1.d0 
		    sumwt[indx] = sumwt[indx] + weight
		    sum[indx] = sum[indx] + weight * arrdata[i]
		    sumsq[indx] = sumsq[indx] + weight * arrdata[i]**2
		}
	    }
	}
end