cbessik.cpp 11.9 KB
``````//  cbessik.cpp -- complex modified Bessel functions.
//  Algorithms and coefficient values from "Computation of Special
//  Functions", Zhang and Jin, John Wiley and Sons, 1996.
//
//
#include <complex>
using namespace std;
#include "bessel.h"

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

double gamma(double x);

int cbessik01(complex<double>z,complex<double>&ci0,complex<double>&ci1,
complex<double>&ck0,complex<double>&ck1,complex<double>&ci0p,
complex<double>&ci1p,complex<double>&ck0p,complex<double>&ck1p)
{
complex<double> z1,z2,zr,zr2,cr,ca,cb,cs,ct,cw;
double a0,w0;
int k,kz;
static double a[] = {
0.125,
7.03125e-2,
7.32421875e-2,
1.1215209960938e-1,
2.2710800170898e-1,
5.7250142097473e-1,
1.7277275025845,
6.0740420012735,
2.4380529699556e1,
1.1001714026925e2,
5.5133589612202e2,
3.0380905109224e3};
static double b[] = {
-0.375,
-1.171875e-1,
-1.025390625e-1,
-1.4419555664063e-1,
-2.7757644653320e-1,
-6.7659258842468e-1,
-1.9935317337513,
-6.8839142681099,
-2.7248827311269e1,
-1.2159789187654e2,
-6.0384407670507e2,
-3.3022722944809e3};
static double a1[] = {
0.125,
0.2109375,
1.0986328125,
1.1775970458984e1,
2.1461706161499e2,
5.9511522710323e3,
2.3347645606175e5,
1.2312234987631e7,
8.401390346421e08,
7.2031420482627e10};

a0 = abs(z);
z2 = z*z;
z1 = z;
if (a0 == 0.0) {
ci0 = cone;
ci1 = czero;
ck0 = complex<double> (1e308,0);
ck1 = complex<double> (1e308,0);
ci0p = czero;
ci1p = complex<double>(0.5,0.0);
ck0p = complex<double>(-1e308,0);
ck1p = complex<double>(-1e308,0);
return 0;
}
if (real(z) < 0.0) z1 = -z;
if (a0 <= 18.0) {
ci0 = cone;
cr = cone;
for (k=1;k<=50;k++) {
cr *= 0.25*z2/(double)(k*k);
ci0 += cr;
if (abs(cr/ci0) < eps) break;
}
ci1 = cone;
cr = cone;
for (k=1;k<=50;k++) {
cr *= 0.25*z2/(double)(k*(k+1.0));
ci1 += cr;
if (abs(cr/ci1) < eps) break;
}
ci1 *= 0.5*z1;
}
else {
if (a0 >= 50.0) kz = 7;
else if (a0 >= 35.0) kz = 9;
else kz = 12;
ca = exp(z1)/sqrt(2.0*M_PI*z1);
ci0 = cone;
zr = 1.0/z1;
for (k=0;k<kz;k++) {
ci0 += a[k]*pow(zr,k+1.0);
}
ci0 *= ca;
ci1 = cone;
for (k=0;k<kz;k++) {
ci1 += b[k]*pow(zr,k+1.0);
}
ci1 *= ca;
}
if (a0 <= 9.0) {
cs = czero;
ct = -log(0.5*z1)-el;
w0 = 0.0;
cr = cone;
for (k=1;k<=50;k++) {
w0 += 1.0/k;
cr *= 0.25*z2/(double)(k*k);
cs += cr*(w0+ct);
if (abs((cs-cw)/cs) < eps) break;
cw = cs;
}
ck0 = ct+cs;
}
else {
cb = 0.5/z1;
zr2 = 1.0/z2;
ck0 = cone;
for (k=0;k<10;k++) {
ck0 += a1[k]*pow(zr2,k+1.0);
}
ck0 *= cb/ci0;
}
ck1 = (1.0/z1 - ci1*ck0)/ci0;
if (real(z) < 0.0) {
if (imag(z) < 0.0) {
ck0 += cii*M_PI*ci0;
ck1 = -ck1+cii*M_PI*ci1;
}
else if (imag(z) > 0.0) {
ck0 -= cii*M_PI*ci0;
ck1 = -ck1-cii*M_PI*ci1;
}
ci1 = -ci1;
}
ci0p = ci1;
ci1p = ci0-1.0*ci1/z;
ck0p = -ck1;
ck1p = -ck0-1.0*ck1/z;
return 0;
}
int cbessikna(int n,complex<double> z,int &nm,complex<double> *ci,
complex<double> *ck,complex<double> *cip,complex<double> *ckp)
{
complex<double> ci0,ci1,ck0,ck1,ckk,cf,cf1,cf2,cs;
double a0;
int k,m,ecode;
a0 = abs(z);
nm = n;
if (a0 < 1.0e-100) {
for (k=0;k<=n;k++) {
ci[k] = czero;
ck[k] = complex<double>(-1e308,0);
cip[k] = czero;
ckp[k] = complex<double>(1e308,0);
}
ci[0] = cone;
cip[1] = complex<double>(0.5,0.0);
return 0;
}
ecode = cbessik01(z,ci[0],ci[1],ck[0],ck[1],cip[0],cip[1],ckp[0],ckp[1]);
if (n < 2) return 0;
ci0 = ci[0];
ci1 = ci[1];
ck0 = ck[0];
ck1 = ck[1];
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) ci[k] = cf;
cf2 = cf1;
cf1 = cf;
}
cs = ci0/cf;
for (k=0;k<=nm;k++) {
ci[k] *= cs;
}
for (k=2;k<=nm;k++) {
if (abs(ci[k-1]) > abs(ci[k-2])) {
ckk = (1.0/z-ci[k]*ck[k-1])/ci[k-1];
}
else {
ckk = (ci[k]*ck[k-2]+2.0*(k-1.0)/(z*z))/ci[k-2];
}
ck[k] = ckk;
}
for (k=2;k<=nm;k++) {
cip[k] = ci[k-1]-(double)k*ci[k]/z;
ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
}
return 0;
}
int cbessiknb(int n,complex<double> z,int &nm,complex<double> *ci,
complex<double> *ck,complex<double> *cip,complex<double> *ckp)
{
complex<double> z1,cbs,csk0,cf,cf0,cf1,ca0,cbkl;
complex<double> cg,cg0,cg1,cs0,cs,cr;
double a0,vt,fac;
int k,kz,l,m;

a0 = abs(z);
nm = n;
if (a0 < 1.0e-100) {
for (k=0;k<=n;k++) {
ci[k] = czero;
ck[k] = complex<double>(1e308,0);
cip[k] = czero;
ckp[k] = complex<double>(-1e308,0);
}
ci[0] = complex<double>(1.0,0.0);
cip[1] = complex<double>(0.5,0.0);
return 0;
}
z1 = z;
if (real(z) < 0.0) z1 = -z;
if (n == 0) nm = 1;
m = msta1(a0,200);
if (m < nm) nm = m;
else m = msta2(a0,nm,15);
cbs = czero;
csk0 = czero;
cf0 = czero;
cf1 = complex<double>(1.0e-100,0.0);
for (k=m;k>=0;k--) {
cf = 2.0*(k+1.0)*cf1/z1+cf0;
if (k <=nm) ci[k] = cf;
if ((k != 0) && (k == 2*(k>>1)))  csk0 += 4.0*cf/(double)k;
cbs += 2.0*cf;
cf0 = cf1;
cf1 = cf;
}
cs0 = exp(z1)/(cbs-cf);
for (k=0;k<=nm;k++) {
ci[k] *= cs0;
}
if (a0 <= 9.0) {
ck[0] = -(log(0.5*z1)+el)*ci[0]+cs0*csk0;
ck[1] = (1.0/z1-ci[1]*ck[0])/ci[0];
}
else {
ca0 = sqrt(M_PI_2/z1)*exp(-z1);
if (a0 >= 200.0) kz = 6;
else if (a0 >= 80.0) kz = 8;
else if (a0 >= 25.0) kz = 10;
else kz = 16;
for (l=0;l<2;l++) {
cbkl = cone;
vt = 4.0*l;
cr = cone;
for (k=1;k<=kz;k++) {
cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
cbkl += cr;
}
ck[l] = ca0*cbkl;
}
}
cg0 = ck[0];
cg1 = ck[1];
for (k=2;k<=nm;k++) {
cg = 2.0*(k-1.0)*cg1/z1+cg0;
ck[k] = cg;
cg0 = cg1;
cg1 = cg;
}
if (real(z) < 0.0) {
fac = 1.0;
for (k=0;k<=nm;k++) {
if (imag(z) < 0.0) {
ck[k] = fac*ck[k]+cii*M_PI*ci[k];
}
else {
ck[k] = fac*ck[k]-cii*M_PI*ci[k];
}
ci[k] *= fac;
fac = -fac;
}
}
cip[0] = ci[1];
ckp[0] = -ck[1];
for (k=1;k<=nm;k++) {
cip[k] = ci[k-1]-(double)k*ci[k]/z;
ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
}
return 0;
}
int cbessikv(double v,complex<double>z,double &vm,complex<double> *civ,
complex<double> *ckv,complex<double> *civp,complex<double> *ckvp)
{
complex<double> z1,z2,ca1,ca,cs,cr,ci0,cbi0,cf,cf1,cf2;
complex<double> ct,cp,cbk0,ca2,cr1,cr2,csu,cws,cb;
complex<double> cg0,cg1,cgk,cbk1,cvk;
double a0,v0,v0p,v0n,vt,w0,piv,gap,gan;
int m,n,k,kz;

a0 = abs(z);
z1 = z;
z2 = z*z;
n = (int)v;
v0 = v-n;
piv = M_PI*v0;
vt = 4.0*v0*v0;
if (n == 0) n = 1;
if (a0 < 1e-100) {
for (k=0;k<=n;k++) {
civ[k] = czero;
ckv[k] = complex<double>(-1e308,0);
civp[k] = czero;
ckvp[k] = complex<double>(1e308,0);
}
if (v0 == 0.0) {
civ[0] = cone;
civp[1] = complex<double> (0.5,0.0);
}
vm = v;
return 0;
}
if (a0 >= 50.0) kz = 8;
else if (a0 >= 35.0) kz = 10;
else kz = 14;
if (real(z) <= 0.0) z1 = -z;
if (a0 < 18.0) {
if (v0 == 0.0) {
ca1 = cone;
}
else {
v0p = 1.0+v0;
gap = gamma(v0p);
ca1 = pow(0.5*z1,v0)/gap;
}
ci0 = cone;
cr = cone;
for (k=1;k<=50;k++) {
cr *= 0.25*z2/(k*(k+v0));
ci0 += cr;
if (abs(cr/ci0) < eps) break;
}
cbi0 = ci0*ca1;
}
else {
ca = exp(z1)/sqrt(2.0*M_PI*z1);
cs = cone;
cr = cone;
for (k=1;k<=kz;k++) {
cr *= -0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
cs += cr;
}
cbi0 = ca*cs;
}
m = msta1(a0,200);
if (m < n) n = 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*(v0+k+1.0)*cf1/z1+cf2;
if (k <= n) civ[k] = cf;
cf2 = cf1;
cf1 = cf;
}
cs = cbi0/cf;
for (k=0;k<=n;k++) {
civ[k] *= cs;
}
if (a0 <= 9.0) {
if (v0 == 0.0) {
ct = -log(0.5*z1)-el;
cs = czero;
w0 = 0.0;
cr = cone;
for (k=1;k<=50;k++) {
w0 += 1.0/k;
cr *= 0.25*z2/(double)(k*k);
cp = cr*(w0+ct);
cs += cp;
if ((k >= 10) && (abs(cp/cs) < eps)) break;
}
cbk0 = ct+cs;
}
else {
v0n = 1.0-v0;
gan = gamma(v0n);
ca2 = 1.0/(gan*pow(0.5*z1,v0));
ca1 = pow(0.5*z1,v0)/gap;
csu = ca2-ca1;
cr1 = cone;
cr2 = cone;
cws = czero;
for (k=1;k<=50;k++) {
cr1 *= 0.25*z2/(k*(k-v0));
cr2 *= 0.25*z2/(k*(k+v0));
csu += ca2*cr1-ca1*cr2;
if ((k >= 10) && (abs((cws-csu)/csu) < eps)) break;
cws = csu;
}
cbk0 = csu*M_PI_2/sin(piv);
}
}
else {
cb = exp(-z1)*sqrt(M_PI_2/z1);
cs = cone;
cr = cone;
for (k=1;k<=kz;k++) {
cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
cs += cr;
}
cbk0 = cb*cs;
}
cbk1 = (1.0/z1-civ[1]*cbk0)/civ[0];
ckv[0] = cbk0;
ckv[1] = cbk1;
cg0 = cbk0;
cg1 = cbk1;
for (k=2;k<=n;k++) {
cgk = 2.0*(v0+k-1.0)*cg1/z1+cg0;
ckv[k] = cgk;
cg0 = cg1;
cg1 = cgk;
}
if (real(z) < 0.0) {
for (k=0;k<=n;k++) {
cvk = exp((k+v0)*M_PI*cii);
if (imag(z) < 0.0) {
ckv[k] = cvk*ckv[k]+M_PI*cii*civ[k];
civ[k] /= cvk;
}
else if (imag(z) > 0.0) {
ckv[k] = ckv[k]/cvk-M_PI*cii*civ[k];
civ[k] *= cvk;
}
}
}
civp[0] = v0*civ[0]/z+civ[1];
ckvp[0] = v0*ckv[0]/z-ckv[1];
for (k=1;k<=n;k++) {
civp[k] = -(k+v0)*civ[k]/z+civ[k-1];
ckvp[k] = -(k+v0)*ckv[k]/z-ckv[k-1];
}
vm = n+v0;
return 0;
}``````