include	"dqm.h"		# data quality mask values

#  SUBGLOBAL -- Subtract an average of all background data
#
#  Description:
#  ------------
#  Do a low-order polynomial fit to one or two files of background (sky)
#  data -- or a subset of background data -- and subtract from the object+sky 
#  data.  There may or may not be a data mask for either the object+sky data 
#  and for the sky data.
#
#  Background data may have been taken before and after the object+sky data,
#  which is a reason why a global curve fit involving two background files
#  may be appropriate.
#  
#  Date		Author			Description
#  ----		------			-----------
#  11-Dec-1984	C. D. Biemesderfer	original
#  20-Mar-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure subglobal (indata, intime, sky, skytime, inmask, skymask, 
			outdata, outmask, order, maskflag, bmaskflag, 
			nptsdata, nptssky)

real	indata[nptsdata]	# input: input science data
double	intime[nptsdata]	# input: input data obs time
real	sky[nptssky]		# input: input background data
double	skytime[nptssky]	# input: input sky obs time
real	inmask[nptsdata]	# input: input data quality mask
real	skymask[nptssky]	# input: input background data quality mask
real	outdata[nptsdata]	# output: output science data
real	outmask[nptsdata]	# output: output data quality mask
int	order			# input: order of the fitting polynomial
bool	maskflag		# input: is there input mask file?
bool	bmaskflag		# input: is there input background mask file?
int	nptsdata		# input: number of input data points
int	nptssky			# input: number of input sky points

int	ngood			# number of good sky points 
double	chisq
real	skyvalue
int	i
pointer	sp
pointer	bgdata, bgtime, coeff
pointer	a, sig, ortho
double	tmin, tmax, normtime
double	polyval()
#==============================================================================
begin

	call smark (sp)
	call salloc (bgdata, nptssky, TY_DOUBLE)
	call salloc (bgtime, nptssky, TY_DOUBLE)
	call salloc (coeff, order+1, TY_DOUBLE)
	call salloc (a, order+1, TY_DOUBLE)
	call salloc (sig, order+1, TY_DOUBLE)
	call salloc (ortho, (order+1)**2, TY_DOUBLE)

	# Let's make sure that the observation times make some sense
	# if (skytime[1] > intime[nptsdata] || skytime[nptssky] < intime[1]) 
	    # call error (1, "No object obs in sky obs interval")

	# Initialize good bg points counter and background value
	ngood  = 0
	skyvalue = 0.

	# Collect ok sky data into one pair of arrays (obs and times)
	if (bmaskflag) {

	    # Look at background data mask
	    do i = 1 , nptssky {
		if (skymask[i] == DQM_SKY_OKVAL) {
		    ngood = ngood + 1
		    Memd[bgdata+ngood-1] = sky[i]
		    Memd[bgtime+ngood-1] = skytime[i]
		}
	    }		
    	} else {

	    # if there is no sky mask
	    ngood = nptssky
	    call achtrd (sky, Memd[bgdata], ngood)
	    call amovd (skytime, Memd[bgtime], ngood)
	}	

	# find extrema of the sky times
	call alimd (Memd[bgtime], ngood, tmin, tmax)

	# Make sure there are enough data to do a sensible fit
	if (ngood < order+2) 
	    call error (1, "Too few good sky points")

	if (tmax == tmin)
	    call error (1, "All points are oberved at the same time")

	# rescale the independent variable (time) to be within [0,1]
	do i = 1, ngood
	    Memd[bgtime+i-1] = (Memd[bgtime+i-1] - tmin) / (tmax - tmin)

	# Fit to background data using a polynomial fitting routine
	call cdpfit(Memd[bgtime], Memd[bgdata], Memd[bgdata], ngood, 0.,
			order, Memd[ortho], Memd[a], Memd[sig], chisq)

	# calculate the "regular polynomial" coefficients from the orthogonal
	# polynomial coefficients
	call cdpbas (Memd[a], Memd[sig], order, Memd[ortho], Memd[coeff])

	# Loop over observations and remove background flux
	do i = 1, nptsdata {

	    # Check input data mask, if necessary
	    if (maskflag && inmask[i] != DQM_OBJ_OKVAL) {
		outdata[i] = indata[i]
		outmask[i] = inmask[i]
	    } else {

		# Convert the obs time to a point in the interval [0, 1]
		normtime = (intime[i] - tmin) / (tmax - tmin)
		if (normtime < 0. || normtime > 1.) {
		    outdata[i] = indata[i]
		    outmask[i] = DQM_NOSKY
		}

		# Evaluate the polynomial to obtain background value
		skyvalue = real(polyval (normtime, Memd[coeff], order))

		# Correct observation
		outdata[i] = indata[i] - skyvalue
		outmask[i] = DQM_OBJ_OKVAL
	    }
	}

	call sfree (sp)

end