/* psffilt.c
********************************************************************************
* psffilt.c -- Support routines to evaluate filter coefficients.
*
* User routines
*	Pro_PSFFilter()		returns coefficients of a filter based on a
*				desired effective point spread function.
*
* Programmer routines
*	XPro_PSFInfo()		returns information about a point spread
*				function.
*	XPro_PSFInteg()		returns the integral of the product of two point
*				spread functions.
*
* NOTE: For most purposes, the absolute amplitude of the XPro_PSFInteg()
* integral is unimportant; what really matters is that its amplitude be correct
* relative to integrals calculated with alternative values of parameter b.  By
* default, this function returns values that are unnormalized and so not
* necessarily correct in absolute amplitude.  Re-compile this routine with the
* symbol ZPRO_NORMALIZED_INTEGRALS defined to produce properly normalized
* results.
*
* Mark Showalter, PDS Rings Node, March 1998
*******************************************************************************/
#include 
#include 
#include 
#include "profile.h"
#include "fortran.h"

#ifndef M_PI
#define M_PI	3.141592653589793238462643
#endif

/**********************
 * Internal constants *
 **********************/

#define SINC_TYPE	1
#define BOXCAR_TYPE	2
#define TRIANGLE_TYPE	3

/********************************
 * Internal function prototypes *
 ********************************/

static RL_FLT8 ZPro_SincSinc RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 s1, RL_FLT8 s2));
static RL_FLT8 ZPro_SincBox  RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 s1, RL_FLT8 h2));
static RL_FLT8 ZPro_SincTri  RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 s1, RL_FLT8 h2));
static RL_FLT8 ZPro_BoxSinc  RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 s2));
static RL_FLT8 ZPro_BoxBox   RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 h2));
static RL_FLT8 ZPro_BoxTri   RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 h2));
static RL_FLT8 ZPro_TriSinc  RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 s2));
static RL_FLT8 ZPro_TriBox   RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 h2));
static RL_FLT8 ZPro_TriTri   RL_PROTO((RL_FLT8 a,  RL_FLT8 b,
                                       RL_FLT8 h1, RL_FLT8 h2));

static RL_FLT8 ZPro_LinLin   RL_PROTO((RL_FLT8 xa1, RL_FLT8 ya1,
                                       RL_FLT8 xa2, RL_FLT8 ya2,
                                       RL_FLT8 xb1, RL_FLT8 yb1,
                                       RL_FLT8 xb2, RL_FLT8 yb2));

static RL_FLT8 ZPro_Si RL_PROTO((RL_FLT8 x));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_PSFFilter (psffilt.c)
*$ Abstract:
*	This function returns the set of filter coefficents needed to evaluate
*	a sample in a series given its new point spread function (PSF), based on
*	the original PSF and the transformation from old sample index to new.
*$ Keywords:
*	PROFILE, FILTER
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	RL_INT4		Pro_PSFFilter(oldpsf, newpsf, center, expand, filter,
*			              nfilter, f1, f2)
*	PRO_OBJECT	*oldpsf, *newpsf;
*	RL_FLT8		center, expand, *filter;
*	RL_INT4		nfilter, *f1, *f2;
*$ Inputs:
*	oldpsf		pointer to the PSF of the original series.
*	newpsf		pointer to the intended PSF of each new sample.
*	center		sample index in original series corresponding to
*			midpoint of new PSF.
*	expand		expansion factor from original series sampling interval
*			to new samples.
*	nfilter		size of filter array.
*$ Outputs:
*	filter[]	array of filter coefficients
*	*f1,*f2		indices in original series corresponding to first and
*			last filter coefficients returned.
*$ Returns:
*	size of new filter, f2-f1+1, or 0 on error.  If value > nfilter, then
*	the array provided is too small so the coefficients are not returned.
*$ Detailed_description:
*	This function returns the set of filter coefficents needed to evaluate
*	a sample in a series given its new point spread function (PSF), based on
*	the original PSF and the transformation from old sample index to new.
*	Note that these weights are unnormalized.
*
*	The x-axis of each PSF is assumed to be in units of the series sample
*	index.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	none
*$ Examples:
*	This snippet of code makes a filter appropriate for filtering a series
*	so that each new sample is the average of the three nearest samples in
*	another series.
*
*	PRO_OBJECT	*oldpsf, *newpsf;
*	RL_FLT8		filter[101], f1, f2;
*
*	// Create unit-width boxcar functions
*	oldpsf = Pro_BoxcarFunc(1., 0.5);
*	newpsf = Pro_BoxcarFunc(1., 0.5);
*
*	// Create filter for 3x expansion
*	Pro_PSFFilter(oldpsf, newpsf, 0., 3., filter, 101, &f1, &f2)
*
*	At this point, f1=-1, f2=1, and filter[0] = filter[1] = filter[2] != 0.
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if either PSF is NULL or is not one of the
*				allowed function types.
*	PRO_SETUP_FAILURE	if the array is too small to hold all the
*				necessary filter coefficients.
*$ Limitations:
*	Only boxcar, triangle and sinc PSFs are supported.
*
*	Coefficients when both PSFs are sinc functions are based on the
*	assumption that both sinc functions have limited domains.  This
*	approximation is acceptable when the width of both sinc functions is at
*	least several cycles.
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_INT4		Pro_PSFFilter(oldpsf, newpsf, center, expand, filter, nfilter,
		              f1, f2)
PRO_OBJECT	*oldpsf, *newpsf;
RL_FLT8		center, expand, *filter;
RL_INT4		nfilter, *f1, *f2;
{
RL_INT4		oldtype, newtype, fsize, n;
RL_FLT8		oldpar, oldhw, newpar, newhw, sumhw;

    /* Extract PSF parameters */
    oldtype = XPro_PSFInfo(oldpsf, &oldhw, &oldpar);
    if (oldtype == 0) return 0;

    newtype = XPro_PSFInfo(newpsf, &newhw, &newpar);
    if (newtype == 0) return 0;

    /* Determine filter size and check it */
    sumhw = oldhw + newhw*expand - 1.;
    *f1 = floor(center - sumhw);
    *f2 = ceil (center + sumhw);
    fsize = *f2 - *f1 + 1;

    if (fsize > nfilter) {
	(void) sprintf(xpro_message, "PSF filter setup failure\n\
array too small to hold %d coefficients", fsize);
	RL_RaiseError("PRO_SETUP_FAILURE", xpro_message);
	return fsize;
    }

/*
 * The integral we wish to return is
 *	f[n] = Integral[ oldpsf(x-(n-center)) newpsf(x/expand) dx ]
 * for
 *	*f1 <= n <= *f2
 *
 * The integral rearranges to
 *	f[n] = Integral[ oldpsf(x) newpsf((x+(n-center)/expand) dx ]
 */

    for (n = *f1; n <= *f2; n++) {
	filter[n-*f1] = XPro_PSFInteg(oldtype, newtype,
			1./expand, (n-center)/expand, oldpar, newpar);
    }

    return fsize;
}

/*
********************************************************************************
* EXPORTED PROGRAMMER ROUTINES
********************************************************************************
* XPro_PSFInfo(psf, halfwidth, param)
*
* This programmer routine checks the that given function is a valid point
* spread function and returns information about it.
*
* Input:
*	psf		point spread function object
*
* Output:
*	*halfwidth	halfwidth of function.
*	*param		parameter of function: halfwidth for boxcar and triangle
*			functions; step size for sinc functions.
*
* Return:		type of PSF: 1=sincp, 2=boxcar, 3=triangle.
*
* Errors:
*	PRO_CLASS_ERROR	if the psf is NULL or is not one of the supported
*			function types.
*******************************************************************************/

RL_INT4		XPro_PSFInfo(psf, halfwidth, param)
PRO_OBJECT	*psf;
RL_FLT8		*halfwidth, *param;
{
    if (XPro_FuncPtr(psf) == NULL) return 0;

    if (Pro_SincInfo(psf, NULL, param, halfwidth, NULL)) {
	return 1;
    }

    if (Pro_BoxcarInfo(psf, NULL, halfwidth)) {
	*param = *halfwidth;
	return 2;
    }

    if (Pro_TriangleInfo(psf, NULL, halfwidth)) {
	*param = *halfwidth;
	return 3;
    }

    XPro_ClassError("supported point spread function type", psf);
    return 0;
}

/*
********************************************************************************
* XPro_PSFInteg(type1, type2, a, b, par1, par2)
*
* This programmer function returns the integral of the product of two point
* spread functions (PSFs), after the second PSF has been translated and scaled:
*	Integral(-infinity,infinity) [psf1(x) psf2(ax+b) dx]
*
* Input:
*	type1, type2	type of each PSF: 1=sincp, 2=boxcar, 3=triangle.
*	a, b		scaling coefficients in argument to second PSF.
*	par1, par2	step size for sincp functions, halfwidth values for
*			boxcar and triangle PSFs; must be positive.
*
* Return:		value of integral, possibly unnormalized.
*
* Three different types of PSFs are supported:
*
* Type #1: sincp(x) == sin(pi*x/step) / (pi*x/step)
*
* Type #2: boxcar(x) == { 1			if -halfwidth <= x <= halfwidth.
*                       { 0			otherwise.
*
* Type #3: triangle(x) == { 1 + x/halfwidth	if -halfwidth <= x <= 0.
*                         { 1 - x/halfwidth	if 0 <= x <= halfwidth.
*                         { 0			otherwise.
*
* A handy transformation used below...
*	Let I(a,b) == Integral[ f(x) g(ax+b) dx ]
*	Let J(a,b) == Integral[ g(x) f(ax+b) dx ]
* Let u = ax+b, du = a dx, x = u/a - b/a.
*	J(a,b) = Integral[ g(u/a - b/a) f(u) du/a ]
*	       = I(1/a, -b/a) / abs(a)
*******************************************************************************/

RL_FLT8		XPro_PSFInteg(type1, type2, a, b, par1, par2)
RL_INT4		type1, type2;
RL_FLT8		a, b, par1, par2;
{
    switch (type1) {

    case SINC_TYPE:
      switch (type2) {
	case SINC_TYPE:     return ZPro_SincSinc(a, b, par1, par2);
	case BOXCAR_TYPE:   return ZPro_SincBox (a, b, par1, par2);
	case TRIANGLE_TYPE: return ZPro_SincTri (a, b, par1, par2);
	default:            return 0.;
      }

    case BOXCAR_TYPE:
      switch (type2) {
	case SINC_TYPE:     return ZPro_BoxSinc(a, b, par1, par2);
	case BOXCAR_TYPE:   return ZPro_BoxBox (a, b, par1, par2);
	case TRIANGLE_TYPE: return ZPro_BoxTri (a, b, par1, par2);
	default:            return 0.;
      }

    case TRIANGLE_TYPE:
      switch (type2) {
	case SINC_TYPE:     return ZPro_TriSinc(a, b, par1, par2);
	case BOXCAR_TYPE:   return ZPro_TriBox (a, b, par1, par2);
	case TRIANGLE_TYPE: return ZPro_TriTri (a, b, par1, par2);
	default:            return 0.;
      }
    }
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* Integrate sincp(x) * sincp(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_SincSinc(a, b, step1, step2)
RL_FLT8		a, b, step1, step2;
{
RL_FLT8		arg, sinc_arg;

/* The integral of sinc(u) sinc(au+b) du can be shown to equal
 *	pi/amax sinc(b/amax)
 * where amax = max(a,1)
 *
 * Let I(a,b) == Integral[ sincp(x/step1) sincp((ax+b)/step2) dx ]
 *	= integral[ sinc(pi/step1 * x) sinc(pi*a/step2 * x + pi*b/step2) dx ]
 *
 * Let u = pi/step1 * x.  dx = du * step1/pi.
 *    I = integral[ sinc(u) sinc(a*step1/step2 * u + pi*b/step2) du * step1/pi]
 *	= step1/pi integral[ sinc(u) sinc(a*step1/step2 * u + pi*b/step2) du ]
 *
 *	= step1/amax sinc(pi*b/step2/amax)
 * where amax = max(a*step1/step2,1)
 */

    a *= step1/step2;
    if (a < 0.) a = -a;
    if (a < 1.) a = 1.;

    arg = M_PI * b/(step2*a);
    sinc_arg = (arg == 0. ? 1. : sin(arg)/arg);

    return sinc_arg
#ifdef ZPRO_NORMALIZED_INTEGRALS
           * step1/a
#endif
           ;
}

/*
********************************************************************************
* Integrate boxcar(x) * boxcar(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_BoxBox(a, b, halfwidth1, halfwidth2)
RL_FLT8		a, b, halfwidth1, halfwidth2;
{
RL_FLT8		xlo, xhi, dx, xtemp;

    /* Find domain of second function */
    xlo = (-halfwidth2 - b) / a;
    xhi = ( halfwidth2 - b) / a;
    if (xlo > xhi) {xtemp = xlo; xlo = xhi; xhi = xtemp;}

    /* Find intersection of domains */
    if (xlo < -halfwidth1) xlo = -halfwidth1;
    if (xhi >  halfwidth1) xhi =  halfwidth1;
    dx = xhi - xlo;

    /* If domains do not intersect, return zero */
    if (dx < 0.) return 0.;

    /* Return integral, equivalent to size of intersection */
    return dx;
}

/*
********************************************************************************
* Integrate sincp(x) * boxcar(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_SincBox(a, b, step1, halfwidth2)
RL_FLT8		a, b, step1, halfwidth2;
{
RL_FLT8		xlo, xhi, xtemp, factor;

    /* Find domain of boxcar function */
    xlo = (-halfwidth2 - b) / a;
    xhi = ( halfwidth2 - b) / a;
    if (xlo > xhi) {xtemp = xlo; xlo = xhi; xhi = xtemp;}

/* The integral from 0 to x of sinc(u) du is the sine integral, Si(x).
 *
 * I = integral(xlo,xhi) [ sincp(u/step1) du ]
 *	= integral(xlo,xhi) [ sinc(pi/step1 * x) dx ]
 * Let u = pi/step1 * x
 *	= integral(pi/step1 * xlo, pi/step1 * xhi) [ sinc(u) du * step1/pi ]
 *	= step1/pi (Si(pi/step1 * xhi) - Si(pi/step1 * xlo))
 */

    factor = M_PI/step1;

    return (ZPro_Si(xhi*factor) - ZPro_Si(xlo*factor))
#ifdef ZPRO_NORMALIZED_INTEGRALS
           / factor
#endif
           ;
}

/*
********************************************************************************
* Integrate boxcar(x) * sincp(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_BoxSinc(a, b, halfwidth1, step2)
RL_FLT8		a, b, halfwidth1, step2;
{

    return ZPro_SincBox(1./a, -b/a, step2, halfwidth1)
#ifdef ZPRO_NORMALIZED_INTEGRALS
           / fabs(a)
#endif
           ;
}

/*
********************************************************************************
* Integrate triangle(x) * triangle(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_TriTri(a, b, halfwidth1, halfwidth2)
RL_FLT8		a, b, halfwidth1, halfwidth2;
{
RL_FLT8		xlo, xhi, xmid, xtemp;

    xlo = (-halfwidth2 - b) / a;
    xmid =             - b  / a;
    xhi = ( halfwidth2 - b) / a;
    if (xlo > xhi) {xtemp = xlo; xlo = xhi; xhi = xtemp;}

    return ZPro_LinLin(-halfwidth1, 0., 0., 1., xlo,  0., xmid, 1.) +
           ZPro_LinLin(-halfwidth1, 0., 0., 1., xmid, 1., xhi,  0.) +
           ZPro_LinLin(0., 1.,  halfwidth1, 0., xlo,  0., xmid, 1.) +
           ZPro_LinLin(0., 1.,  halfwidth1, 0., xmid, 1., xhi,  0.);
}

/*
********************************************************************************
* Integrate triangle(x) * boxcar(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_TriBox(a, b, halfwidth1, halfwidth2)
RL_FLT8		a, b, halfwidth1, halfwidth2;
{
RL_FLT8		xlo, xhi, xtemp;

    xlo = (-halfwidth2 - b) / a;
    xhi = ( halfwidth2 - b) / a;
    if (xlo > xhi) {xtemp = xlo; xlo = xhi; xhi = xtemp;}

    return ZPro_LinLin(-halfwidth1, 0., 0., 1., xlo, 1., xhi, 1.) +
           ZPro_LinLin(0., 1.,  halfwidth1, 0., xlo, 1., xhi, 1.);
}

/*
********************************************************************************
* Integrate boxcar(x) * triangle(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_BoxTri(a, b, halfwidth1, halfwidth2)
RL_FLT8		a, b, halfwidth1, halfwidth2;
{

    return ZPro_TriBox(1./a, -b/a, halfwidth2, halfwidth1)
#ifdef ZPRO_NORMALIZED_INTEGRALS
           / fabs(a)
#endif
           ;
}

/*
********************************************************************************
* Integrate sincp(x) * triangle(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_SincTri(a, b, step1, halfwidth2)
RL_FLT8		a, b, step1, halfwidth2;
{
RL_FLT8		ulo, umid, uhi, utemp, factor;

/* We need to evaluate the integral of a two linear functions times a sinc.
 * I(a,b) = Integral[sinc(pi/step1*x) * triangle(a*x + b) dx]
 *
 * Let u = pi/step1 * x
 * I(a,b) = Integral[sinc(u) * triangle(a*step1/pi*u + b) du * step1/pi]
 *
 * Replace the triangle function with its definition
 * I(a,b) = step1/pi * (Integral(ulo,umid)[sinc(u) * (u-ulo)/(uhi-umid) du]
 *                    - Integral(umid,uhi)[sinc(u) * (u-uhi)/(uhi-umid) du])
 * where
 * ulo   = (-hw2 - b)/a * pi/step1
 * umid  = (     - b)/a * pi/step1
 * uhi   = ( hw2 - b)/a * pi/step1
 *
 * Rearrange the two integrals so their form is identical
 * I(a,b) = (step1/pi)^2 * a/hw2 * (Integral(ulo,umid)[sinc(u) * (u-ulo) du]
 *                                + Integral(uhi,umid)[sinc(u) * (u-uhi) du])
 * I(a,b) = (J(ulo,umid) + J(uhi,umid)) * (step1/pi)^2 * a/hw2
 * where
 * J(s,t) = Integral(s,t)[sinc(u) * (u-s) du]
 *
 * Expand the new integral
 * J(s,t) = Integral(s,t)[sin(u) du] - s*Integral(s,t)[sinc(u) du]
 *        = cos(s) - cos(t) - s*Si(t) + s*Si(s)
 *
 * Back-substitute
 * I(a,b) = (cos(ulo) - cos(umid) - ulo*Si(umid) + ulo*Si(ulo) +
 *           cos(uhi) - cos(umid) - uhi*Si(umid) + uhi*Si(uhi))
 *          * (step1/pi)^2 * a/hw2
 * I(a,b) = (cos(ulo) + cos(uhi) - 2*cos(umid) +
 *           ulo*Si(ulo) + uhi*Si(uhi) - (ulo+uhi)*Si(umid))
 *          * (step1/pi)^2 * a/hw2
 */

    factor = M_PI / (a*step1);
    ulo  = (-halfwidth2 - b) * factor;
    umid =              - b  * factor;
    uhi  = ( halfwidth2 - b) * factor;

    if (ulo > uhi) {utemp = ulo; ulo = uhi; uhi = utemp;}

    return (cos(ulo) + cos(uhi) - 2.*cos(umid) +
            uhi*ZPro_Si(uhi) + ulo*ZPro_Si(ulo) - (ulo+uhi)*ZPro_Si(umid))
#ifdef ZPRO_NORMALIZED_INTEGRALS
           * step1 / (M_PI * (uhi - umid));
#endif
           ;
}

/*
********************************************************************************
* Integrate triangle(x) * sincp(ax+b) dx
*******************************************************************************/

static RL_FLT8	ZPro_TriSinc(a, b, halfwidth1, step2)
RL_FLT8		a, b, halfwidth1, step2;
{

    return ZPro_SincTri(1./a, -b/a, step2, halfwidth1)
#ifdef ZPRO_NORMALIZED_INTEGRALS
           / fabs(a)
#endif
           ;
}

/*
********************************************************************************
* Internal workhorse function to evaluate a product of two linear functions over
* their common domain.
*******************************************************************************/

static RL_FLT8	ZPro_LinLin(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2)
RL_FLT8		xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2;
{
RL_FLT8		x1, x2, slope_a, yintercept_a, slope_b, yintercept_b,
		aa, bb, cc;

    /* Find common domain */
    x1 = (xa1 > xb1 ? xa1 : xb1);
    x2 = (xa2 < xb2 ? xa2 : xb2);

    /* If it's empty, return zero */
    if (x1 >= x2) return 0.;

    /* Derive coefficients of function A: y = slope*x + intercept */
    slope_a = (ya2 - ya1) / (xa2 - xa1);
    yintercept_a = ya1 - slope_a * xa1;

    /* Derive coefficients of function B: y = slope*x + intercept */
    slope_b = (yb2 - yb1) / (xb2 - xb1);
    yintercept_b = yb1 - slope_b * xb1;

    /* Derive coefficients of integral of product: ax^3 + bx^2 + cx */
    aa = slope_a * slope_b / 3.;
    bb = (slope_a * yintercept_b + slope_b * yintercept_a) / 2.;
    cc = yintercept_a * yintercept_b;

    /* Evaluate definite integral */
    return ((aa*x2 + bb)*x2 + cc)*x2 - ((aa*x1 + bb)*x1 + cc)*x1;
}

/*
********************************************************************************
* ZPro_Si(x) -- Approximate sine integral
*
* Si(x) == Integral{0...x} [ sin(u)/u du ]
*
* This routine uses a fast but approximate expression for x > 1.  It should be
* accurate to 6 digits, good enough for most purposes.
*******************************************************************************/
#define CUTOFF	1.

static RL_FLT8	ZPro_Si(x)
RL_FLT8		x;
{
RL_FLT8		value, absx, x2, term, sum, prev, fx, gx;
RL_INT4		j;

RL_FLT8	a1 = 38.027264, a2 = 265.187033, a3 =  335.677320, a4 =  38.102495;
RL_FLT8	b1 = 40.021433, b2 = 322.624911, b3 =  570.236280, b4 = 157.105423;
RL_FLT8	c1 = 42.242855, c2 = 302.757865, c3 =  352.018498, c4 =  21.821899;
RL_FLT8	d1 = 48.196927, d2 = 482.485984, d3 = 1114.978885, d4 = 449.690326;

	absx = fabs(x);
	x2 = x * x;

/*******************************************************************************
* Handle the easy case first
*******************************************************************************/

	if (x == 0.) {
	    value = 0.;
	}

/*******************************************************************************
* Small x solution: si(x) = x - x^3/3*3! + x^5/5*5! - x^7/7*7! ...
*
* In the algorithm below, term = (-1)^((j-1)/2) x^j / j! at each step:
*	term(1) =  x
*	term(3) = -x^3 / 3!
*	term(5) =  x^5 / 5!
* The series is a sum of term(j)/j.  The value returned is exact.
*******************************************************************************/

	else if (absx <= CUTOFF) {
	    j    = 1;
	    term = x;
	    sum  = x;
	    do {
		j += 2;
		term *= x2 / (double) ((1-j) * j);

		prev = sum;
		sum += term / j;
	    } while (prev != sum);

	    value = sum;
	}

/*******************************************************************************
* Large x solution: si(x) = PI/2 - f(x)cos(x) -  g(x)sin(x)
*
* where f(x) and g(x) are auxiliary functions defined by
*	Abramowitz, M., and I. A. Stegun 1966.  Handbook of Mathematical
*	Functions.  National Bureau of Standards Applied Mathematics Series
*	55, pp. 231--233.
* See Section 5.2.
*
* Expressions for f(x) and g(x) have errors of magnitude < 5e-7 for x >= 1.
*******************************************************************************/

	else {

	    /* Evaluate approximate auxiliary functions */
	    fx =  ((((x2 + a1)*x2 + a2)*x2 + a3)*x2 + a4) /
	         (((((x2 + b1)*x2 + b2)*x2 + b3)*x2 + b4)*absx);

	    gx =  ((((x2 + c1)*x2 + c2)*x2 + c3)*x2 + c4) /
	         (((((x2 + d1)*x2 + d2)*x2 + d3)*x2 + d4)*x2);

	    /* Solve for sine integral */
	    value = M_PI/2. - fx*cos(absx) - gx*sin(absx);

	    if (x < 0.) value = -value;
	}

	return value;
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
*	FPro_PSFFilter (psffilt.c)
*$ Abstract:
*	This function returns the set of filter coefficents needed to evaluate
*	a sample in a series given its new point spread function (PSF), based on
*	the original PSF and the transformation from old sample index to new.
*$ Keywords:
*	PROFILE, FILTER
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	integer*4 function FPro_PSFFilter(oldpsf, newpsf, center, expand, filter,
*                                   nfilter, f1, f2)
*	integer*4	oldpsf, newpsf, nfilter, f1, f2
*	real*8		center, expand, filter(*)
*$ Inputs:
*	oldpsf		FORTRAN pointer to the PSF of original series.
*	newpsf		FORTRAN pointer to the intended PSF of each new sample.
*	center		sample index in original series corresponding to
*			midpoint of new PSF.
*	expand		expansion factor from original series sampling interval
*			to new samples.
*	nfilter		size of filter array.
*$ Outputs:
*	filter(*)	array of filter coefficients
*	f1, f2		indices in original series corresponding to first and
*			last filter coefficients returned.
*$ Returns:
*	size of new filter, f2-f1+1, or 0 on error.  If value > nfilter, then
*	the array provided is too small so the coefficients are not returned.
*$ Detailed_description:
*	This function returns the set of filter coefficents needed to evaluate
*	a sample in a series given its new point spread function (PSF), based on
*	the original PSF and the transformation from old sample index to new.
*	Note that these weights are unnormalized.
*
*	The x-axis of each PSF is assumed to be in units of the series sample
*	index.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	none
*$ Examples:
*	This snippet of code makes a filter appropriate for filtering a series
*	so that each new sample is the average of the three nearest samples in
*	another series.
*
*	integer*4	oldpsf, newpsf
*	real*8		filter(101), f1, f2
*
*	! Create unit-width boxcar functions
*	oldpsf = FPro_BoxcarFunc(1.d0, 0.5d0)
*	newpsf = FPro_BoxcarFunc(1.d0, 0.5d0)
*
*	! Create filter for 3x expansion
*	call Pro_PSFFilter(oldpsf, newpsf, 0.d0, 3.d0, filter, 101, f1, f2)
*
*	At this point, f1=-1, f2=1, and filter(1) = filter(2) = filter(3) != 0.
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if either PSF is NULL or is not one of the
*				allowed function types.
*	PRO_SETUP_FAILURE	if the array is too small to hold all the
*				necessary filter coefficients.
*	FORTRAN_POINTER_ERROR	if any object is not a valid FORTRAN pointer.
*$ Limitations:
*	Only boxcar, triangle and sinc PSFs are supported.
*
*	Coefficients when both PSFs are sinc functions are based on the
*	assumption that both sinc functions have limited domains.  This
*	approximation is acceptable when the width of both sinc functions is at
*	least several cycles.
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_INT4	FORTRAN_NAME(fpro_psffilter) (oldpsf, newpsf, center, expand, filter,
					nfilter, f1, f2)
RL_INT4	*oldpsf, *newpsf, *nfilter, *f1, *f2;
RL_FLT8	*center, *expand, *filter;
{
RL_VOID *oldptr, *newptr;

    /* Look up function pointers */
    oldptr = FORT_GetPointer(*oldpsf);
    if (oldptr == NULL) return 0;

    newptr = FORT_GetPointer(*newpsf);
    if (newptr == NULL) return 0;

    /* call function */
    return Pro_PSFFilter((PRO_OBJECT *) oldptr, (PRO_OBJECT *) newptr,
		*center, *expand, filter, *nfilter, f1, f2);
}

/*******************************************************************************
*/