/* 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 Ring-Moon Systems Node, March 1998
*******************************************************************************/
#include <stdio.h>
#include <string.h>
#include <math.h>
#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);
}
/*******************************************************************************
*/