// cbessik.cpp -- complex modified 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 using namespace std; #include "bessel.h" static complex cii(0.0,1.0); static complex czero(0.0,0.0); static complex cone(1.0,0.0); double gamma(double x); int cbessik01(complexz,complex&ci0,complex&ci1, complex&ck0,complex&ck1,complex&ci0p, complex&ci1p,complex&ck0p,complex&ck1p) { complex 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 (1e308,0); ck1 = complex (1e308,0); ci0p = czero; ci1p = complex(0.5,0.0); ck0p = complex(-1e308,0); ck1p = complex(-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 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 z,int &nm,complex *ci, complex *ck,complex *cip,complex *ckp) { complex 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(-1e308,0); cip[k] = czero; ckp[k] = complex(1e308,0); } ci[0] = cone; cip[1] = complex(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(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 z,int &nm,complex *ci, complex *ck,complex *cip,complex *ckp) { complex z1,cbs,csk0,cf,cf0,cf1,ca0,cbkl; complex 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(1e308,0); cip[k] = czero; ckp[k] = complex(-1e308,0); } ci[0] = complex(1.0,0.0); cip[1] = complex(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(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,complexz,double &vm,complex *civ, complex *ckv,complex *civp,complex *ckvp) { complex z1,z2,ca1,ca,cs,cr,ci0,cbi0,cf,cf1,cf2; complex ct,cp,cbk0,ca2,cr1,cr2,csu,cws,cb; complex 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(-1e308,0); civp[k] = czero; ckvp[k] = complex(1e308,0); } if (v0 == 0.0) { civ[0] = cone; civp[1] = complex (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(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; }