cbessjy.cpp

Go to the documentation of this file.

//  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 <lifev/navier_stokes/function/bessel/bessel.hpp>

namespace bessel
{
double gamma (double);

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

int cbessjy01 (std::complex<double> z, std::complex<double>& cj0, std::complex<double>& cj1,
               std::complex<double>& cy0, std::complex<double>& cy1, std::complex<double>& cj0p,
               std::complex<double>& cj1p, std::complex<double>& cy0p, std::complex<double>& cy1p)
{
    std::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 = std::abs (z);
    z2 = z * z;
    z1 = z;
    if (a0 == 0.0)
    {
        cj0 = cone;
        cj1 = czero;
        cy0 = std::complex<double> (-1e308, 0);
        cy1 = std::complex<double> (-1e308, 0);
        cj0p = czero;
        cj1p = std::complex<double> (0.5, 0.0);
        cy0p = std::complex<double> (1e308, 0);
        cy1p = std::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 (std::abs (cr) < std::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 (std::abs (cr) < std::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 (std::abs (cp) < std::abs (cs) *eps)
            {
                break;
            }
        }
        cy0 = M_2_PI * ( (std::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 (std::abs (cp) < std::abs (cs) *eps)
            {
                break;
            }
        }
        cy1 = M_2_PI * ( (std::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] * std::pow (z1, -2.0 * k - 2.0);
        }
        cq0 = -0.125 / z1;
        for (k = 0; k < kz; k++)
        {
            cq0 += b[k] * std::pow (z1, -2.0 * k - 3.0);
        }
        cu = std::sqrt (M_2_PI / z1);
        cj0 = cu * (cp0 * std::cos (ct1) - cq0 * std::sin (ct1) );
        cy0 = cu * (cp0 * std::sin (ct1) + cq0 * std::cos (ct1) );
        ct2 = z1 - 0.75 * M_PI;
        cp1 = cone;
        for (k = 0; k < kz; k++)
        {
            cp1 += a1[k] * std::pow (z1, -2.0 * k - 2.0);
        }
        cq1 = 0.375 / z1;
        for (k = 0; k < kz; k++)
        {
            cq1 += b1[k] * std::pow (z1, -2.0 * k - 3.0);
        }
        cj1 = cu * (cp1 * std::cos (ct2) - cq1 * std::sin (ct2) );
        cy1 = cu * (cp1 * std::sin (ct2) + cq1 * std::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, std::complex<double> z, int& nm, std::complex<double>* cj,
               std::complex<double>* cy, std::complex<double>* cjp, std::complex<double>* cyp)
{
    std::complex<double> cbj0, cbj1, cby0, cby1, cj0, cjk, cj1, cf, cf1, cf2;
    std::complex<double> cs, cg0, cg1, cyk, cyl1, cyl2, cylk, cp11, cp12, cp21, cp22;
    std::complex<double> ch0, ch1, ch2;
    double a0, yak, ya1, ya0, wa;
    int m, k, lb, lb0;

    if (n < 0)
    {
        return 1;
    }
    a0 = std::abs (z);
    nm = n;
    if (a0 < 1.0e-100)
    {
        for (k = 0; k <= n; k++)
        {
            cj[k] = czero;
            cy[k] = std::complex<double> (-1e308, 0);
            cjp[k] = czero;
            cyp[k] = std::complex<double> (1e308, 0);
        }
        cj[0] = cone;
        cjp[1] = std::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 = std::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 (std::abs (cbj0) > std::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 = std::abs (cby0);
    lb = 0;
    cg0 = cby0;
    cg1 = cby1;
    for (k = 2; k <= nm; k++)
    {
        cyk = 2.0 * (k - 1.0) * cg1 / z - cg0;
        yak = std::abs (cyk);
        ya1 = std::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 (std::abs (cj[0]) > std::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 = std::abs (cy[k]);
                if (wa < std::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, std::complex<double> z, int& nm, std::complex<double>* cj,
               std::complex<double>* cy, std::complex<double>* cjp, std::complex<double>* cyp)
{
    std::complex<double> cf, cf0, cf1, cf2, cbs, csu, csv, cs0, ce;
    std::complex<double> ct1, cp0, cq0, cp1, cq1, cu, cbj0, cby0, cbj1, cby1;
    std::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 = std::abs (imag (z) );
    a0 = std::abs (z);
    nm = n;
    if (a0 < 1.0e-100)
    {
        for (k = 0; k <= n; k++)
        {
            cj[k] = czero;
            cy[k] = std::complex<double> (-1e308, 0);
            cjp[k] = czero;
            cyp[k] = std::complex<double> (1e308, 0);
        }
        cj[0] = cone;
        cjp[1] = std::complex<double> (0.5, 0.0);
        return 0;
    }
    if ( (a0 <= 300.0) || (n > (int) (0.25 * a0) ) )
    {
        if (n == 0)
        {
            nm = 1;
        }
        m = msta1 (a0, 200);
        if (m < nm)
        {
            nm = m;
        }
        else
        {
            m = msta2 (a0, nm, 15);
        }
        cbs = czero;
        csu = czero;
        csv = czero;
        cf2 = czero;
        cf1 = std::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;
            }
            if ( ( (k & 1) == 0) && (k != 0) )
            {
                if (y0 <= 1.0)
                {
                    cbs += 2.0 * cf;
                }
                else
                {
                    cbs += (-1) * ( (k & 2) - 1) * 2.0 * cf;
                }
                csu += (double) ( (-1) * ( (k & 2) - 1) ) * cf / (double) k;
            }
            else if (k > 1)
            {
                csv += (double) ( (-1) * ( (k & 2) - 1) * k) * cf / (double) (k * k - 1.0);
            }
            cf2 = cf1;
            cf1 = cf;
        }
        if (y0 <= 1.0)
        {
            cs0 = cbs + cf;
        }
        else
        {
            cs0 = (cbs + cf) / std::cos (z);
        }
        for (k = 0; k <= nm; k++)
        {
            cj[k] /= cs0;
        }
        ce = std::log (0.5 * z) + el;
        cy[0] = M_2_PI * (ce * cj[0] - 4.0 * csu / cs0);
        cy[1] = M_2_PI * (-cj[0] / z + (ce - 1.0) * cj[1] - 4.0 * csv / cs0);
    }
    else
    {
        ct1 = z - M_PI_4;
        cp0 = cone;
        for (k = 0; k < 4; k++)
        {
            cp0 += a[k] * std::pow (z, -2.0 * k - 2.0);
        }
        cq0 = -0.125 / z;
        for (k = 0; k < 4; k++)
        {
            cq0 += b[k] * std::pow (z, -2.0 * k - 3.0);
        }
        cu = std::sqrt (M_2_PI / z);
        cbj0 = cu * (cp0 * std::cos (ct1) - cq0 * std::sin (ct1) );
        cby0 = cu * (cp0 * std::sin (ct1) + cq0 * std::cos (ct1) );
        cj[0] = cbj0;
        cy[0] = cby0;
        ct2 = z - 0.75 * M_PI;
        cp1 = cone;
        for (k = 0; k < 4; k++)
        {
            cp1 += a1[k] * std::pow (z, -2.0 * k - 2.0);
        }
        cq1 = 0.375 / z;
        for (k = 0; k < 4; k++)
        {
            cq1 += b1[k] * std::pow (z, -2.0 * k - 3.0);
        }
        cbj1 = cu * (cp1 * std::cos (ct2) - cq1 * std::sin (ct2) );
        cby1 = cu * (cp1 * std::sin (ct2) + cq1 * std::cos (ct2) );
        cj[1] = cbj1;
        cy[1] = cby1;
        for (k = 2; k <= n; k++)
        {
            cbjk = 2.0 * (k - 1.0) * cbj1 / z - cbj0;
            cj[k] = cbjk;
            cbj0 = cbj1;
            cbj1 = cbjk;
        }
    }
    cjp[0] = -cj[1];
    for (k = 1; k <= nm; k++)
    {
        cjp[k] = cj[k - 1] - (double) k * cj[k] / z;
    }
    if (std::abs (cj[0]) > 1.0)
    {
        cy[1] = (cj[1] * cy[0] - 2.0 / (M_PI * z) ) / cj[0];
    }
    for (k = 2; k <= nm; k++)
    {
        if (std::abs (cj[k - 1]) >= std::abs (cj[k - 2]) )
        {
            cyy = (cj[k] * cy[k - 1] - 2.0 / (M_PI * z) ) / cj[k - 1];
        }
        else
        {
            cyy = (cj[k] * cy[k - 2] - 4.0 * (k - 1.0) / (M_PI * z * z) ) / cj[k - 2];
        }
        cy[k] = cyy;
    }
    cyp[0] = -cy[1];
    for (k = 1; k <= nm; k++)
    {
        cyp[k] = cy[k - 1] - (double) k * cy[k] / z;
    }

    return 0;
}

int cbessjyva (double v, std::complex<double> z, double& vm, std::complex<double>* cjv,
               std::complex<double>* cyv, std::complex<double>* cjvp, std::complex<double>* cyvp)
{
    std::complex<double> z1, z2, zk, cjvl, cr, ca, cjv0, cjv1, cpz, crp;
    std::complex<double> cqz, crq, ca0, cck, csk, cyv0, cyv1, cju0, cju1, cb;
    std::complex<double> cs, cs0, cr0, cs1, cr1, cec, cf, cf0, cf1, cf2;
    std::complex<double> cfac0, cfac1, cg0, cg1, cyk, cp11, cp12, cp21, cp22;
    std::complex<double> ch0, ch1, ch2, cyl1, cyl2, cylk;

    double a0, v0, pv0, pv1, vl, ga, gb, vg, vv, w0, w1, ya0, yak, ya1, wa;
    int j, n, k, kz, l, lb, lb0, m;

    a0 = std::abs (z);
    z1 = z;
    z2 = z * z;
    n = (int) v;


    v0 = v - n;

    pv0 = M_PI * v0;
    pv1 = M_PI * (1.0 + v0);
    if (a0 < 1.0e-100)
    {
        for (k = 0; k <= n; k++)
        {
            cjv[k] = czero;
            cyv[k] = std::complex<double> (-1e308, 0);
            cjvp[k] = czero;
            cyvp[k] = std::complex<double> (1e308, 0);

        }
        if (v0 == 0.0)
        {
            cjv[0] = cone;
            cjvp[1] = std::complex<double> (0.5, 0.0);
        }
        else
        {
            cjvp[0] = std::complex<double> (1e308, 0);
        }
        vm = v;
        return 0;
    }
    if (real (z1) < 0.0)
    {
        z1 = -z;
    }
    if (a0 <= 12.0)
    {
        for (l = 0; l < 2; l++)
        {
            vl = v0 + l;
            cjvl = cone;
            cr = cone;
            for (k = 1; k <= 40; k++)
            {
                cr *= -0.25 * z2 / (k * (k + vl) );
                cjvl += cr;
                if (std::abs (cr) < std::abs (cjvl) *eps)
                {
                    break;
                }
            }
            vg = 1.0 + vl;
            ga = gamma (vg);
            ca = std::pow (0.5 * z1, vl) / ga;
            if (l == 0)
            {
                cjv0 = cjvl * ca;
            }
            else
            {
                cjv1 = cjvl * ca;
            }
        }
    }
    else
    {
        if (a0 >= 50.0)
        {
            kz = 8;
        }
        else if (a0 >= 35.0)
        {
            kz = 10;
        }
        else
        {
            kz = 11;
        }
        for (j = 0; j < 2; j++)
        {
            vv = 4.0 * (j + v0) * (j + v0);
            cpz = cone;
            crp = cone;
            for (k = 1; k <= kz; k++)
            {
                crp = -0.78125e-2 * crp * (vv - std::pow (4.0 * k - 3.0, 2.0) ) *
                      (vv - std::pow (4.0 * k - 1.0, 2.0) ) / (k * (2.0 * k - 1.0) * z2);
                cpz += crp;
            }
            cqz = cone;
            crq = cone;
            for (k = 1; k <= kz; k++)
            {
                crq = -0.78125e-2 * crq * (vv - std::pow (4.0 * k - 1.0, 2.0) ) *
                      (vv - std::pow (4.0 * k + 1.0, 2.0) ) / (k * (2.0 * k + 1.0) * z2);
                cqz += crq;
            }
            cqz *= 0.125 * (vv - 1.0) / z1;
            zk = z1 - (0.5 * (j + v0) + 0.25) * M_PI;
            ca0 = std::sqrt (M_2_PI / z1);
            cck = std::cos (zk);
            csk = std::sin (zk);
            if (j == 0)
            {
                cjv0 = ca0 * (cpz * cck - cqz * csk);
                cyv0 = ca0 * (cpz * csk + cqz + cck);
            }
            else
            {
                cjv1 = ca0 * (cpz * cck - cqz * csk);
                cyv1 = ca0 * (cpz * csk + cqz * cck);
            }
        }
    }
    if (a0 <= 12.0)
    {
        if (v0 != 0.0)
        {
            for (l = 0; l < 2; l++)
            {
                vl = v0 + l;
                cjvl = cone;
                cr = cone;
                for (k = 1; k <= 40; k++)
                {
                    cr *= -0.25 * z2 / (k * (k - vl) );
                    cjvl += cr;
                    if (std::abs (cr) < std::abs (cjvl) *eps)
                    {
                        break;
                    }
                }
                vg = 1.0 - vl;
                gb = gamma (vg);
                cb = std::pow (2.0 / z1, vl) / gb;
                if (l == 0)
                {
                    cju0 = cjvl * cb;
                }
                else
                {
                    cju1 = cjvl * cb;
                }
            }
            cyv0 = (cjv0 * std::cos (pv0) - cju0) / std::sin (pv0);
            cyv1 = (cjv1 * std::cos (pv1) - cju1) / std::sin (pv1);
        }
        else
        {
            cec = std::log (0.5 * z1) + el;
            cs0 = czero;
            w0 = 0.0;
            cr0 = cone;
            for (k = 1; k <= 30; k++)
            {
                w0 += 1.0 / k;
                cr0 *= -0.25 * z2 / (double) (k * k);
                cs0 += cr0 * w0;
            }
            cyv0 = M_2_PI * (cec * cjv0 - cs0);
            cs1 = cone;
            w1 = 0.0;
            cr1 = cone;
            for (k = 1; k <= 30; k++)
            {
                w1 += 1.0 / k;
                cr1 *= -0.25 * z2 / (k * (k + 1.0) );
                cs1 += cr1 * (2.0 * w1 + 1.0 / (k + 1.0) );
            }
            cyv1 = M_2_PI * (cec * cjv1 - 1.0 / z1 - 0.25 * z1 * cs1);
        }
    }
    if (real (z) < 0.0)
    {
        cfac0 = std::exp (pv0 * cii);
        cfac1 = std::exp (pv1 * cii);
        if (imag (z) < 0.0)
        {
            cyv0 = cfac0 * cyv0 - 2.0 * cii * std::cos (pv0) * cjv0;
            cyv1 = cfac1 * cyv1 - 2.0 * cii * std::cos (pv1) * cjv1;
            cjv0 /= cfac0;
            cjv1 /= cfac1;
        }
        else if (imag (z) > 0.0)
        {
            cyv0 = cyv0 / cfac0 + 2.0 * cii * std::cos (pv0) * cjv0;
            cyv1 = cyv1 / cfac1 + 2.0 * cii * std::cos (pv1) * cjv1;
            cjv0 *= cfac0;
            cjv1 *= cfac1;
        }
    }
    cjv[0] = cjv0;
    cjv[1] = cjv1;
    if ( (n >= 2) && (n <= (int) (0.25 * a0) ) )
    {
        cf0 = cjv0;
        cf1 = cjv1;
        for (k = 2; k <= n; k++)
        {
            cf = 2.0 * (k + v0 - 1.0) * cf1 / z - cf0;
            cjv[k] = cf;
            cf0 = cf1;
            cf1 = cf;
        }
    }
    else if (n >= 2)
    {
        m = msta1 (a0, 200);
        if (m < n)
        {
            n = m;
        }
        else
        {
            m = msta2 (a0, n, 15);
        }
        cf2 = czero;
        cf1 = std::complex<double> (1.0e-100, 0.0);
        for (k = m; k >= 0; k--)
        {
            cf = 2.0 * (v0 + k + 1.0) * cf1 / z - cf2;
            if (k <= n)
            {
                cjv[k] = cf;
            }
            cf2 = cf1;
            cf1 = cf;
        }
        if (std::abs (cjv0) > std::abs (cjv1) )
        {
            cs = cjv0 / cf;
        }
        else
        {
            cs = cjv1 / cf2;
        }
        for (k = 0; k <= n; k++)
        {
            cjv[k] *= cs;
        }
    }
    cjvp[0] = v0 * cjv[0] / z - cjv[1];
    for (k = 1; k <= n; k++)
    {
        cjvp[k] = - (k + v0) * cjv[k] / z + cjv[k - 1];
    }
    cyv[0] = cyv0;
    cyv[1] = cyv1;
    ya0 = std::abs (cyv0);
    lb = 0;
    cg0 = cyv0;
    cg1 = cyv1;
    for (k = 2; k <= n; k++)
    {
        cyk = 2.0 * (v0 + k - 1.0) * cg1 / z - cg0;
        yak = std::abs (cyk);
        ya1 = std::abs (cg0);
        if ( (yak < ya0) && (yak < ya1) )
        {
            lb = k;
        }
        cyv[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 + v0) * 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 + v0) * ch1 / z - ch2;
                ch2 = ch1;
                ch1 = ch0;
            }
            cp11 = ch0;
            cp21 = ch2;
            if (lb == n)
            {
                cjv[lb + 1] = 2.0 * (lb + v0) * cjv[lb] / z - cjv[lb - 1];
            }
            if (std::abs (cjv[0]) > std::abs (cjv[1]) )
            {
                cyv[lb + 1] = (cjv[lb + 1] * cyv0 - 2.0 * cp11 / (M_PI * z) ) / cjv[0];
                cyv[lb] = (cjv[lb] * cyv0 + 2.0 * cp12 / (M_PI * z) ) / cjv[0];
            }
            else
            {
                cyv[lb + 1] = (cjv[lb + 1] * cyv1 - 2.0 * cp21 / (M_PI * z) ) / cjv[1];
                cyv[lb] = (cjv[lb] * cyv1 + 2.0 * cp22 / (M_PI * z) ) / cjv[1];
            }
            cyl2 = cyv[lb + 1];
            cyl1 = cyv[lb];
            for (k = lb - 1; k >= 0; k--)
            {
                cylk = 2.0 * (k + v0 + 1.0) * cyl1 / z - cyl2;
                cyv[k] = cylk;
                cyl2 = cyl1;
                cyl1 = cylk;
            }
            cyl1 = cyv[lb];
            cyl2 = cyv[lb + 1];
            for (k = lb + 1; k < n; k++)
            {
                cylk = 2.0 * (k + v0) * cyl2 / z - cyl1;
                cyv[k + 1] = cylk;
                cyl1 = cyl2;
                cyl2 = cylk;
            }
            for (k = 2; k <= n; k++)
            {
                wa = std::abs (cyv[k]);
                if (wa < std::abs (cyv[k - 1]) )
                {
                    lb = k;
                }
            }
        }
    }
    cyvp[0] = v0 * cyv[0] / z - cyv[1];
    for (k = 1; k <= n; k++)
    {
        cyvp[k] = cyv[k - 1] - (k + v0) * cyv[k] / z;
    }
    vm = n + v0;
    return 0;
}
}