// 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. // // (C) 2003, C. Bond. All rights reserved. // #define _USE_MATH_DEFINES #include #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 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; }