include	"dqm.h"		# data quality mask values
define	SZ_FITTYPE	8

#  FITNSUB -- Subtract an interpolated value of sky from object + sky data
#
#  Description:
#  ------------
#  Subtract sky observations from object + sky observations taking data masks
#  into consideration and producing a data mask for the resulting data.
#  The sky observations which are subtracted are obtained by interpolation
#  or by curve fitting.
#  
#  Date		Author			Description
#  ----		------			-----------
#  11-Dec-1984	C. D. Biemesderfer	original
#  20-Mar-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure fitnsub (indata, intime, sky, skytime, inmask, skymask, 
			outdata, outmask, fittype, nptsfit, 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
char	fittype[SZ_FITTYPE]	# input: fitting scheme
int	nptsfit			# input: number of fitting points
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 
int	first, last
real	skyvalue
int	i, j
bool	streq()
#==============================================================================
begin

	# make sure it is simultaneous or interleaved observation
	if (nptsdata != nptssky)
	    call error (1, "object and sky do not have same number of pixels")

	# 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 {

		# Find sky observations closest to the data point
		# since this routine only deals with the scheme other than 
		# OVERALL, the sky obs must be simultaneous or interleaved, i.e.
		# there is a corresponding sky pixel for each data pixel
		
		# find the first pixel, consider both cases of odd and even 
		# fitting points
		first = i - nptsfit/2
		last = first + nptsfit - 1

		# shift fitting interval right if falls off left and vice versa
		if (first < 1) {
		    first = 1
		    last = nptsfit
		} else if (last > nptssky) {
		    first = nptssky - nptsfit + 1
		    last = nptssky
		}

		# check sky mask to make sure there is enough sky points
		ngood = 0
		do j = first , last {
		    if (bmaskflag && skymask[j] != DQM_SKY_OKVAL)
			;
		    else 
			ngood = ngood + 1
		}

		# least square fitting case:
		if (streq(fittype, "leastsq")) {
		    if (ngood < order+2) {
		        outdata[i] = indata[i]
		        outmask[i] = DQM_NOSKY
		    } else 
			call subglobal(indata[i], intime[i], sky[first],
				skytime[first], inmask[i], skymask[first],
				outdata[i], outmask[i], order, maskflag, 
				bmaskflag, 1, nptsfit)

		} else {
		    if (ngood < nptsfit) {
		        outdata[i] = indata[i]
		        outmask[i] = DQM_NOSKY
		    } else {

			# nearest neighbor case:
			if (streq(fittype, "nearest")) {
			    outdata[i] = indata[i] - sky[i]
			    outmask[i] = DQM_OBJ_OKVAL

			# linear fitting case:
			} else if (streq(fittype, "linear")) {
			    skyvalue = (intime[i] - skytime[first]) *
				(sky[last] - sky[first]) / 
				(skytime[last] - skytime[first])
			    outdata[i] = indata[i] - skyvalue
			    outmask[i] = DQM_OBJ_OKVAL
			}
		    }
		}
	    }
	}

end