bessjy.cpp

Individual files are available in the following links: Bessjy -- Bessel functions Jn and Yn for r

//  bessjy.cpp -- computation of Bessel functions Jn, Yn and their
//      derivatives. Algorithms and coefficient values from
//      "Computation of Special Functions", Zhang and Jin, John
//      Wiley and Sons, 1996.
//
//  (C) 2003, C. Bond. All rights reserved.
//
// Note that 'math.h' provides (or should provide) values for:
//      pi      M_PI
//      2/pi    M_2_PI
//      pi/4    M_PI_4
//      pi/2    M_PI_2
//
#include <math.h>
#include "bessel.h"

double gamma(double x);
//
//  INPUT:
//      double x    -- argument of Bessel function
//
//  OUTPUT (via address pointers):
//      double j0   -- Bessel function of 1st kind, 0th order
//      double j1   -- Bessel function of 1st kind, 1st order
//      double y0   -- Bessel function of 2nd kind, 0th order
//      double y1   -- Bessel function of 2nd kind, 1st order
//      double j0p  -- derivative of Bessel function of 1st kind, 0th order
//      double j1p  -- derivative of Bessel function of 1st kind, 1st order
//      double y0p  -- derivative of Bessel function of 2nd kind, 0th order
//      double y1p  -- derivative of Bessel function of 2nd kind, 1st order
//
//  RETURN:
//      int error code: 0 = OK, 1 = error
//
//  This algorithm computes the above functions using series expansions.
//
int bessjy01a(double x,double &j0,double &j1,double &y0,double &y1,
    double &j0p,double &j1p,double &y0p,double &y1p)
{
    double x2,r,ec,w0,w1,r0,r1,cs0,cs1;
    double cu,p0,q0,p1,q1,t1,t2;
    int k,kz;
    static double a[] = {
        -7.03125e-2,
         0.112152099609375,
        -0.5725014209747314,
         6.074042001273483,
        -1.100171402692467e2,
         3.038090510922384e3,
        -1.188384262567832e5,
         6.252951493434797e6,
        -4.259392165047669e8,
         3.646840080706556e10,
        -3.833534661393944e12,
         4.854014686852901e14,
        -7.286857349377656e16,
         1.279721941975975e19};
    static double b[] = {
         7.32421875e-2,
        -0.2271080017089844,
         1.727727502584457,
        -2.438052969955606e1,
         5.513358961220206e2,
        -1.825775547429318e4,
         8.328593040162893e5,
        -5.006958953198893e7,
         3.836255180230433e9,
        -3.649010818849833e11,
         4.218971570284096e13,
        -5.827244631566907e15,
         9.476288099260110e17,
        -1.792162323051699e20};
    static double a1[] = {
         0.1171875,
        -0.1441955566406250,
         0.6765925884246826,
        -6.883914268109947,
         1.215978918765359e2,
        -3.302272294480852e3,
         1.276412726461746e5,
        -6.656367718817688e6,
         4.502786003050393e8,
        -3.833857520742790e10,
         4.011838599133198e12,
        -5.060568503314727e14,
         7.572616461117958e16,
        -1.326257285320556e19};
    static double b1[] = {
        -0.1025390625,
         0.2775764465332031,
        -1.993531733751297,
         2.724882731126854e1,
        -6.038440767050702e2,
         1.971837591223663e4,
        -8.902978767070678e5,
         5.310411010968522e7,
        -4.043620325107754e9,
         3.827011346598605e11,
        -4.406481417852278e13,
         6.065091351222699e15,
        -9.833883876590679e17,
         1.855045211579828e20};

    if (x < 0.0) return 1;
    if (x == 0.0) {
        j0 = 1.0;
        j1 = 0.0;
        y0 = -1e308;
        y1 = -1e308;
        j0p = 0.0;
        j1p = 0.5;
        y0p = 1e308;
        y1p = 1e308;
        return 0;
    }
    x2 = x*x;
    if (x <= 12.0) {
        j0 = 1.0;
        r = 1.0;
        for (k=1;k<=30;k++) {
            r *= -0.25*x2/(k*k);
            j0 += r;
            if (fabs(r) < fabs(j0)*1e-15) break;
        }
        j1 = 1.0;
        r = 1.0;
        for (k=1;k<=30;k++) {
            r *= -0.25*x2/(k*(k+1));
            j1 += r;
            if (fabs(r) < fabs(j1)*1e-15) break;
        }
        j1 *= 0.5*x;
        ec = log(0.5*x)+el;
        cs0 = 0.0;
        w0 = 0.0;
        r0 = 1.0;
        for (k=1;k<=30;k++) {
            w0 += 1.0/k;
            r0 *= -0.25*x2/(k*k);
            r = r0 * w0;
            cs0 += r;
            if (fabs(r) < fabs(cs0)*1e-15) break;
        }
        y0 = M_2_PI*(ec*j0-cs0);
        cs1 = 1.0;
        w1 = 0.0;
        r1 = 1.0;
        for (k=1;k<=30;k++) {
            w1 += 1.0/k;
            r1 *= -0.25*x2/(k*(k+1));
            r = r1*(2.0*w1+1.0/(k+1));
            cs1 += r;
            if (fabs(r) < fabs(cs1)*1e-15) break;
        }
        y1 = M_2_PI * (ec*j1-1.0/x-0.25*x*cs1);
    }
    else {
        if (x >= 50.0) kz = 8;          // Can be changed to 10
        else if (x >= 35.0) kz = 10;    //  "       "        12
        else kz = 12;                   //  "       "        14
        t1 = x-M_PI_4;
        p0 = 1.0;
        q0 = -0.125/x;
        for (k=0;k<kz;k++) {
            p0 += a[k]*pow(x,-2*k-2);
            q0 += b[k]*pow(x,-2*k-3);
        }
        cu = sqrt(M_2_PI/x);
        j0 = cu*(p0*cos(t1)-q0*sin(t1));
        y0 = cu*(p0*sin(t1)+q0*cos(t1));
        t2 = x-0.75*M_PI;
        p1 = 1.0;
        q1 = 0.375/x;
        for (k=0;k<kz;k++) {
            p1 += a1[k]*pow(x,-2*k-2);
            q1 += b1[k]*pow(x,-2*k-3);
        }
        j1 = cu*(p1*cos(t2)-q1*sin(t2));
        y1 = cu*(p1*sin(t2)+q1*cos(t2));
    }
    j0p = -j1;
    j1p = j0-j1/x;
    y0p = -y1;
    y1p = y0-y1/x;
    return 0;
}
//
//  INPUT:
//      double x    -- argument of Bessel function
//
//  OUTPUT:
//      double j0   -- Bessel function of 1st kind, 0th order
//      double j1   -- Bessel function of 1st kind, 1st order
//      double y0   -- Bessel function of 2nd kind, 0th order
//      double y1   -- Bessel function of 2nd kind, 1st order
//      double j0p  -- derivative of Bessel function of 1st kind, 0th order
//      double j1p  -- derivative of Bessel function of 1st kind, 1st order
//      double y0p  -- derivative of Bessel function of 2nd kind, 0th order
//      double y1p  -- derivative of Bessel function of 2nd kind, 1st order
//
//  RETURN:
//      int error code: 0 = OK, 1 = error
//
//  This algorithm computes the functions using polynomial approximations.
//
int bessjy01b(double x,double &j0,double &j1,double &y0,double &y1,
    double &j0p,double &j1p,double &y0p,double &y1p)
{
    double t,t2,dtmp,a0,p0,q0,p1,q1,ta0,ta1;
    if (x < 0.0) return 1;
    if (x == 0.0) {
        j0 = 1.0;
        j1 = 0.0;
        y0 = -1e308;
        y1 = -1e308;
        j0p = 0.0;
        j1p = 0.5;
        y0p = 1e308;
        y1p = 1e308;
        return 0;
    }
    if(x <= 4.0) {
        t = x/4.0;
        t2 = t*t;
        j0 = ((((((-0.5014415e-3*t2+0.76771853e-2)*t2-0.0709253492)*t2+
            0.4443584263)*t2-1.7777560599)*t2+3.9999973021)*t2
            -3.9999998721)*t2+1.0;
        j1 = t*(((((((-0.1289769e-3*t2+0.22069155e-2)*t2-0.0236616773)*t2+
            0.1777582922)*t2-0.8888839649)*t2+2.6666660544)*t2-
            3.999999971)*t2+1.9999999998);
        dtmp = (((((((-0.567433e-4*t2+0.859977e-3)*t2-0.94855882e-2)*t2+
            0.0772975809)*t2-0.4261737419)*t2+1.4216421221)*t2-
            2.3498519931)*t2+1.0766115157)*t2+0.3674669052;
        y0 = M_2_PI*log(0.5*x)*j0+dtmp;
        dtmp = (((((((0.6535773e-3*t2-0.0108175626)*t2+0.107657607)*t2-
            0.7268945577)*t2+3.1261399273)*t2-7.3980241381)*t2+
            6.8529236342)*t2+0.3932562018)*t2-0.6366197726;
        y1 = M_2_PI*log(0.5*x)*j1+dtmp/x;
    }
    else {
        t = 4.0/x;
        t2 = t*t;
        a0 = sqrt(M_2_PI/x);
        p0 = ((((-0.9285e-5*t2+0.43506e-4)*t2-0.122226e-3)*t2+
             0.434725e-3)*t2-0.4394275e-2)*t2+0.999999997;
        q0 = t*(((((0.8099e-5*t2-0.35614e-4)*t2+0.85844e-4)*t2-
            0.218024e-3)*t2+0.1144106e-2)*t2-0.031249995);
        ta0 = x-M_PI_4;
        j0 = a0*(p0*cos(ta0)-q0*sin(ta0));
        y0 = a0*(p0*sin(ta0)+q0*cos(ta0));
        p1 = ((((0.10632e-4*t2-0.50363e-4)*t2+0.145575e-3)*t2
            -0.559487e-3)*t2+0.7323931e-2)*t2+1.000000004;
        q1 = t*(((((-0.9173e-5*t2+0.40658e-4)*t2-0.99941e-4)*t2
            +0.266891e-3)*t2-0.1601836e-2)*t2+0.093749994);
        ta1 = x-0.75*M_PI;
        j1 = a0*(p1*cos(ta1)-q1*sin(ta1));
        y1 = a0*(p1*sin(ta1)+q1*cos(ta1));
    }
    j0p = -j1;
    j1p = j0-j1/x;
    y0p = -y1;
    y1p = y0-y1/x;
    return 0;
}
int msta1(double x,int mp)
{
    double a0,f0,f1,f;
    int i,n0,n1,nn;

    a0 = fabs(x);
    n0 = (int)(1.1*a0)+1;
    f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-mp;
    n1 = n0+5;
    f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp;
    for (i=0;i<20;i++) {
        nn = n1-(n1-n0)/(1.0-f0/f1);
        f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp;
        if (abs(nn-n1) < 1) break;
        n0 = n1;
        f0 = f1;
        n1 = nn;
        f1 = f;
    }
    return nn;
}
int msta2(double x,int n,int mp)
{
    double a0,ejn,hmp,f0,f1,f,obj;
    int i,n0,n1,nn;

    a0 = fabs(x);
    hmp = 0.5*mp;
    ejn = 0.5*log10(6.28*n)-n*log10(1.36*a0/n);
    if (ejn <= hmp) {
        obj = mp;
        n0 = (int)(1.1*a0);
        if (n0 < 1) n0 = 1;
    }
    else {
        obj = hmp+ejn;
        n0 = n;
    }
    f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-obj;
    n1 = n0+5;
    f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj;
    for (i=0;i<20;i++) {
        nn = n1-(n1-n0)/(1.0-f0/f1);
        f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj;
        if (abs(nn-n1) < 1) break;
        n0 = n1;
        f0 = f1;
        n1 = nn;
        f1 = f;
    }
    return nn+10;
}
//
//  INPUT:
//  double x    -- argument of Bessel function of 1st and 2nd kind.
//  int n       -- order
//
//  OUPUT:
//
//  int nm      -- highest order actually computed (nm <= n)
//  double jn[] -- Bessel function of 1st kind, orders from 0 to nm
//  double yn[] -- Bessel function of 2nd kind, orders from 0 to nm
//  double j'n[]-- derivative of Bessel function of 1st kind,
//                      orders from 0 to nm
//  double y'n[]-- derivative of Bessel function of 2nd kind,
//                      orders from 0 to nm
//
//  Computes Bessel functions of all order up to 'n' using recurrence
//  relations. If 'nm' < 'n' only 'nm' orders are returned.
//
int bessjyna(int n,double x,int &nm,double *jn,double *yn,
    double *jnp,double *ynp)
{
    double bj0,bj1,f,f0,f1,f2,cs;
    int i,k,m,ecode;

    nm = n;
    if ((x < 0.0) || (n < 0)) return 1;
    if (x < 1e-15) {
        for (i=0;i<=n;i++) {
            jn[i] = 0.0;
            yn[i] = -1e308;
            jnp[i] = 0.0;
            ynp[i] = 1e308;
        }
        jn[0] = 1.0;
        jnp[1] = 0.5;
        return 0;
    }