bessik.CPP 12.1 KB
``````//  bessik.cpp -- computation of modified Bessel functions In, Kn
//      and their derivatives. Algorithms and coefficient values from
//      "Computation of Special Functions", Zhang and Jin, John
//      Wiley and Sons, 1996.
//
//
#define _USE_MATH_DEFINES
#include <math.h>
#include "bessel.h"

double gamma(double x);

int bessik01a(double x,double &i0,double &i1,double &k0,double &k1,
double &i0p,double &i1p,double &k0p,double &k1p)
{
double r,x2,ca,cb,ct,ww,w0,xr,xr2;
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};

if (x < 0.0) return 1;
if (x == 0.0) {
i0 = 1.0;
i1 = 0.0;
k0 = 1e308;
k1 = 1e308;
i0p = 0.0;
i1p = 0.5;
k0p = -1e308;
k1p = -1e308;
return 0;
}
x2 = x*x;
if (x <= 18.0) {
i0 = 1.0;
r = 1.0;
for (k=1;k<=50;k++) {
r *= 0.25*x2/(k*k);
i0 += r;
if (fabs(r/i0) < eps) break;
}
i1 = 1.0;
r = 1.0;
for (k=1;k<=50;k++) {
r *= 0.25*x2/(k*(k+1));
i1 += r;
if (fabs(r/i1) < eps) break;
}
i1 *= 0.5*x;
}
else {
if (x >= 50.0) kz = 7;
else if (x >= 35.0) kz = 9;
else kz = 12;
ca = exp(x)/sqrt(2.0*M_PI*x);
i0 = 1.0;
xr = 1.0/x;
for (k=0;k<kz;k++) {
i0 += a[k]*pow(xr,k+1);
}
i0 *= ca;
i1 = 1.0;
for (k=0;k<kz;k++) {
i1 += b[k]*pow(xr,k+1);
}
i1 *= ca;
}
if (x <= 9.0) {
ct = -(log(0.5*x)+el);
k0 = 0.0;
w0 = 0.0;
r = 1.0;
ww = 0.0;
for (k=1;k<=50;k++) {
w0 += 1.0/k;
r *= 0.25*x2/(k*k);
k0 += r*(w0+ct);
if (fabs((k0-ww)/k0) < eps) break;
ww = k0;
}
k0 += ct;
}
else {
cb = 0.5/x;
xr2 = 1.0/x2;
k0 = 1.0;
for (k=0;k<8;k++) {
k0 += a1[k]*pow(xr2,k+1);
}
k0 *= cb/i0;
}
k1 = (1.0/x - i1*k0)/i0;
i0p = i1;
i1p = i0-i1/x;
k0p = -k1;
k1p = -k0-k1/x;
return 0;
}

int bessik01b(double x,double &i0,double &i1,double &k0,double &k1,
double &i0p,double &i1p,double &k0p,double &k1p)
{
double t,t2,dtmp,dtmp1;

if (x < 0.0) return 1;
if (x == 0.0) {
i0 = 1.0;
i1 = 0.0;
k0 = 1e308;
k1 = 1e308;
i0p = 0.0;
i1p = 0.5;
k0p = -1e308;
k1p = -1e308;
return 0;
}
if (x < 3.75) {
t = x/3.75;
t2 = t*t;
i0 = (((((0.0045813*t2+0.0360768)*t2+0.2659732)*t2+
1.2067492)*t2+3.0899424)*t2+3.5156229)*t2+1.0;
i1 = x*(((((0.00032411*t2+0.00301532)*t2+0.02658733*t2+
0.15084934)*t2+0.51498869)*t2+0.87890594)*t2+0.5);
}
else {
t = 3.75/x;
dtmp1 = exp(x)/sqrt(x);
dtmp = (((((((0.00392377*t-0.01647633)*t+0.026355537)*t-0.02057706)*t+
0.00916281)*t-0.00157565)*t+0.00225319)*t+0.01328592)*t+0.39894228;
i0 = dtmp*dtmp1;
dtmp = (((((((-0.00420059*t+0.01787654)*t-0.02895312)*t+0.02282967)*t-
0.01031555)*t+0.00163801)*t-0.00362018)*t-0.03988024)*t+0.39894228;
i1 = dtmp*dtmp1;
}
if (x < 2.0) {
t = 0.5*x;
t2 = t*t;       // already calculated above
dtmp = (((((0.0000074*t2+0.0001075)*t2+0.00262698)*t2+0.0348859)*t2+
0.23069756)*t2+0.4227842)*t2-0.57721566;
k0 = dtmp - i0*log(t);
dtmp = (((((-0.00004686*t2-0.00110404)*t2-0.01919402)*t2-
0.18156897)*t2-0.67278578)*t2+0.15443144)*t2+1.0;
k1 = dtmp/x + i1*log(t);
}
else {
t = 2.0/x;
dtmp1 = exp(-x)/sqrt(x);
dtmp = (((((0.00053208*t-0.0025154)*t+0.00587872)*t-0.01062446)*t+
0.02189568)*t-0.07832358)*t+1.25331414;
k0 = dtmp*dtmp1;
dtmp = (((((-0.00068245*t+0.00325614)*t-0.00780353)*t+0.01504268)*t-
0.0365562)*t+0.23498619)*t+1.25331414;
k1 = dtmp*dtmp1;
}
i0p = i1;
i1p = i0 - i1/x;
k0p = -k1;
k1p = -k0 - k1/x;
return 0;
}
int bessikna(int n,double x,int &nm,double *in,double *kn,
double *inp,double *knp)
{
double bi0,bi1,bk0,bk1,g,g0,g1,f,f0,f1,h,h0,h1,s0;
int k,m,ecode;

if ((x < 0.0) || (n < 0)) return 1;
if (x < eps) {
for (k=0;k<=n;k++) {
in[k] = 0.0;
kn[k] = 1e308;
inp[k] = 0.0;
knp[k] = -1e308;
}
in[0] = 1.0;
inp[1] = 0.5;
return 0;
}
nm = n;
ecode = bessik01a(x,in[0],in[1],kn[0],kn[1],inp[0],inp[1],knp[0],knp[1]);
if (n < 2) return 0;
bi0 = in[0];
bi1 = in[1];
bk0 = kn[0];
bk1 = kn[1];
if ((x > 40.0) && (n < (int)(0.25*x))) {
h0 = bi0;
h1 = bi1;
for (k=2;k<=n;k++) {
h = -2.0*(k-1.0)*h1/x+h0;
in[k] = h;
h0 = h1;
h1 = h;
}
}
else {
m = msta1(x,200);
if (m < n) nm = m;
else m = msta2(x,n,15);
f0 = 0.0;
f1 = 1.0e-100;
for (k=m;k>=0;k--) {
f = 2.0*(k+1.0)*f1/x+f0;
if (x <= nm) in[k] = f;
f0 = f1;
f1 = f;
}
s0 = bi0/f;
for (k=0;k<=m;k++) {
in[k] *= s0;
}
}
g0 = bk0;
g1 = bk1;
for (k=2;k<=nm;k++) {
g = 2.0*(k-1.0)*g1/x+g0;
kn[k] = g;
g0 = g1;
g1 = g;
}
for (k=2;k<=nm;k++) {
inp[k] = in[k-1]-k*in[k]/x;
knp[k] = -kn[k-1]-k*kn[k]/x;
}
return 0;
}
int bessiknb(int n,double x,int &nm,double *in,double *kn,
double *inp,double *knp)
{
double s0,bs,f,f0,f1,sk0,a0,bkl,vt,r,g,g0,g1;
int k,kz,m,l;

if ((x < 0.0) || (n < 0)) return 1;
if (x < eps) {
for (k=0;k<=n;k++) {
in[k] = 0.0;
kn[k] = 1e308;
inp[k] = 0.0;
knp[k] = -1e308;
}
in[0] = 1.0;
inp[1] = 0.5;
return 0;
}
nm = n;
if (n == 0) nm = 1;
m = msta1(x,200);
if (m < nm) nm = m;
else m = msta2(x,nm,15);
bs = 0.0;
sk0 = 0.0;
f0 = 0.0;
f1 = 1.0e-100;
for (k=m;k>=0;k--) {
f = 2.0*(k+1.0)*f1/x+f0;
if (k <= nm) in[k] = f;
if ((k != 0) && (k == 2*(int)(k/2))) {
sk0 += 4.0*f/k;
}
bs += 2.0*f;
f0 = f1;
f1 = f;
}
s0 = exp(x)/(bs-f);
for (k=0;k<=nm;k++) {
in[k] *= s0;
}
if (x <= 8.0) {
kn[0] = -(log(0.5*x)+el)*in[0]+s0*sk0;
kn[1] = (1.0/x-in[1]*kn[0])/in[0];
}
else {
a0 = sqrt(M_PI_2/x)*exp(-x);
if (x >= 200.0) kz = 6;
else if (x >= 80.0) kz = 8;
else if (x >= 25.0) kz = 10;
else kz = 16;
for (l=0;l<2;l++) {
bkl = 1.0;
vt = 4.0*l;
r = 1.0;
for (k=1;k<=kz;k++) {
r *= 0.125*(vt-pow(2.0*k-1.0,2))/(k*x);
bkl += r;
}
kn[l] = a0*bkl;
}
}
g0 = kn[0];
g1 = kn[1];
for (k=2;k<=nm;k++) {
g = 2.0*(k-1.0)*g1/x+g0;
kn[k] = g;
g0 = g1;
g1 = g;
}
inp[0] = in[1];
knp[0] = -kn[1];
for (k=1;k<=nm;k++) {
inp[k] = in[k-1]-k*in[k]/x;
knp[k] = -kn[k-1]-k*kn[k]/x;
}
return 0;
}

//  The following program computes the modified Bessel functions
//  Iv(x) and Kv(x) for arbitrary positive order. For negative
//  order use:
//
//          I-v(x) = Iv(x) + 2/pi sin(v pi) Kv(x)
//          K-v(x) = Kv(x)
//
int bessikv(double v,double x,double &vm,double *iv,double *kv,
double *ivp,double *kvp)
{
double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum;
double f,f1,f2,ct,cs,wa,gan,ww,w0,v0n;
double r1,r2,bk0,bk1,bk2,a2,cb;
int n,k,kz,m;

if ((v < 0.0) || (x < 0.0)) return 1;
x2 = x*x;
n = (int)v;
v0 = v-n;
if (n == 0) n = 1;
if (x == 0.0) {
for (k=0;k<=n;k++) {
iv[k] = 0.0;
kv[k] = -1e308;
ivp[k] = 0.0;
kvp[k] = 1e308;
}
if (v0 == 0.0) {
iv[0] = 1.0;
ivp[1] = 0.5;
}
vm = v;
return 0;
}
piv = M_PI*v0;
vt = 4.0*v0*v0;
if (v0 == 0.0) {
a1 = 1.0;
}
else {
v0p = 1.0+v0;
gap = gamma(v0p);
a1 = pow(0.5*x,v0)/gap;
}
if (x >= 50.0) kz = 8;
else if (x >= 35.0) kz = 10;
else kz = 14;
if (x <= 18.0) {
bi0 = 1.0;
r = 1.0;
for (k=1;k<=30;k++) {
r *= 0.25*x2/(k*(k+v0));
bi0 += r;
if (fabs(r/bi0) < eps) break;
}
bi0 *= a1;
}
else {
ca = exp(x)/sqrt(2.0*M_PI*x);
sum = 1.0;
r = 1.0;
for (k=1;k<=kz;k++) {
r *= -0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
sum += r;
}
bi0 = ca*sum;
}
m = msta1(x,200);
if (m < n) n = m;
else m = msta2(x,n,15);
f2 = 0.0;
f1 = 1.0e-100;
for (k=m;k>=0;k--) {
f = 2.0*(v0+k+1.0)*f1/x+f2;
if (k <= n) iv[k] = f;
f2 = f1;
f1 = f;
}
cs = bi0/f;
for (k=0;k<=n;k++) {
iv[k] *= cs;
}
ivp[0] = v0*iv[0]/x+iv[1];
for (k=1;k<=n;k++) {
ivp[k] = -(k+v0)*iv[k]/x+iv[k-1];
}
ww = 0.0;
if (x <= 9.0) {
if (v0 == 0.0) {
ct = -log(0.5*x)-el;
cs = 0.0;
w0 = 0.0;
r = 1.0;
for (k=1;k<=50;k++) {
w0 += 1.0/k;
r *= 0.25*x2/(k*k);
cs += r*(w0+ct);
wa = fabs(cs);
if (fabs((wa-ww)/wa) < eps) break;
ww = wa;
}
bk0 = ct+cs;
}
else {
v0n = 1.0-v0;
gan = gamma(v0n);
a2 = 1.0/(gan*pow(0.5*x,v0));
a1 = pow(0.5*x,v0)/gap;
sum = a2-a1;
r1 = 1.0;
r2 = 1.0;
for (k=1;k<=120;k++) {
r1 *= 0.25*x2/(k*(k-v0));
r2 *= 0.25*x2/(k*(k+v0));
sum += a2*r1-a1*r2;
wa = fabs(sum);
if (fabs((wa-ww)/wa) < eps) break;
ww = wa;
}
bk0 = M_PI_2*sum/sin(piv);
}
}
else {
cb = exp(-x)*sqrt(M_PI_2/x);
sum = 1.0;
r = 1.0;
for (k=1;k<=kz;k++) {
r *= 0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
sum += r;
}
bk0 = cb*sum;
}
bk1 = (1.0/x-iv[1]*bk0)/iv[0];
kv[0] = bk0;
kv[1] = bk1;
for (k=2;k<=n;k++) {
bk2 = 2.0*(v0+k-1.0)*bk1/x+bk0;
kv[k] = bk2;
bk0 = bk1;
bk1 = bk2;
}
kvp[0] = v0*kv[0]/x-kv[1];
for (k=1;k<=n;k++) {
kvp[k] = -(k+v0)*kv[k]/x-kv[k-1];
}
vm = n+v0;
return 0;
}``````