functions.cpp

Go to the documentation of this file.

/* fuctions.C
   Collection of special functions

   NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)

   Revision history
   Date          Description                                     Programmer
   ------------  --------------------------------------------    ------------
   30/04/1992    Copied from "numerical recipies"                Rene Bakker

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "Algorithms/bet/math/functions.h"

/* Bessel function
   ======================================================================= */

static double bessj0(double x) {
    double ax, z;
    double xx, y, ans, ans1, ans2;

    if((ax = fabs(x)) < 8.0) {
        y = x * x;
        ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
                                    + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
        ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
                                    + y * (59272.64853 + y * (267.8532712 + y * 1.0))));
        ans = ans1 / ans2;
    } else {
        z = 8.0 / ax;
        y = z * z;
        xx = ax - 0.785398164;
        ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                          + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
        ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                                       + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                                               - y * 0.934935152e-7)));
        ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
    }
    return (double) ans;
} /* bessj0 */

static double bessj1(double x) {
    double ax, z;
    double xx, y, ans, ans1, ans2;

    if((ax = fabs(x)) < 8.0) {
        y = x * x;
        ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
                                         + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
        ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
                                     + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
        ans = ans1 / ans2;
    } else {
        z = 8.0 / ax;
        y = z * z;
        xx = ax - 2.356194491;
        ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
                                             + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
        ans2 = 0.04687499995 + y * (-0.2002690873e-3
                                    + y * (0.8449199096e-5 + y * (-0.88228987e-6
                                            + y * 0.105787412e-6)));
        ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
        if(x < 0.0) ans = -ans;
    }
    return (double) ans;
} /* bessj1 */

/*************************\
*                         *
*   External functions    *
*                         *
\*************************/

#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

double bessj(int n, double x) {
    int j, jsum, m;
    double ax, bj, bjm, bjp, sum, tox, ans;

    if(n < 0) {
        fprintf(stderr, "Illegal argument bessel function.\n");
        exit(1);
    }
    if(n == 0) return bessj0(x);
    if(n == 1) return bessj1(x);

    ax = fabs(x);
    if(ax == 0.0) return 0.0;
    else if(ax > (double) n) {
        tox = 2.0 / ax;
        bjm = bessj0(ax);
        bj = bessj1(ax);
        for(j = 1; j < n; j++) {
            bjp = j * tox * bj - bjm;
            bjm = bj;
            bj = bjp;
        }
        ans = bj;
    } else {
        tox = 2.0 / ax;
        m = 2 * ((n + (int) sqrt(ACC * n)) / 2);
        jsum = 0;
        bjp = ans = sum = 0.0;
        bj = 1.0;
        for(j = m; j > 0; j--) {
            bjm = j * tox * bj - bjp;
            bjp = bj;
            bj = bjm;
            if(fabs(bj) > BIGNO) {
                bj *= BIGNI;
                bjp *= BIGNI;
                ans *= BIGNI;
                sum *= BIGNI;
            }
            if(jsum) sum += bj;
            jsum = !jsum;
            if(j == n) ans = bjp;
        }
        sum = 2.0 * sum - bj;
        ans /= sum;
    }
    return x < 0.0 && n % 2 == 1 ? -ans : ans;
} /* bessj */



#undef ACC
#undef BIGNO
#undef BIGNI


/* Error function
   ======================================================================= */

/*
   errfc()
   Returns the complementary error function erfc(x) with fractional error
   everywhere less than 1.2e-07. */
double errfc(double x) {
    double t, z, ans;

    z = fabs(x);
    t = 1.0 / (1.0 + 0.5 * z);
    ans = t * exp(-z * z - 1.26551223 +
                  t * (1.00002368 +
                       t * (0.37409196 +
                            t * (0.09678418 +
                                 t * (-0.18628806 +
                                      t * (0.27886807 +
                                           t * (-1.13520398 +
                                                   t * (1.48851587 +
                                                           t * (-0.82215223 +
                                                                   t * 0.17087277)))))))));
    return x >= 0.0 ? ans : 2.0 - ans;
} /*erfc() */