TABLE OF CONTENTS

CsplineUtils/cSplineConvolution [ Functions ]

NAME

    cSplineConvolution --- compute the convolution of two basis functions

FUNCTION

This function computes the convolution of two spline basis functions, that is, 

       int_0^infty ( x^k B_1(x) * B_2(x) ) dx

where the splines may be of different orders, but are defined on the same set of knots. See also splineconv.

SYNOPSIS

341 double cSplineConvolution(int k, int n1,int j1, int n2, int j2, double *knots, int splord)

INPUTS

    k      the power of x used in the convolution
    n1     order of the first spline
    j1     largest knot of the first spline
    n2     order of the second spline
    j2     largest knot of the second spline
    knots  set of knots on which the splines are defined
    splord order of the splines

OUTPUTS

    The k-th order convolution of the splines defined by the input parameters

SOURCE

345 {
346     double out=0;
347     // if the splines don't overlap, the convolution is 0
348     if((j1 -n1>=j2) | (j2 -n2>=j1)) return out;
349     // if both splines are first-order, the convolution is trivial
350     if((n1==1) & (n2==1)){
351         out= 1.0/(k+1.0)*(pow(knots[j1+splord],k+1.0)-pow(knots[j1-1+splord],k+1.0));
352         return out;
353     }
354     // assume that n1>n2 wlog, if this is not true, switch the indices
355     if(n2>n1){
356         int n3=n1; n1=n2; n2=n3;
357         int j3=j1; j1=j2; j2=j3;
358     }
359     // use a magic recursive formula!
360     out=0;
361     double denom1 = knots[j1-1+splord]-knots[j1-n1+splord];
362     double denom2 = knots[j1+splord] - knots[j1-n1+1+splord];
363     if(denom1>0){
364         out += 1.0/denom1 * cSplineConvolution(k+1,n1-1,j1-1,n2,j2,knots,splord);
365         out -= knots[j1-n1+splord]/denom1 * cSplineConvolution(k,n1-1,j1-1,n2,j2,knots,splord);
366     }
367     if(denom2>0){
368         out += knots[j1+splord]/denom2* cSplineConvolution(k,n1-1,j1,n2,j2,knots,splord);
369         out -= 1.0/denom2 * cSplineConvolution(k+1,n1-1,j1,n2,j2,knots, splord);
370     }
371     return out;
372 
373 }