cbessjy.cpp

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

//  cbessjy.cpp -- complex Bessel functions.
//  Algorithms and coefficient values from "Computation of Special
//  Functions", Zhang and Jin, John Wiley and Sons, 1996.
//
//  (C) 2003, C. Bond. All rights reserved.
//
#include <complex.h>
#include "bessel.h"
double gamma(double);

static complex<double> cii(0.0,1.0);
static complex<double> cone(1.0,0.0);
static complex<double> czero(0.0,0.0);

int cbessjy01(complex<double> z,complex<double> &cj0,complex<double> &cj1,
    complex<double> &cy0,complex<double> &cy1,complex<double> &cj0p,
    complex<double> &cj1p,complex<double> &cy0p,complex<double> &cy1p)
{
    complex<double> z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
    double a0,w0,w1;
    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};

    a0 = abs(z);
    z2 = z*z;
    z1 = z;
    if (a0 == 0.0) {
        cj0 = cone;
        cj1 = czero;
        cy0 = complex<double>(-1e308,0);
        cy1 = complex<double>(-1e308,0);
        cj0p = czero;
        cj1p = complex<double>(0.5,0.0);
        cy0p = complex<double>(1e308,0);
        cy1p = complex<double>(1e308,0);
        return 0;
    }
    if (real(z) < 0.0) z1 = -z;
    if (a0 <= 12.0) {
        cj0 = cone;
        cr = cone;
        for (k=1;k<=40;k++) {
            cr *= -0.25*z2/(double)(k*k);
            cj0 += cr;
            if (abs(cr) < abs(cj0)*eps) break;
        }
        cj1 = cone;
        cr = cone;
        for (k=1;k<=40;k++) {
            cr *= -0.25*z2/(k*(k+1.0));
            cj1 += cr;
            if (abs(cr) < abs(cj1)*eps) break;
        }
        cj1 *= 0.5*z1;
        w0 = 0.0;
        cr = cone;
        cs = czero;
        for (k=1;k<=40;k++) {
            w0 += 1.0/k;
            cr *= -0.25*z2/(double)(k*k);
            cp = cr*w0;
            cs += cp;
            if (abs(cp) < abs(cs)*eps) break;
        }
        cy0 = M_2_PI*((log(0.5*z1)+el)*cj0-cs);
        w1 = 0.0;
        cr = cone;
        cs = cone;
        for (k=1;k<=40;k++) {
            w1 += 1.0/k;
            cr *= -0.25*z2/(k*(k+1.0));
            cp = cr*(2.0*w1+1.0/(k+1.0));
            cs += cp;
            if (abs(cp) < abs(cs)*eps) break;
        }
        cy1 = M_2_PI*((log(0.5*z1)+el)*cj1-1.0/z1-0.25*z1*cs);
    }
    else {
        if (a0 >= 50.0) kz = 8;         // can be changed to 10
        else if (a0 >= 35.0) kz = 10;   //   "      "     "  12
        else kz = 12;                   //   "      "     "  14
        ct1 = z1 - M_PI_4;
        cp0 = cone;
        for (k=0;k<kz;k++) {
            cp0 += a[k]*pow(z1,-2.0*k-2.0);
        }
        cq0 = -0.125/z1;
        for (k=0;k<kz;k++) {
            cq0 += b[k]*pow(z1,-2.0*k-3.0);
        }
        cu = sqrt(M_2_PI/z1);
        cj0 = cu*(cp0*cos(ct1)-cq0*sin(ct1));
        cy0 = cu*(cp0*sin(ct1)+cq0*cos(ct1));
        ct2 = z1 - 0.75*M_PI;
        cp1 = cone;
        for (k=0;k<kz;k++) {
            cp1 += a1[k]*pow(z1,-2.0*k-2.0);
        }
        cq1 = 0.375/z1;
        for (k=0;k<kz;k++) {
            cq1 += b1[k]*pow(z1,-2.0*k-3.0);
        }
        cj1 = cu*(cp1*cos(ct2)-cq1*sin(ct2));
        cy1 = cu*(cp1*sin(ct2)+cq1*cos(ct2));
    }
    if (real(z) < 0.0) {
        if (imag(z) < 0.0) {
            cy0 -= 2.0*cii*cj0;
            cy1 = -(cy1-2.0*cii*cj1);
        }
        else if (imag(z) > 0.0) {
            cy0 += 2.0*cii*cj0;
            cy1 = -(cy1+2.0*cii*cj1);
        }
        cj1 = -cj1;
    }
    cj0p = -cj1;
    cj1p = cj0-cj1/z;
    cy0p = -cy1;
    cy1p = cy0-cy1/z;
    return 0;
}

int cbessjyna(int n,complex<double> z,int &nm,complex<double> *cj,
    complex<double> *cy,complex<double> *cjp,complex<double> *cyp)
{
    complex<double> cbj0,cbj1,cby0,cby1,cj0,cjk,cj1,cf,cf1,cf2;
    complex<double> cs,cg0,cg1,cyk,cyl1,cyl2,cylk,cp11,cp12,cp21,cp22;
    complex<double> ch0,ch1,ch2;
    double a0,yak,ya1,ya0,wa;
    int m,k,kz,lb,lb0;

    if (n < 0) return 1;
    a0 = abs(z);
    nm = n;
    if (a0 < 1.0e-100) {
        for (k=0;k<=n;k++) {
            cj[k] = czero;
            cy[k] = complex<double> (-1e308,0);
            cjp[k] = czero;
            cyp[k] = complex<double>(1e308,0);
        }
        cj[0] = cone;
        cjp[1] = complex<double>(0.5,0.0);
        return 0;
    }
    cbessjy01(z,cj[0],cj[1],cy[0],cy[1],cjp[0],cjp[1],cyp[0],cyp[1]);
    cbj0 = cj[0];
    cbj1 = cj[1];
    cby0 = cy[0];
    cby1 = cy[1];
    if (n <= 1) return 0;
    if (n < (int)0.25*a0) {
        cj0 = cbj0;
        cj1 = cbj1;
        for (k=2;k<=n;k++) {
            cjk = 2.0*(k-1.0)*cj1/z-cj0;
            cj[k] = cjk;
            cj0 = cj1;
            cj1 = cjk;
        }
    }
    else {
        m = msta1(a0,200);
        if (m < n) nm = m;
        else m = msta2(a0,n,15);
        cf2 = czero;
        cf1 = complex<double> (1.0e-100,0.0);
        for (k=m;k>=0;k--) {
            cf = 2.0*(k+1.0)*cf1/z-cf2;
            if (k <=nm) cj[k] = cf;
            cf2 = cf1;
            cf1 = cf;
        }
        if (abs(cbj0) > abs(cbj1)) cs = cbj0/cf;
        else cs = cbj1/cf2;
        for (k=0;k<=nm;k++) {
            cj[k] *= cs;
        }
    }
    for (k=2;k<=nm;k++) {
        cjp[k] = cj[k-1]-(double)k*cj[k]/z;
    }
    ya0 = abs(cby0);
    lb = 0;
    cg0 = cby0;
    cg1 = cby1;
    for (k=2;k<=nm;k++) {
        cyk = 2.0*(k-1.0)*cg1/z-cg0;
        yak = abs(cyk);
        ya1 = abs(cg0);
        if ((yak < ya0) && (yak < ya1)) lb = k;
        cy[k] = cyk;
        cg0 = cg1;
        cg1 = cyk;
    }
    lb0 = 0;
    if ((lb > 4) && (imag(z) != 0.0)) {
        while (lb != lb0) {
            ch2 = cone;
            ch1 = czero;
            lb0 = lb;
            for (k=lb;k>=1;k--) {
                ch0 = 2.0*k*ch1/z-ch2;
                ch2 = ch1;
                ch1 = ch0;
            }
            cp12 = ch0;
            cp22 = ch2;
            ch2 = czero;
            ch1 = cone;
            for (k=lb;k>=1;k--) {
                ch0 = 2.0*k*ch1/z-ch2;
                ch2 = ch1;
                ch1 = ch0;
            }
            cp11 = ch0;
            cp21 = ch2;
            if (lb == nm)
                cj[lb+1] = 2.0*lb*cj[lb]/z-cj[lb-1];
            if (abs(cj[0]) > abs(cj[1])) {
                cy[lb+1] = (cj[lb+1]*cby0-2.0*cp11/(M_PI*z))/cj[0];
                cy[lb] = (cj[lb]*cby0+2.0*cp12/(M_PI*z))/cj[0];
            }
            else {
                cy[lb+1] = (cj[lb+1]*cby1-2.0*cp21/(M_PI*z))/cj[1];
                cy[lb] = (cj[lb]*cby1+2.0*cp22/(M_PI*z))/cj[1];
            }
            cyl2 = cy[lb+1];
            cyl1 = cy[lb];
            for (k=lb-1;k>=0;k--) {
                cylk = 2.0*(k+1.0)*cyl1/z-cyl2;
                cy[k] = cylk;
                cyl2 = cyl1;
                cyl1 = cylk;
            }
            cyl1 = cy[lb];
            cyl2 = cy[lb+1];
            for (k=lb+1;k<n;k++) {
                cylk = 2.0*k*cyl2/z-cyl1;
                cy[k+1] = cylk;
                cyl1 = cyl2;
                cyl2 = cylk;
            }
            for (k=2;k<=nm;k++) {
                wa = abs(cy[k]);
                if (wa < abs(cy[k-1])) lb = k;
            }
        }
    }
    for (k=2;k<=nm;k++) {
        cyp[k] = cy[k-1]-(double)k*cy[k]/z;
    }
    return 0;
}

int cbessjynb(int n,complex<double> z,int &nm,complex<double> *cj,
    complex<double> *cy,complex<double> *cjp,complex<double> *cyp)
{
    complex<double> cf,cf0,cf1,cf2,cbs,csu,csv,cs0,ce;
    complex<double> ct1,cp0,cq0,cp1,cq1,cu,cbj0,cby0,cbj1,cby1;
    complex<double> cyy,cbjk,ct2;
    double a0,y0;
    int k,m;
    static double a[] = {
        -0.7031250000000000e-1,
         0.1121520996093750,
        -0.5725014209747314,
         6.074042001273483};
    static double b[] = {
         0.7324218750000000e-1,
        -0.2271080017089844,
         1.727727502584457,
        -2.438052969955606e1};
    static double a1[] = {
         0.1171875,
        -0.1441955566406250,
         0.6765925884246826,
        -6.883914268109947};
    static double b1[] = {
       -0.1025390625,
        0.2775764465332031,
       -1.993531733751297,
        2.724882731126854e1};

    y0 = abs(imag(z));
    a0 = abs(z);
    nm = n;
    if (a0 < 1.0e-100) {
        for (k=0;k<=n;k++) {
            cj[k] = czero;
            cy[k] = complex<double> (-1e308,0);
            cjp[k] = czero;
            cyp[k] = complex<double>(1e308,0);
        }
        cj[0] = cone;
        cjp[1] = complex<double>(0.5,0.0);
        return 0;