#---------------------------------------------------------------------------
.help ab_integrate May94 source
.ih
NAME
ab_integrate -- Integrate a function over a list of ranges.
.endhelp
#---------------------------------------------------------------------------
procedure ab_integrate (x, y, err, errflg, n_points, start, endr, n_ranges,
                        integ, integ_err, start_out, end_out)

double	x[n_points]		# I:  The x-values of the function.
double	y[n_points]		# I:  The y-values of the function.
double	err[n_points]		# I:  Absolute error in y-values.
bool	errflg			# I:  TRUE if there is error information.
int	n_points		# I:  Length of the function.
double	start[n_ranges]		# I:  The start of the integration.
double	endr[n_ranges]		# I:  The end of the integration.
int	n_ranges		# I:  Number of integrations.
double	integ[n_ranges]		# O:  Integrated function.
double	integ_err[n_ranges]	# O:  Error in integration, INDEF if none.
double	start_out[n_ranges]	# O:  Pixel where integration started.
double 	end_out[n_ranges]	# O:  Pixel where integration ended.

# Declarations.
int	i, i1, i2, j
double	r1, r2, sum, y1, y2, v1, v2, terror,dx

begin
        # Loop on points in start and end
        do i=1, n_ranges {

            # Compute indices of start(i) and endr(i)
            call ab_index(x, n_points, start[i], i1, r1)
            call ab_index(x, n_points, endr[i], i2, r2)

	    # Integrate between i1 and i2 (compute twice the integral)
	    # we will divide by two at the end.
	    sum=0.d0
	    terror = 0.d0
	    if (i2 != i1)
		do j=i1, i2-1 {
		    dx = x[j+1]-x[j]
		    sum=sum+(y[j]+y[j+1])*dx
			if (errflg)
			    terror = terror + err[j]*dx + err[j+1]*dx
		}

	    # Compute function value at r1 and r2.
	    y1 = y[i1] + (r1-i1)*(y[i1+1]-y[i1])
	    y2 = y[i2] + (r2-i2)*(y[i2+1]-y[i2])

	    # Compute integral from i1 to r1
	    v1 = (y[i1]+y1)*(start[i]-x[i1])
	    
	    # Compute integral from i2 to r2
	    v2 = (y[i2]+y2)*(endr[i]-x[i2])

	    # Integral from r1 to r2  integral[i1-i2] +
	    # integral[i2-r2] - integral[i1, r1]

	    sum=sum+v2-v1
	    integ[i]=sum/2.d0
	    start_out[i]=r1
	    end_out[i]=r2

            # Handle error propogation.
            if (errflg) {

                # terror contains the error from the integration of
                # the whole values.  Now calculate the error from
                # the partial values.
                y1 = err[i1] + (r1-i1)*(err[i1+1]-err[i1])
                y2 = err[i2] + (r2-i2)*(err[i2+1]-err[i2])

                v1 = (err[i1]+y1)*(start[i]-x[i1])
                v2 = (err[i2]+y2)*(endr[i]-x[i2])

                # Now figure complete error.
                terror = v1+v2
		integ_err[i] = terror/2.d0
                
            } else
                integ_err[i] = INDEFD
	}
end
#---------------------------------------------------------------------------
# end of ab_integrate
#---------------------------------------------------------------------------