optics / imie

Commit da3d4e0e3dfabd920e06f7e3efc307a3d4b9059c

Authored by dmayerich
0 parents

Initial commit.

Showing 29 changed files with 12860 additions and 0 deletions   Show diff stats
Inline Side-by-side
BESSEL.h 0 → 100644
  Show/Hide comments   View file @da3d4e0
  1 +++ a/BESSEL.h
  1 +#ifndef bessH
  2 +#define bessH
  3 +#define _USE_MATH_DEFINES
  4 +#include <math.h>
  5 +#include <complex>
  6 +using namespace std;
  7 +#define eps 1e-15
  8 +#define el 0.5772156649015329
  9 +
  10 +int msta1(double x,int mp);
  11 +int msta2(double x,int n,int mp);
  12 +int bessjy01a(double x,double &j0,double &j1,double &y0,double &y1,
  13 + double &j0p,double &j1p,double &y0p,double &y1p);
  14 +int bessjy01b(double x,double &j0,double &j1,double &y0,double &y1,
  15 + double &j0p,double &j1p,double &y0p,double &y1p);
  16 +int bessjyna(int n,double x,int &nm,double *jn,double *yn,
  17 + double *jnp,double *ynp);
  18 +int bessjynb(int n,double x,int &nm,double *jn,double *yn,
  19 + double *jnp,double *ynp);
  20 +int bessjyv(double v,double x,double &vm,double *jv,double *yv,
  21 + double *jvp,double *yvp);
  22 +int bessik01a(double x,double &i0,double &i1,double &k0,double &k1,
  23 + double &i0p,double &i1p,double &k0p,double &k1p);
  24 +int bessik01b(double x,double &i0,double &i1,double &k0,double &k1,
  25 + double &i0p,double &i1p,double &k0p,double &k1p);
  26 +int bessikna(int n,double x,int &nm,double *in,double *kn,
  27 + double *inp,double *knp);
  28 +int bessiknb(int n,double x,int &nm,double *in,double *kn,
  29 + double *inp,double *knp);
  30 +int bessikv(double v,double x,double &vm,double *iv,double *kv,
  31 + double *ivp,double *kvp);
  32 +int cbessjy01(complex<double> z,complex<double> &cj0,complex<double> &cj1,
  33 + complex<double> &cy0,complex<double> &cy1,complex<double> &cj0p,
  34 + complex<double> &cj1p,complex<double> &cy0p,complex<double> &cy1p);
  35 +int cbessjyna(int n,complex<double> z,int &nm,complex<double> *cj,
  36 + complex<double> *cy,complex<double> *cjp,complex<double> *cyp);
  37 +int cbessjynb(int n,complex<double> z,int &nm,complex<double> *cj,
  38 + complex<double> *cy,complex<double> *cjp,complex<double> *cyp);
  39 +int cbessik01(complex<double>z,complex<double>&ci0,complex<double>&ci1,
  40 + complex<double>&ck0,complex<double>&ck1,complex<double>&ci0p,
  41 + complex<double>&ci1p,complex<double>&ck0p,complex<double>&ck1p);
  42 +int cbessikna(int n,complex<double> z,int &nm,complex<double> *ci,
  43 + complex<double> *ck,complex<double> *cip,complex<double> *ckp);
  44 +int cbessiknb(int n,complex<double> z,int &nm,complex<double> *ci,
  45 + complex<double> *ck,complex<double> *cip,complex<double> *ckp);
  46 +int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
  47 + complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
  48 +int cbessikv(double v,complex<double>z,double &vm,complex<double> *civ,
  49 + complex<double> *ckv,complex<double> *civp,complex<double> *ckvp);
  50 +
  51 +#endif
... ...
BESSIK.CPP 0 → 100644
  Show/Hide comments   View file @da3d4e0
  1 +++ a/BESSIK.CPP
  1 +// bessik.cpp -- computation of modified Bessel functions In, Kn
  2 +// and their derivatives. Algorithms and coefficient values from
  3 +// "Computation of Special Functions", Zhang and Jin, John
  4 +// Wiley and Sons, 1996.
  5 +//
  6 +// (C) 2003, C. Bond. All rights reserved.
  7 +//
  8 +#define _USE_MATH_DEFINES
  9 +#include <math.h>
  10 +#include "bessel.h"
  11 +
  12 +double gamma(double x);
  13 +
  14 +int bessik01a(double x,double &i0,double &i1,double &k0,double &k1,
  15 + double &i0p,double &i1p,double &k0p,double &k1p)
  16 +{
  17 + double r,x2,ca,cb,ct,ww,w0,xr,xr2;
  18 + int k,kz;
  19 + static double a[] = {
  20 + 0.125,
  21 + 7.03125e-2,
  22 + 7.32421875e-2,
  23 + 1.1215209960938e-1,
  24 + 2.2710800170898e-1,
  25 + 5.7250142097473e-1,
  26 + 1.7277275025845,
  27 + 6.0740420012735,
  28 + 2.4380529699556e1,
  29 + 1.1001714026925e2,
  30 + 5.5133589612202e2,
  31 + 3.0380905109224e3};
  32 + static double b[] = {
  33 + -0.375,
  34 + -1.171875e-1,
  35 + -1.025390625e-1,
  36 + -1.4419555664063e-1,
  37 + -2.7757644653320e-1,
  38 + -6.7659258842468e-1,
  39 + -1.9935317337513,
  40 + -6.8839142681099,
  41 + -2.7248827311269e1,
  42 + -1.2159789187654e2,
  43 + -6.0384407670507e2,
  44 + -3.3022722944809e3};
  45 + static double a1[] = {
  46 + 0.125,
  47 + 0.2109375,
  48 + 1.0986328125,
  49 + 1.1775970458984e1,
  50 + 2.1461706161499e2,
  51 + 5.9511522710323e3,
  52 + 2.3347645606175e5,
  53 + 1.2312234987631e7};
  54 +
  55 + if (x < 0.0) return 1;
  56 + if (x == 0.0) {
  57 + i0 = 1.0;
  58 + i1 = 0.0;
  59 + k0 = 1e308;
  60 + k1 = 1e308;
  61 + i0p = 0.0;
  62 + i1p = 0.5;
  63 + k0p = -1e308;
  64 + k1p = -1e308;
  65 + return 0;
  66 + }
  67 + x2 = x*x;
  68 + if (x <= 18.0) {
  69 + i0 = 1.0;
  70 + r = 1.0;
  71 + for (k=1;k<=50;k++) {
  72 + r *= 0.25*x2/(k*k);
  73 + i0 += r;
  74 + if (fabs(r/i0) < eps) break;
  75 + }
  76 + i1 = 1.0;
  77 + r = 1.0;
  78 + for (k=1;k<=50;k++) {
  79 + r *= 0.25*x2/(k*(k+1));
  80 + i1 += r;
  81 + if (fabs(r/i1) < eps) break;
  82 + }
  83 + i1 *= 0.5*x;
  84 + }
  85 + else {
  86 + if (x >= 50.0) kz = 7;
  87 + else if (x >= 35.0) kz = 9;
  88 + else kz = 12;
  89 + ca = exp(x)/sqrt(2.0*M_PI*x);
  90 + i0 = 1.0;
  91 + xr = 1.0/x;
  92 + for (k=0;k<kz;k++) {
  93 + i0 += a[k]*pow(xr,k+1);
  94 + }
  95 + i0 *= ca;
  96 + i1 = 1.0;
  97 + for (k=0;k<kz;k++) {
  98 + i1 += b[k]*pow(xr,k+1);
  99 + }
  100 + i1 *= ca;
  101 + }
  102 + if (x <= 9.0) {
  103 + ct = -(log(0.5*x)+el);
  104 + k0 = 0.0;
  105 + w0 = 0.0;
  106 + r = 1.0;
  107 + ww = 0.0;
  108 + for (k=1;k<=50;k++) {
  109 + w0 += 1.0/k;
  110 + r *= 0.25*x2/(k*k);
  111 + k0 += r*(w0+ct);
  112 + if (fabs((k0-ww)/k0) < eps) break;
  113 + ww = k0;
  114 + }
  115 + k0 += ct;
  116 + }
  117 + else {
  118 + cb = 0.5/x;
  119 + xr2 = 1.0/x2;
  120 + k0 = 1.0;
  121 + for (k=0;k<8;k++) {
  122 + k0 += a1[k]*pow(xr2,k+1);
  123 + }
  124 + k0 *= cb/i0;
  125 + }
  126 + k1 = (1.0/x - i1*k0)/i0;
  127 + i0p = i1;
  128 + i1p = i0-i1/x;
  129 + k0p = -k1;
  130 + k1p = -k0-k1/x;
  131 + return 0;
  132 +}
  133 +
  134 +int bessik01b(double x,double &i0,double &i1,double &k0,double &k1,
  135 + double &i0p,double &i1p,double &k0p,double &k1p)
  136 +{
  137 + double t,t2,dtmp,dtmp1;
  138 +
  139 + if (x < 0.0) return 1;
  140 + if (x == 0.0) {
  141 + i0 = 1.0;
  142 + i1 = 0.0;
  143 + k0 = 1e308;
  144 + k1 = 1e308;
  145 + i0p = 0.0;
  146 + i1p = 0.5;
  147 + k0p = -1e308;
  148 + k1p = -1e308;
  149 + return 0;
  150 + }
  151 + if (x < 3.75) {
  152 + t = x/3.75;
  153 + t2 = t*t;
  154 + i0 = (((((0.0045813*t2+0.0360768)*t2+0.2659732)*t2+
  155 + 1.2067492)*t2+3.0899424)*t2+3.5156229)*t2+1.0;
  156 + i1 = x*(((((0.00032411*t2+0.00301532)*t2+0.02658733*t2+
  157 + 0.15084934)*t2+0.51498869)*t2+0.87890594)*t2+0.5);
  158 + }
  159 + else {
  160 + t = 3.75/x;
  161 + dtmp1 = exp(x)/sqrt(x);
  162 + dtmp = (((((((0.00392377*t-0.01647633)*t+0.026355537)*t-0.02057706)*t+
  163 + 0.00916281)*t-0.00157565)*t+0.00225319)*t+0.01328592)*t+0.39894228;
  164 + i0 = dtmp*dtmp1;
  165 + dtmp = (((((((-0.00420059*t+0.01787654)*t-0.02895312)*t+0.02282967)*t-
  166 + 0.01031555)*t+0.00163801)*t-0.00362018)*t-0.03988024)*t+0.39894228;
  167 + i1 = dtmp*dtmp1;
  168 + }
  169 + if (x < 2.0) {
  170 + t = 0.5*x;
  171 + t2 = t*t; // already calculated above
  172 + dtmp = (((((0.0000074*t2+0.0001075)*t2+0.00262698)*t2+0.0348859)*t2+
  173 + 0.23069756)*t2+0.4227842)*t2-0.57721566;
  174 + k0 = dtmp - i0*log(t);
  175 + dtmp = (((((-0.00004686*t2-0.00110404)*t2-0.01919402)*t2-
  176 + 0.18156897)*t2-0.67278578)*t2+0.15443144)*t2+1.0;
  177 + k1 = dtmp/x + i1*log(t);
  178 + }
  179 + else {
  180 + t = 2.0/x;
  181 + dtmp1 = exp(-x)/sqrt(x);
  182 + dtmp = (((((0.00053208*t-0.0025154)*t+0.00587872)*t-0.01062446)*t+
  183 + 0.02189568)*t-0.07832358)*t+1.25331414;
  184 + k0 = dtmp*dtmp1;
  185 + dtmp = (((((-0.00068245*t+0.00325614)*t-0.00780353)*t+0.01504268)*t-
  186 + 0.0365562)*t+0.23498619)*t+1.25331414;
  187 + k1 = dtmp*dtmp1;
  188 + }
  189 + i0p = i1;
  190 + i1p = i0 - i1/x;
  191 + k0p = -k1;
  192 + k1p = -k0 - k1/x;
  193 + return 0;
  194 +}
  195 +int bessikna(int n,double x,int &nm,double *in,double *kn,
  196 + double *inp,double *knp)
  197 +{
  198 + double bi0,bi1,bk0,bk1,g,g0,g1,f,f0,f1,h,h0,h1,s0;
  199 + int k,m,ecode;
  200 +
  201 + if ((x < 0.0) || (n < 0)) return 1;
  202 + if (x < eps) {
  203 + for (k=0;k<=n;k++) {
  204 + in[k] = 0.0;
  205 + kn[k] = 1e308;
  206 + inp[k] = 0.0;
  207 + knp[k] = -1e308;
  208 + }
  209 + in[0] = 1.0;
  210 + inp[1] = 0.5;
  211 + return 0;
  212 + }
  213 + nm = n;
  214 + ecode = bessik01a(x,in[0],in[1],kn[0],kn[1],inp[0],inp[1],knp[0],knp[1]);
  215 + if (n < 2) return 0;
  216 + bi0 = in[0];
  217 + bi1 = in[1];
  218 + bk0 = kn[0];
  219 + bk1 = kn[1];
  220 + if ((x > 40.0) && (n < (int)(0.25*x))) {
  221 + h0 = bi0;
  222 + h1 = bi1;
  223 + for (k=2;k<=n;k++) {
  224 + h = -2.0*(k-1.0)*h1/x+h0;
  225 + in[k] = h;
  226 + h0 = h1;
  227 + h1 = h;
  228 + }
  229 + }
  230 + else {
  231 + m = msta1(x,200);
  232 + if (m < n) nm = m;
  233 + else m = msta2(x,n,15);
  234 + f0 = 0.0;
  235 + f1 = 1.0e-100;
  236 + for (k=m;k>=0;k--) {
  237 + f = 2.0*(k+1.0)*f1/x+f0;
  238 + if (x <= nm) in[k] = f;
  239 + f0 = f1;
  240 + f1 = f;
  241 + }
  242 + s0 = bi0/f;
  243 + for (k=0;k<=m;k++) {
  244 + in[k] *= s0;
  245 + }
  246 + }
  247 + g0 = bk0;
  248 + g1 = bk1;
  249 + for (k=2;k<=nm;k++) {
  250 + g = 2.0*(k-1.0)*g1/x+g0;
  251 + kn[k] = g;
  252 + g0 = g1;
  253 + g1 = g;
  254 + }
  255 + for (k=2;k<=nm;k++) {
  256 + inp[k] = in[k-1]-k*in[k]/x;
  257 + knp[k] = -kn[k-1]-k*kn[k]/x;
  258 + }
  259 + return 0;
  260 +}
  261 +int bessiknb(int n,double x,int &nm,double *in,double *kn,
  262 + double *inp,double *knp)
  263 +{
  264 + double s0,bs,f,f0,f1,sk0,a0,bkl,vt,r,g,g0,g1;
  265 + int k,kz,m,l;
  266 +
  267 + if ((x < 0.0) || (n < 0)) return 1;
  268 + if (x < eps) {
  269 + for (k=0;k<=n;k++) {
  270 + in[k] = 0.0;
  271 + kn[k] = 1e308;
  272 + inp[k] = 0.0;
  273 + knp[k] = -1e308;
  274 + }
  275 + in[0] = 1.0;
  276 + inp[1] = 0.5;
  277 + return 0;
  278 + }
  279 + nm = n;
  280 + if (n == 0) nm = 1;
  281 + m = msta1(x,200);
  282 + if (m < nm) nm = m;
  283 + else m = msta2(x,nm,15);
  284 + bs = 0.0;
  285 + sk0 = 0.0;
  286 + f0 = 0.0;
  287 + f1 = 1.0e-100;
  288 + for (k=m;k>=0;k--) {
  289 + f = 2.0*(k+1.0)*f1/x+f0;
  290 + if (k <= nm) in[k] = f;
  291 + if ((k != 0) && (k == 2*(int)(k/2))) {
  292 + sk0 += 4.0*f/k;
  293 + }
  294 + bs += 2.0*f;
  295 + f0 = f1;
  296 + f1 = f;
  297 + }
  298 + s0 = exp(x)/(bs-f);
  299 + for (k=0;k<=nm;k++) {
  300 + in[k] *= s0;
  301 + }
  302 + if (x <= 8.0) {
  303 + kn[0] = -(log(0.5*x)+el)*in[0]+s0*sk0;
  304 + kn[1] = (1.0/x-in[1]*kn[0])/in[0];
  305 + }
  306 + else {
  307 + a0 = sqrt(M_PI_2/x)*exp(-x);
  308 + if (x >= 200.0) kz = 6;
  309 + else if (x >= 80.0) kz = 8;
  310 + else if (x >= 25.0) kz = 10;
  311 + else kz = 16;
  312 + for (l=0;l<2;l++) {
  313 + bkl = 1.0;
  314 + vt = 4.0*l;
  315 + r = 1.0;
  316 + for (k=1;k<=kz;k++) {
  317 + r *= 0.125*(vt-pow(2.0*k-1.0,2))/(k*x);
  318 + bkl += r;
  319 + }
  320 + kn[l] = a0*bkl;
  321 + }
  322 + }
  323 + g0 = kn[0];
  324 + g1 = kn[1];
  325 + for (k=2;k<=nm;k++) {
  326 + g = 2.0*(k-1.0)*g1/x+g0;
  327 + kn[k] = g;
  328 + g0 = g1;
  329 + g1 = g;
  330 + }
  331 + inp[0] = in[1];
  332 + knp[0] = -kn[1];
  333 + for (k=1;k<=nm;k++) {
  334 + inp[k] = in[k-1]-k*in[k]/x;
  335 + knp[k] = -kn[k-1]-k*kn[k]/x;
  336 + }
  337 + return 0;
  338 +}
  339 +
  340 +// The following program computes the modified Bessel functions
  341 +// Iv(x) and Kv(x) for arbitrary positive order. For negative
  342 +// order use:
  343 +//
  344 +// I-v(x) = Iv(x) + 2/pi sin(v pi) Kv(x)
  345 +// K-v(x) = Kv(x)
  346 +//
  347 +int bessikv(double v,double x,double &vm,double *iv,double *kv,
  348 + double *ivp,double *kvp)
  349 +{
  350 + double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum;
  351 + double f,f0,f1,f2,ct,cs,wa,gan,ww,w0,v0n;
  352 + double r1,r2,bk0,bk1,bk2,a2,cb;
  353 + int n,k,kz,m;
  354 +
  355 + if ((v < 0.0) || (x < 0.0)) return 1;
  356 + x2 = x*x;
  357 + n = (int)v;
  358 + v0 = v-n;
  359 + if (n == 0) n = 1;
  360 + if (x == 0.0) {
  361 + for (k=0;k<=n;k++) {
  362 + iv[k] = 0.0;
  363 + kv[k] = -1e308;
  364 + ivp[k] = 0.0;
  365 + kvp[k] = 1e308;
  366 + }
  367 + if (v0 == 0.0) {
  368 + iv[0] = 1.0;
  369 + ivp[1] = 0.5;
  370 + }
  371 + vm = v;
  372 + return 0;
  373 + }
  374 + piv = M_PI*v0;
  375 + vt = 4.0*v0*v0;
  376 + if (v0 == 0.0) {
  377 + a1 = 1.0;
  378 + }
  379 + else {
  380 + v0p = 1.0+v0;
  381 + gap = gamma(v0p);
  382 + a1 = pow(0.5*x,v0)/gap;
  383 + }
  384 + if (x >= 50.0) kz = 8;
  385 + else if (x >= 35.0) kz = 10;
  386 + else kz = 14;
  387 + if (x <= 18.0) {
  388 + bi0 = 1.0;
  389 + r = 1.0;
  390 + for (k=1;k<=30;k++) {
  391 + r *= 0.25*x2/(k*(k+v0));
  392 + bi0 += r;
  393 + if (fabs(r/bi0) < eps) break;
  394 + }
  395 + bi0 *= a1;
  396 + }
  397 + else {
  398 + ca = exp(x)/sqrt(2.0*M_PI*x);
  399 + sum = 1.0;
  400 + r = 1.0;
  401 + for (k=1;k<=kz;k++) {
  402 + r *= -0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
  403 + sum += r;
  404 + }
  405 + bi0 = ca*sum;
  406 + }
  407 + m = msta1(x,200);
  408 + if (m < n) n = m;
  409 + else m = msta2(x,n,15);
  410 + f2 = 0.0;
  411 + f1 = 1.0e-100;
  412 + for (k=m;k>=0;k--) {
  413 + f = 2.0*(v0+k+1.0)*f1/x+f2;
  414 + if (k <= n) iv[k] = f;
  415 + f2 = f1;
  416 + f1 = f;
  417 + }
  418 + cs = bi0/f;
  419 + for (k=0;k<=n;k++) {
  420 + iv[k] *= cs;
  421 + }
  422 + ivp[0] = v0*iv[0]/x+iv[1];
  423 + for (k=1;k<=n;k++) {
  424 + ivp[k] = -(k+v0)*iv[k]/x+iv[k-1];
  425 + }
  426 + ww = 0.0;
  427 + if (x <= 9.0) {
  428 + if (v0 == 0.0) {
  429 + ct = -log(0.5*x)-el;
  430 + cs = 0.0;
  431 + w0 = 0.0;
  432 + r = 1.0;
  433 + for (k=1;k<=50;k++) {
  434 + w0 += 1.0/k;
  435 + r *= 0.25*x2/(k*k);
  436 + cs += r*(w0+ct);
  437 + wa = fabs(cs);
  438 + if (fabs((wa-ww)/wa) < eps) break;
  439 + ww = wa;
  440 + }
  441 + bk0 = ct+cs;
  442 + }
  443 + else {
  444 + v0n = 1.0-v0;
  445 + gan = gamma(v0n);
  446 + a2 = 1.0/(gan*pow(0.5*x,v0));
  447 + a1 = pow(0.5*x,v0)/gap;
  448 + sum = a2-a1;
  449 + r1 = 1.0;
  450 + r2 = 1.0;
  451 + for (k=1;k<=120;k++) {
  452 + r1 *= 0.25*x2/(k*(k-v0));
  453 + r2 *= 0.25*x2/(k*(k+v0));
  454 + sum += a2*r1-a1*r2;
  455 + wa = fabs(sum);
  456 + if (fabs((wa-ww)/wa) < eps) break;
  457 + ww = wa;
  458 + }
  459 + bk0 = M_PI_2*sum/sin(piv);
  460 + }
  461 + }
  462 + else {
  463 + cb = exp(-x)*sqrt(M_PI_2/x);
  464 + sum = 1.0;
  465 + r = 1.0;
  466 + for (k=1;k<=kz;k++) {
  467 + r *= 0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
  468 + sum += r;
  469 + }
  470 + bk0 = cb*sum;
  471 + }
  472 + bk1 = (1.0/x-iv[1]*bk0)/iv[0];
  473 + kv[0] = bk0;
  474 + kv[1] = bk1;
  475 + for (k=2;k<=n;k++) {
  476 + bk2 = 2.0*(v0+k-1.0)*bk1/x+bk0;
  477 + kv[k] = bk2;
  478 + bk0 = bk1;
  479 + bk1 = bk2;
  480 + }
  481 + kvp[0] = v0*kv[0]/x-kv[1];
  482 + for (k=1;k<=n;k++) {
  483 + kvp[k] = -(k+v0)*kv[k]/x-kv[k-1];
  484 + }
  485 + vm = n+v0;
  486 + return 0;
  487 +}
... ...
BESSJY.CPP 0 → 100644
  Show/Hide comments   View file @da3d4e0
  1 +++ a/BESSJY.CPP
  1 +// bessjy.cpp -- computation of Bessel functions Jn, Yn and their
  2 +// derivatives. Algorithms and coefficient values from
  3 +// "Computation of Special Functions", Zhang and Jin, John
  4 +// Wiley and Sons, 1996.
  5 +//
  6 +// (C) 2003, C. Bond. All rights reserved.
  7 +//
  8 +// Note that 'math.h' provides (or should provide) values for:
  9 +// pi M_PI
  10 +// 2/pi M_2_PI
  11 +// pi/4 M_PI_4
  12 +// pi/2 M_PI_2
  13 +//
  14 +#define _USE_MATH_DEFINES
  15 +#include <math.h>
  16 +#include "bessel.h"
  17 +
  18 +double gamma(double x);
  19 +//
  20 +// INPUT:
  21 +// double x -- argument of Bessel function
  22 +//
  23 +// OUTPUT (via address pointers):
  24 +// double j0 -- Bessel function of 1st kind, 0th order
  25 +// double j1 -- Bessel function of 1st kind, 1st order
  26 +// double y0 -- Bessel function of 2nd kind, 0th order
  27 +// double y1 -- Bessel function of 2nd kind, 1st order
  28 +// double j0p -- derivative of Bessel function of 1st kind, 0th order
  29 +// double j1p -- derivative of Bessel function of 1st kind, 1st order
  30 +// double y0p -- derivative of Bessel function of 2nd kind, 0th order
  31 +// double y1p -- derivative of Bessel function of 2nd kind, 1st order
  32 +//
  33 +// RETURN:
  34 +// int error code: 0 = OK, 1 = error
  35 +//
  36 +// This algorithm computes the above functions using series expansions.
  37 +//
  38 +int bessjy01a(double x,double &j0,double &j1,double &y0,double &y1,
  39 + double &j0p,double &j1p,double &y0p,double &y1p)
  40 +{
  41 + double x2,r,ec,w0,w1,r0,r1,cs0,cs1;
  42 + double cu,p0,q0,p1,q1,t1,t2;
  43 + int k,kz;
  44 + static double a[] = {
  45 + -7.03125e-2,
  46 + 0.112152099609375,
  47 + -0.5725014209747314,
  48 + 6.074042001273483,
  49 + -1.100171402692467e2,
  50 + 3.038090510922384e3,
  51 + -1.188384262567832e5,
  52 + 6.252951493434797e6,
  53 + -4.259392165047669e8,
  54 + 3.646840080706556e10,
  55 + -3.833534661393944e12,
  56 + 4.854014686852901e14,
  57 + -7.286857349377656e16,
  58 + 1.279721941975975e19};
  59 + static double b[] = {
  60 + 7.32421875e-2,
  61 + -0.2271080017089844,
  62 + 1.727727502584457,
  63 + -2.438052969955606e1,
  64 + 5.513358961220206e2,
  65 + -1.825775547429318e4,
  66 + 8.328593040162893e5,
  67 + -5.006958953198893e7,
  68 + 3.836255180230433e9,
  69 + -3.649010818849833e11,
  70 + 4.218971570284096e13,
  71 + -5.827244631566907e15,
  72 + 9.476288099260110e17,
  73 + -1.792162323051699e20};
  74 + static double a1[] = {
  75 + 0.1171875,
  76 + -0.1441955566406250,
  77 + 0.6765925884246826,
  78 + -6.883914268109947,
  79 + 1.215978918765359e2,
  80 + -3.302272294480852e3,
  81 + 1.276412726461746e5,
  82 + -6.656367718817688e6,
  83 + 4.502786003050393e8,
  84 + -3.833857520742790e10,
  85 + 4.011838599133198e12,
  86 + -5.060568503314727e14,
  87 + 7.572616461117958e16,
  88 + -1.326257285320556e19};
  89 + static double b1[] = {
  90 + -0.1025390625,
  91 + 0.2775764465332031,
  92 + -1.993531733751297,
  93 + 2.724882731126854e1,
  94 + -6.038440767050702e2,
  95 + 1.971837591223663e4,
  96 + -8.902978767070678e5,
  97 + 5.310411010968522e7,
  98 + -4.043620325107754e9,
  99 + 3.827011346598605e11,
  100 + -4.406481417852278e13,
  101 + 6.065091351222699e15,
  102 + -9.833883876590679e17,
  103 + 1.855045211579828e20};
  104 +
  105 + if (x < 0.0) return 1;
  106 + if (x == 0.0) {
  107 + j0 = 1.0;
  108 + j1 = 0.0;
  109 + y0 = -1e308;
  110 + y1 = -1e308;
  111 + j0p = 0.0;
  112 + j1p = 0.5;
  113 + y0p = 1e308;
  114 + y1p = 1e308;
  115 + return 0;
  116 + }
  117 + x2 = x*x;
  118 + if (x <= 12.0) {
  119 + j0 = 1.0;
  120 + r = 1.0;
  121 + for (k=1;k<=30;k++) {
  122 + r *= -0.25*x2/(k*k);
  123 + j0 += r;
  124 + if (fabs(r) < fabs(j0)*1e-15) break;
  125 + }
  126 + j1 = 1.0;
  127 + r = 1.0;
  128 + for (k=1;k<=30;k++) {
  129 + r *= -0.25*x2/(k*(k+1));
  130 + j1 += r;
  131 + if (fabs(r) < fabs(j1)*1e-15) break;
  132 + }
  133 + j1 *= 0.5*x;
  134 + ec = log(0.5*x)+el;
  135 + cs0 = 0.0;
  136 + w0 = 0.0;
  137 + r0 = 1.0;
  138 + for (k=1;k<=30;k++) {
  139 + w0 += 1.0/k;
  140 + r0 *= -0.25*x2/(k*k);
  141 + r = r0 * w0;
  142 + cs0 += r;
  143 + if (fabs(r) < fabs(cs0)*1e-15) break;
  144 + }
  145 + y0 = M_2_PI*(ec*j0-cs0);
  146 + cs1 = 1.0;
  147 + w1 = 0.0;
  148 + r1 = 1.0;
  149 + for (k=1;k<=30;k++) {
  150 + w1 += 1.0/k;
  151 + r1 *= -0.25*x2/(k*(k+1));
  152 + r = r1*(2.0*w1+1.0/(k+1));
  153 + cs1 += r;
  154 + if (fabs(r) < fabs(cs1)*1e-15) break;
  155 + }
  156 + y1 = M_2_PI * (ec*j1-1.0/x-0.25*x*cs1);
  157 + }
  158 + else {
  159 + if (x >= 50.0) kz = 8; // Can be changed to 10
  160 + else if (x >= 35.0) kz = 10; // " " 12
  161 + else kz = 12; // " " 14
  162 + t1 = x-M_PI_4;
  163 + p0 = 1.0;
  164 + q0 = -0.125/x;
  165 + for (k=0;k<kz;k++) {
  166 + p0 += a[k]*pow(x,-2*k-2);
  167 + q0 += b[k]*pow(x,-2*k-3);
  168 + }
  169 + cu = sqrt(M_2_PI/x);
  170 + j0 = cu*(p0*cos(t1)-q0*sin(t1));
  171 + y0 = cu*(p0*sin(t1)+q0*cos(t1));
  172 + t2 = x-0.75*M_PI;
  173 + p1 = 1.0;
  174 + q1 = 0.375/x;
  175 + for (k=0;k<kz;k++) {
  176 + p1 += a1[k]*pow(x,-2*k-2);
  177 + q1 += b1[k]*pow(x,-2*k-3);
  178 + }
  179 + j1 = cu*(p1*cos(t2)-q1*sin(t2));
  180 + y1 = cu*(p1*sin(t2)+q1*cos(t2));
  181 + }
  182 + j0p = -j1;
  183 + j1p = j0-j1/x;
  184 + y0p = -y1;
  185 + y1p = y0-y1/x;
  186 + return 0;
  187 +}
  188 +//
  189 +// INPUT:
  190 +// double x -- argument of Bessel function
  191 +//
  192 +// OUTPUT:
  193 +// double j0 -- Bessel function of 1st kind, 0th order
  194 +// double j1 -- Bessel function of 1st kind, 1st order
  195 +// double y0 -- Bessel function of 2nd kind, 0th order
  196 +// double y1 -- Bessel function of 2nd kind, 1st order
  197 +// double j0p -- derivative of Bessel function of 1st kind, 0th order
  198 +// double j1p -- derivative of Bessel function of 1st kind, 1st order
  199 +// double y0p -- derivative of Bessel function of 2nd kind, 0th order
  200 +// double y1p -- derivative of Bessel function of 2nd kind, 1st order
  201 +//
  202 +// RETURN:
  203 +// int error code: 0 = OK, 1 = error
  204 +//
  205 +// This algorithm computes the functions using polynomial approximations.
  206 +//
  207 +int bessjy01b(double x,double &j0,double &j1,double &y0,double &y1,
  208 + double &j0p,double &j1p,double &y0p,double &y1p)
  209 +{
  210 + double t,t2,dtmp,a0,p0,q0,p1,q1,ta0,ta1;
  211 + if (x < 0.0) return 1;
  212 + if (x == 0.0) {
  213 + j0 = 1.0;
  214 + j1 = 0.0;
  215 + y0 = -1e308;
  216 + y1 = -1e308;
  217 + j0p = 0.0;
  218 + j1p = 0.5;
  219 + y0p = 1e308;
  220 + y1p = 1e308;
  221 + return 0;
  222 + }
  223 + if(x <= 4.0) {
  224 + t = x/4.0;
  225 + t2 = t*t;
  226 + j0 = ((((((-0.5014415e-3*t2+0.76771853e-2)*t2-0.0709253492)*t2+
  227 + 0.4443584263)*t2-1.7777560599)*t2+3.9999973021)*t2
  228 + -3.9999998721)*t2+1.0;
  229 + j1 = t*(((((((-0.1289769e-3*t2+0.22069155e-2)*t2-0.0236616773)*t2+
  230 + 0.1777582922)*t2-0.8888839649)*t2+2.6666660544)*t2-
  231 + 3.999999971)*t2+1.9999999998);
  232 + dtmp = (((((((-0.567433e-4*t2+0.859977e-3)*t2-0.94855882e-2)*t2+
  233 + 0.0772975809)*t2-0.4261737419)*t2+1.4216421221)*t2-
  234 + 2.3498519931)*t2+1.0766115157)*t2+0.3674669052;
  235 + y0 = M_2_PI*log(0.5*x)*j0+dtmp;
  236 + dtmp = (((((((0.6535773e-3*t2-0.0108175626)*t2+0.107657607)*t2-
  237 + 0.7268945577)*t2+3.1261399273)*t2-7.3980241381)*t2+
  238 + 6.8529236342)*t2+0.3932562018)*t2-0.6366197726;
  239 + y1 = M_2_PI*log(0.5*x)*j1+dtmp/x;
  240 + }
  241 + else {
  242 + t = 4.0/x;
  243 + t2 = t*t;
  244 + a0 = sqrt(M_2_PI/x);
  245 + p0 = ((((-0.9285e-5*t2+0.43506e-4)*t2-0.122226e-3)*t2+
  246 + 0.434725e-3)*t2-0.4394275e-2)*t2+0.999999997;
  247 + q0 = t*(((((0.8099e-5*t2-0.35614e-4)*t2+0.85844e-4)*t2-
  248 + 0.218024e-3)*t2+0.1144106e-2)*t2-0.031249995);
  249 + ta0 = x-M_PI_4;
  250 + j0 = a0*(p0*cos(ta0)-q0*sin(ta0));
  251 + y0 = a0*(p0*sin(ta0)+q0*cos(ta0));
  252 + p1 = ((((0.10632e-4*t2-0.50363e-4)*t2+0.145575e-3)*t2
  253 + -0.559487e-3)*t2+0.7323931e-2)*t2+1.000000004;
  254 + q1 = t*(((((-0.9173e-5*t2+0.40658e-4)*t2-0.99941e-4)*t2
  255 + +0.266891e-3)*t2-0.1601836e-2)*t2+0.093749994);
  256 + ta1 = x-0.75*M_PI;
  257 + j1 = a0*(p1*cos(ta1)-q1*sin(ta1));
  258 + y1 = a0*(p1*sin(ta1)+q1*cos(ta1));
  259 + }
  260 + j0p = -j1;
  261 + j1p = j0-j1/x;
  262 + y0p = -y1;
  263 + y1p = y0-y1/x;
  264 + return 0;
  265 +}
  266 +int msta1(double x,int mp)
  267 +{
  268 + double a0,f0,f1,f;
  269 + int i,n0,n1,nn;
  270 +
  271 + a0 = fabs(x);
  272 + n0 = (int)(1.1*a0)+1;
  273 + f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-mp;
  274 + n1 = n0+5;
  275 + f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp;
  276 + for (i=0;i<20;i++) {
  277 + nn = n1-(n1-n0)/(1.0-f0/f1);
  278 + f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp;
  279 + if (abs(nn-n1) < 1) break;
  280 + n0 = n1;
  281 + f0 = f1;
  282 + n1 = nn;
  283 + f1 = f;
  284 + }
  285 + return nn;
  286 +}
  287 +int msta2(double x,int n,int mp)
  288 +{
  289 + double a0,ejn,hmp,f0,f1,f,obj;
  290 + int i,n0,n1,nn;
  291 +
  292 + a0 = fabs(x);
  293 + hmp = 0.5*mp;
  294 + ejn = 0.5*log10(6.28*n)-n*log10(1.36*a0/n);
  295 + if (ejn <= hmp) {
  296 + obj = mp;
  297 + n0 = (int)(1.1*a0);
  298 + if (n0 < 1) n0 = 1;
  299 + }
  300 + else {
  301 + obj = hmp+ejn;
  302 + n0 = n;
  303 + }
  304 + f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-obj;
  305 + n1 = n0+5;
  306 + f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj;
  307 + for (i=0;i<20;i++) {
  308 + nn = n1-(n1-n0)/(1.0-f0/f1);
  309 + f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj;
  310 + if (abs(nn-n1) < 1) break;
  311 + n0 = n1;
  312 + f0 = f1;
  313 + n1 = nn;
  314 + f1 = f;
  315 + }
  316 + return nn+10;
  317 +}
  318 +//
  319 +// INPUT:
  320 +// double x -- argument of Bessel function of 1st and 2nd kind.
  321 +// int n -- order
  322 +//
  323 +// OUPUT:
  324 +//
  325 +// int nm -- highest order actually computed (nm <= n)
  326 +// double jn[] -- Bessel function of 1st kind, orders from 0 to nm
  327 +// double yn[] -- Bessel function of 2nd kind, orders from 0 to nm
  328 +// double j'n[]-- derivative of Bessel function of 1st kind,
  329 +// orders from 0 to nm
  330 +// double y'n[]-- derivative of Bessel function of 2nd kind,
  331 +// orders from 0 to nm
  332 +//
  333 +// Computes Bessel functions of all order up to 'n' using recurrence
  334 +// relations. If 'nm' < 'n' only 'nm' orders are returned.
  335 +//
  336 +int bessjyna(int n,double x,int &nm,double *jn,double *yn,
  337 + double *jnp,double *ynp)
  338 +{
  339 + double bj0,bj1,f,f0,f1,f2,cs;
  340 + int i,k,m,ecode;
  341 +
  342 + nm = n;
  343 + if ((x < 0.0) || (n < 0)) return 1;
  344 + if (x < 1e-15) {
  345 + for (i=0;i<=n;i++) {
  346 + jn[i] = 0.0;
  347 + yn[i] = -1e308;
  348 + jnp[i] = 0.0;
  349 + ynp[i] = 1e308;
  350 + }
  351 + jn[0] = 1.0;
  352 + jnp[1] = 0.5;
  353 + return 0;
  354 + }
  355 + ecode = bessjy01a(x,jn[0],jn[1],yn[0],yn[1],jnp[0],jnp[1],ynp[0],ynp[1]);
  356 + if (n < 2) return 0;
  357 + bj0 = jn[0];
  358 + bj1 = jn[1];
  359 + if (n < (int)0.9*x) {
  360 + for (k=2;k<=n;k++) {
  361 + jn[k] = 2.0*(k-1.0)*bj1/x-bj0;
  362 + bj0 = bj1;
  363 + bj1 = jn[k];
  364 + }
  365 + }
  366 + else {
  367 + m = msta1(x,200);
  368 + if (m < n) nm = m;
  369 + else m = msta2(x,n,15);
  370 + f2 = 0.0;
  371 + f1 = 1.0e-100;
  372 + for (k=m;k>=0;k--) {
  373 + f = 2.0*(k+1.0)/x*f1-f2;
  374 + if (k <= nm) jn[k] = f;
  375 + f2 = f1;
  376 + f1 = f;
  377 + }
  378 + if (fabs(bj0) > fabs(bj1)) cs = bj0/f;
  379 + else cs = bj1/f2;
  380 + for (k=0;k<=nm;k++) {
  381 + jn[k] *= cs;
  382 + }
  383 + }
  384 + for (k=2;k<=nm;k++) {
  385 + jnp[k] = jn[k-1]-k*jn[k]/x;
  386 + }
  387 + f0 = yn[0];
  388 + f1 = yn[1];
  389 + for (k=2;k<=nm;k++) {
  390 + f = 2.0*(k-1.0)*f1/x-f0;
  391 + yn[k] = f;
  392 + f0 = f1;
  393 + f1 = f;
  394 + }
  395 + for (k=2;k<=nm;k++) {
  396 + ynp[k] = yn[k-1]-k*yn[k]/x;
  397 + }
  398 + return 0;
  399 +}
  400 +//
  401 +// Same input and output conventions as above. Different recurrence
  402 +// relations used for 'x' < 300.
  403 +//
  404 +int bessjynb(int n,double x,int &nm,double *jn,double *yn,
  405 + double *jnp,double *ynp)
  406 +{
  407 + double t1,t2,f,f0,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
  408 + double ec,bs,byk,p0,p1,q0,q1;
  409 + static double a[] = {
  410 + -0.7031250000000000e-1,
  411 + 0.1121520996093750,
  412 + -0.5725014209747314,
  413 + 6.074042001273483};
  414 + static double b[] = {
  415 + 0.7324218750000000e-1,
  416 + -0.2271080017089844,
  417 + 1.727727502584457,
  418 + -2.438052969955606e1};
  419 + static double a1[] = {
  420 + 0.1171875,
  421 + -0.1441955566406250,
  422 + 0.6765925884246826,
  423 + -6.883914268109947};
  424 + static double b1[] = {
  425 + -0.1025390625,
  426 + 0.2775764465332031,
  427 + -1.993531733751297,
  428 + 2.724882731126854e1};
  429 +
  430 + int i,k,m;
  431 + nm = n;
  432 + if ((x < 0.0) || (n < 0)) return 1;
  433 + if (x < 1e-15) {
  434 + for (i=0;i<=n;i++) {
  435 + jn[i] = 0.0;
  436 + yn[i] = -1e308;
  437 + jnp[i] = 0.0;
  438 + ynp[i] = 1e308;
  439 + }
  440 + jn[0] = 1.0;
  441 + jnp[1] = 0.5;
  442 + return 0;
  443 + }
  444 + if (x <= 300.0 || n > (int)(0.9*x)) {
  445 + if (n == 0) nm = 1;
  446 + m = msta1(x,200);
  447 + if (m < nm) nm = m;
  448 + else m = msta2(x,nm,15);
  449 + bs = 0.0;
  450 + su = 0.0;
  451 + sv = 0.0;
  452 + f2 = 0.0;
  453 + f1 = 1.0e-100;
  454 + for (k = m;k>=0;k--) {
  455 + f = 2.0*(k+1.0)/x*f1 - f2;
  456 + if (k <= nm) jn[k] = f;
  457 + if ((k == 2*(int)(k/2)) && (k != 0)) {
  458 + bs += 2.0*f;
  459 +// su += pow(-1,k>>1)*f/(double)k;
  460 + su += (-1)*((k & 2)-1)*f/(double)k;
  461 + }
  462 + else if (k > 1) {
  463 +// sv += pow(-1,k>>1)*k*f/(k*k-1.0);
  464 + sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0);
  465 + }
  466 + f2 = f1;
  467 + f1 = f;
  468 + }
  469 + s0 = bs+f;
  470 + for (k=0;k<=nm;k++) {
  471 + jn[k] /= s0;
  472 + }
  473 + ec = log(0.5*x) +0.5772156649015329;
  474 + by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
  475 + yn[0] = by0;
  476 + by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
  477 + yn[1] = by1;
  478 + }
  479 + else {
  480 + t1 = x-M_PI_4;
  481 + p0 = 1.0;
  482 + q0 = -0.125/x;
  483 + for (k=0;k<4;k++) {
  484 + p0 += a[k]*pow(x,-2*k-2);
  485 + q0 += b[k]*pow(x,-2*k-3);
  486 + }
  487 + cu = sqrt(M_2_PI/x);
  488 + bj0 = cu*(p0*cos(t1)-q0*sin(t1));
  489 + by0 = cu*(p0*sin(t1)+q0*cos(t1));
  490 + jn[0] = bj0;
  491 + yn[0] = by0;
  492 + t2 = x-0.75*M_PI;
  493 + p1 = 1.0;
  494 + q1 = 0.375/x;
  495 + for (k=0;k<4;k++) {
  496 + p1 += a1[k]*pow(x,-2*k-2);
  497 + q1 += b1[k]*pow(x,-2*k-3);
  498 + }
  499 + bj1 = cu*(p1*cos(t2)-q1*sin(t2));
  500 + by1 = cu*(p1*sin(t2)+q1*cos(t2));
  501 + jn[1] = bj1;
  502 + yn[1] = by1;
  503 + for (k=2;k<=nm;k++) {
  504 + bjk = 2.0*(k-1.0)*bj1/x-bj0;
  505 + jn[k] = bjk;
  506 + bj0 = bj1;
  507 + bj1 = bjk;
  508 + }
  509 + }
  510 + jnp[0] = -jn[1];
  511 + for (k=1;k<=nm;k++) {
  512 + jnp[k] = jn[k-1]-k*jn[k]/x;
  513 + }
  514 + for (k=2;k<=nm;k++) {
  515 + byk = 2.0*(k-1.0)*by1/x-by0;
  516 + yn[k] = byk;
  517 + by0 = by1;
  518 + by1 = byk;
  519 + }
  520 + ynp[0] = -yn[1];
  521 + for (k=1;k<=nm;k++) {
  522 + ynp[k] = yn[k-1]-k*yn[k]/x;
  523 + }
  524 + return 0;
  525 +
  526 +}
  527 +
  528 +// The following routine computes Bessel Jv(x) and Yv(x) for
  529 +// arbitrary positive order (v). For negative order, use:
  530 +//
  531 +// J-v(x) = Jv(x)cos(v pi) - Yv(x)sin(v pi)
  532 +// Y-v(x) = Jv(x)sin(v pi) + Yv(x)cos(v pi)
  533 +//
  534 +int bessjyv(double v,double x,double &vm,double *jv,double *yv,
  535 + double *djv,double *dyv)
  536 +{
  537 + double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2;
  538 + double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq;
  539 + double b,ec,w0,w1,bjy0,bjy1,bju0,bju1,pv0,pv1,byvk;
  540 + int j,k,l,m,n,kz;
  541 +
  542 + x2 = x*x;
  543 + n = (int)v;
  544 + v0 = v-n;
  545 + if ((x < 0.0) || (v < 0.0)) return 1;
  546 + if (x < 1e-15) {
  547 + for (k=0;k<=n;k++) {
  548 + jv[k] = 0.0;
  549 + yv[k] = -1e308;
  550 + djv[k] = 0.0;
  551 + dyv[k] = 1e308;
  552 + if (v0 == 0.0) {
  553 + jv[0] = 1.0;
  554 + djv[1] = 0.5;
  555 + }
  556 + else djv[0] = 1e308;
  557 + }
  558 + vm = v;
  559 + return 0;
  560 + }
  561 + if (x <= 12.0) {
  562 + for (l=0;l<2;l++) {
  563 + vl = v0 + l;
  564 + bjvl = 1.0;
  565 + r = 1.0;
  566 + for (k=1;k<=40;k++) {
  567 + r *= -0.25*x2/(k*(k+vl));
  568 + bjvl += r;
  569 + if (fabs(r) < fabs(bjvl)*1e-15) break;
  570 + }
  571 + vg = 1.0 + vl;
  572 + a = pow(0.5*x,vl)/gamma(vg);
  573 + if (l == 0) bjv0 = bjvl*a;
  574 + else bjv1 = bjvl*a;
  575 + }
  576 + }
  577 + else {
  578 + if (x >= 50.0) kz = 8;
  579 + else if (x >= 35.0) kz = 10;
  580 + else kz = 11;
  581 + for (j=0;j<2;j++) {
  582 + vv = 4.0*(j+v0)*(j+v0);
  583 + px = 1.0;
  584 + rp = 1.0;
  585 + for (k=1;k<=kz;k++) {
  586 + rp *= (-0.78125e-2)*(vv-pow(4.0*k-3.0,2.0))*
  587 + (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*x2);
  588 + px += rp;
  589 + }
  590 + qx = 1.0;
  591 + rq = 1.0;
  592 + for (k=1;k<=kz;k++) {
  593 + rq *= (-0.78125e-2)*(vv-pow(4.0*k-1.0,2.0))*
  594 + (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*x2);
  595 + qx += rq;
  596 + }
  597 + qx *= 0.125*(vv-1.0)/x;
  598 + xk = x-(0.5*(j+v0)+0.25)*M_PI;
  599 + a0 = sqrt(M_2_PI/x);
  600 + ck = cos(xk);
  601 + sk = sin(xk);
  602 +
  603 + if (j == 0) {
  604 + bjv0 = a0*(px*ck-qx*sk);
  605 + byv0 = a0*(px*sk+qx*ck);
  606 + }
  607 + else if (j == 1) {
  608 + bjv1 = a0*(px*ck-qx*sk);
  609 + byv1 = a0*(px*sk+qx*ck);
  610 + }
  611 + }
  612 + }
  613 + jv[0] = bjv0;
  614 + jv[1] = bjv1;
  615 + djv[0] = v0*jv[0]/x-jv[1];
  616 + djv[1] = -(1.0+v0)*jv[1]/x+jv[0];
  617 + if ((n >= 2) && (n <= (int)(0.9*x))) {
  618 + f0 = bjv0;
  619 + f1 = bjv1;
  620 + for (k=2;k<=n;k++) {
  621 + f = 2.0*(k+v0-1.0)*f1/x-f0;
  622 + jv[k] = f;
  623 + f0 = f1;
  624 + f1 = f;
  625 + }
  626 + }
  627 + else if (n >= 2) {
  628 + m = msta1(x,200);
  629 + if (m < n) n = m;
  630 + else m = msta2(x,n,15);
  631 + f2 = 0.0;
  632 + f1 = 1.0e-100;
  633 + for (k=m;k>=0;k--) {
  634 + f = 2.0*(v0+k+1.0)*f1/x-f2;
  635 + if (k <= n) jv[k] = f;
  636 + f2 = f1;
  637 + f1 = f;
  638 + }
  639 + if (fabs(bjv0) > fabs(bjv1)) cs = bjv0/f;
  640 + else cs = bjv1/f2;
  641 + for (k=0;k<=n;k++) {
  642 + jv[k] *= cs;
  643 + }
  644 + }
  645 + for (k=2;k<=n;k++) {
  646 + djv[k] = -(k+v0)*jv[k]/x+jv[k-1];
  647 + }
  648 + if (x <= 12.0) {
  649 + if (v0 != 0.0) {
  650 + for (l=0;l<2;l++) {
  651 + vl = v0 +l;
  652 + bjvl = 1.0;
  653 + r = 1.0;
  654 + for (k=1;k<=40;k++) {
  655 + r *= -0.25*x2/(k*(k-vl));
  656 + bjvl += r;
  657 + if (fabs(r) < fabs(bjvl)*1e-15) break;
  658 + }
  659 + vg = 1.0-vl;
  660 + b = pow(2.0/x,vl)/gamma(vg);
  661 + if (l == 0) bju0 = bjvl*b;
  662 + else bju1 = bjvl*b;
  663 + }
  664 + pv0 = M_PI*v0;
  665 + pv1 = M_PI*(1.0+v0);
  666 + byv0 = (bjv0*cos(pv0)-bju0)/sin(pv0);
  667 + byv1 = (bjv1*cos(pv1)-bju1)/sin(pv1);
  668 + }
  669 + else {
  670 + ec = log(0.5*x)+el;
  671 + cs0 = 0.0;
  672 + w0 = 0.0;
  673 + r0 = 1.0;
  674 + for (k=1;k<=30;k++) {
  675 + w0 += 1.0/k;
  676 + r0 *= -0.25*x2/(k*k);
  677 + cs0 += r0*w0;
  678 + }
  679 + byv0 = M_2_PI*(ec*bjv0-cs0);
  680 + cs1 = 1.0;
  681 + w1 = 0.0;
  682 + r1 = 1.0;
  683 + for (k=1;k<=30;k++) {
  684 + w1 += 1.0/k;
  685 + r1 *= -0.25*x2/(k*(k+1));
  686 + cs1 += r1*(2.0*w1+1.0/(k+1.0));
  687 + }
  688 + byv1 = M_2_PI*(ec*bjv1-1.0/x-0.25*x*cs1);
  689 + }
  690 + }
  691 + yv[0] = byv0;
  692 + yv[1] = byv1;
  693 + for (k=2;k<=n;k++) {
  694 + byvk = 2.0*(v0+k-1.0)*byv1/x-byv0;
  695 + yv[k] = byvk;
  696 + byv0 = byv1;
  697 + byv1 = byvk;
  698 + }
  699 + dyv[0] = v0*yv[0]/x-yv[1];
  700 + for (k=1;k<=n;k++) {
  701 + dyv[k] = -(k+v0)*yv[k]/x+yv[k-1];
  702 + }
  703 + vm = n + v0;
  704 + return 0;
  705 +}
  706 +
... ...
CBESSIK.CPP 0 → 100644
  Show/Hide comments   View file @da3d4e0
  1 +++ a/CBESSIK.CPP
  1 +// cbessik.cpp -- complex modified Bessel functions.
  2 +// Algorithms and coefficient values from "Computation of Special
  3 +// Functions", Zhang and Jin, John Wiley and Sons, 1996.
  4 +//
  5 +// (C) 2003, C. Bond. All rights reserved.
  6 +//
  7 +#include <complex>
  8 +using namespace std;
  9 +#include "bessel.h"
  10 +
  11 +static complex<double> cii(0.0,1.0);
  12 +static complex<double> czero(0.0,0.0);
  13 +static complex<double> cone(1.0,0.0);
  14 +
  15 +double gamma(double x);
  16 +
  17 +int cbessik01(complex<double>z,complex<double>&ci0,complex<double>&ci1,
  18 + complex<double>&ck0,complex<double>&ck1,complex<double>&ci0p,
  19 + complex<double>&ci1p,complex<double>&ck0p,complex<double>&ck1p)
  20 +{
  21 + complex<double> z1,z2,zr,zr2,cr,ca,cb,cs,ct,cw;
  22 + double a0,w0;
  23 + int k,kz;
  24 + static double a[] = {
  25 + 0.125,
  26 + 7.03125e-2,
  27 + 7.32421875e-2,
  28 + 1.1215209960938e-1,
  29 + 2.2710800170898e-1,
  30 + 5.7250142097473e-1,
  31 + 1.7277275025845,
  32 + 6.0740420012735,
  33 + 2.4380529699556e1,
  34 + 1.1001714026925e2,
  35 + 5.5133589612202e2,
  36 + 3.0380905109224e3};
  37 + static double b[] = {
  38 + -0.375,
  39 + -1.171875e-1,
  40 + -1.025390625e-1,
  41 + -1.4419555664063e-1,
  42 + -2.7757644653320e-1,
  43 + -6.7659258842468e-1,
  44 + -1.9935317337513,
  45 + -6.8839142681099,
  46 + -2.7248827311269e1,
  47 + -1.2159789187654e2,
  48 + -6.0384407670507e2,
  49 + -3.3022722944809e3};
  50 + static double a1[] = {
  51 + 0.125,
  52 + 0.2109375,
  53 + 1.0986328125,
  54 + 1.1775970458984e1,
  55 + 2.1461706161499e2,
  56 + 5.9511522710323e3,
  57 + 2.3347645606175e5,
  58 + 1.2312234987631e7,
  59 + 8.401390346421e08,
  60 + 7.2031420482627e10};
  61 +
  62 + a0 = abs(z);
  63 + z2 = z*z;
  64 + z1 = z;
  65 + if (a0 == 0.0) {
  66 + ci0 = cone;
  67 + ci1 = czero;
  68 + ck0 = complex<double> (1e308,0);
  69 + ck1 = complex<double> (1e308,0);
  70 + ci0p = czero;
  71 + ci1p = complex<double>(0.5,0.0);
  72 + ck0p = complex<double>(-1e308,0);
  73 + ck1p = complex<double>(-1e308,0);
  74 + return 0;
  75 + }
  76 + if (real(z) < 0.0) z1 = -z;
  77 + if (a0 <= 18.0) {
  78 + ci0 = cone;
  79 + cr = cone;
  80 + for (k=1;k<=50;k++) {
  81 + cr *= 0.25*z2/(double)(k*k);
  82 + ci0 += cr;
  83 + if (abs(cr/ci0) < eps) break;
  84 + }
  85 + ci1 = cone;
  86 + cr = cone;
  87 + for (k=1;k<=50;k++) {
  88 + cr *= 0.25*z2/(double)(k*(k+1.0));
  89 + ci1 += cr;
  90 + if (abs(cr/ci1) < eps) break;
  91 + }
  92 + ci1 *= 0.5*z1;
  93 + }
  94 + else {
  95 + if (a0 >= 50.0) kz = 7;
  96 + else if (a0 >= 35.0) kz = 9;
  97 + else kz = 12;
  98 + ca = exp(z1)/sqrt(2.0*M_PI*z1);
  99 + ci0 = cone;
  100 + zr = 1.0/z1;
  101 + for (k=0;k<kz;k++) {
  102 + ci0 += a[k]*pow(zr,k+1.0);
  103 + }
  104 + ci0 *= ca;
  105 + ci1 = cone;
  106 + for (k=0;k<kz;k++) {
  107 + ci1 += b[k]*pow(zr,k+1.0);
  108 + }
  109 + ci1 *= ca;
  110 + }
  111 + if (a0 <= 9.0) {
  112 + cs = czero;
  113 + ct = -log(0.5*z1)-el;
  114 + w0 = 0.0;
  115 + cr = cone;
  116 + for (k=1;k<=50;k++) {
  117 + w0 += 1.0/k;
  118 + cr *= 0.25*z2/(double)(k*k);
  119 + cs += cr*(w0+ct);
  120 + if (abs((cs-cw)/cs) < eps) break;
  121 + cw = cs;
  122 + }
  123 + ck0 = ct+cs;
  124 + }
  125 + else {
  126 + cb = 0.5/z1;
  127 + zr2 = 1.0/z2;
  128 + ck0 = cone;
  129 + for (k=0;k<10;k++) {
  130 + ck0 += a1[k]*pow(zr2,k+1.0);
  131 + }
  132 + ck0 *= cb/ci0;
  133 + }
  134 + ck1 = (1.0/z1 - ci1*ck0)/ci0;
  135 + if (real(z) < 0.0) {
  136 + if (imag(z) < 0.0) {
  137 + ck0 += cii*M_PI*ci0;
  138 + ck1 = -ck1+cii*M_PI*ci1;
  139 + }
  140 + else if (imag(z) > 0.0) {
  141 + ck0 -= cii*M_PI*ci0;
  142 + ck1 = -ck1-cii*M_PI*ci1;
  143 + }
  144 + ci1 = -ci1;
  145 + }
  146 + ci0p = ci1;
  147 + ci1p = ci0-1.0*ci1/z;
  148 + ck0p = -ck1;
  149 + ck1p = -ck0-1.0*ck1/z;
  150 + return 0;
  151 +}
  152 +int cbessikna(int n,complex<double> z,int &nm,complex<double> *ci,
  153 + complex<double> *ck,complex<double> *cip,complex<double> *ckp)
  154 +{
  155 + complex<double> ci0,ci1,ck0,ck1,ckk,cf,cf1,cf2,cs;
  156 + double a0;
  157 + int k,m,ecode;
  158 + a0 = abs(z);
  159 + nm = n;
  160 + if (a0 < 1.0e-100) {
  161 + for (k=0;k<=n;k++) {
  162 + ci[k] = czero;
  163 + ck[k] = complex<double>(-1e308,0);
  164 + cip[k] = czero;
  165 + ckp[k] = complex<double>(1e308,0);
  166 + }
  167 + ci[0] = cone;
  168 + cip[1] = complex<double>(0.5,0.0);
  169 + return 0;
  170 + }
  171 + ecode = cbessik01(z,ci[0],ci[1],ck[0],ck[1],cip[0],cip[1],ckp[0],ckp[1]);
  172 + if (n < 2) return 0;
  173 + ci0 = ci[0];
  174 + ci1 = ci[1];
  175 + ck0 = ck[0];
  176 + ck1 = ck[1];
  177 + m = msta1(a0,200);
  178 + if (m < n) nm = m;
  179 + else m = msta2(a0,n,15);
  180 + cf2 = czero;
  181 + cf1 = complex<double>(1.0e-100,0.0);
  182 + for (k=m;k>=0;k--) {
  183 + cf = 2.0*(k+1.0)*cf1/z+cf2;
  184 + if (k <= nm) ci[k] = cf;
  185 + cf2 = cf1;
  186 + cf1 = cf;
  187 + }
  188 + cs = ci0/cf;
  189 + for (k=0;k<=nm;k++) {
  190 + ci[k] *= cs;
  191 + }
  192 + for (k=2;k<=nm;k++) {
  193 + if (abs(ci[k-1]) > abs(ci[k-2])) {
  194 + ckk = (1.0/z-ci[k]*ck[k-1])/ci[k-1];
  195 + }
  196 + else {
  197 + ckk = (ci[k]*ck[k-2]+2.0*(k-1.0)/(z*z))/ci[k-2];
  198 + }
  199 + ck[k] = ckk;
  200 + }
  201 + for (k=2;k<=nm;k++) {
  202 + cip[k] = ci[k-1]-(double)k*ci[k]/z;
  203 + ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
  204 + }
  205 + return 0;
  206 +}
  207 +int cbessiknb(int n,complex<double> z,int &nm,complex<double> *ci,
  208 + complex<double> *ck,complex<double> *cip,complex<double> *ckp)
  209 +{
  210 + complex<double> z1,cbs,csk0,cf,cf0,cf1,ca0,cbkl;
  211 + complex<double> cg,cg0,cg1,cs0,cs,cr;
  212 + double a0,vt,fac;
  213 + int k,kz,l,m;
  214 +
  215 + a0 = abs(z);
  216 + nm = n;
  217 + if (a0 < 1.0e-100) {
  218 + for (k=0;k<=n;k++) {
  219 + ci[k] = czero;
  220 + ck[k] = complex<double>(1e308,0);
  221 + cip[k] = czero;
  222 + ckp[k] = complex<double>(-1e308,0);
  223 + }
  224 + ci[0] = complex<double>(1.0,0.0);
  225 + cip[1] = complex<double>(0.5,0.0);
  226 + return 0;
  227 + }
  228 + z1 = z;
  229 + if (real(z) < 0.0) z1 = -z;
  230 + if (n == 0) nm = 1;
  231 + m = msta1(a0,200);
  232 + if (m < nm) nm = m;
  233 + else m = msta2(a0,nm,15);
  234 + cbs = czero;
  235 + csk0 = czero;
  236 + cf0 = czero;
  237 + cf1 = complex<double>(1.0e-100,0.0);
  238 + for (k=m;k>=0;k--) {
  239 + cf = 2.0*(k+1.0)*cf1/z1+cf0;
  240 + if (k <=nm) ci[k] = cf;
  241 + if ((k != 0) && (k == 2*(k>>1))) csk0 += 4.0*cf/(double)k;
  242 + cbs += 2.0*cf;
  243 + cf0 = cf1;
  244 + cf1 = cf;
  245 + }
  246 + cs0 = exp(z1)/(cbs-cf);
  247 + for (k=0;k<=nm;k++) {
  248 + ci[k] *= cs0;
  249 + }
  250 + if (a0 <= 9.0) {
  251 + ck[0] = -(log(0.5*z1)+el)*ci[0]+cs0*csk0;
  252 + ck[1] = (1.0/z1-ci[1]*ck[0])/ci[0];
  253 + }
  254 + else {
  255 + ca0 = sqrt(M_PI_2/z1)*exp(-z1);
  256 + if (a0 >= 200.0) kz = 6;
  257 + else if (a0 >= 80.0) kz = 8;
  258 + else if (a0 >= 25.0) kz = 10;
  259 + else kz = 16;
  260 + for (l=0;l<2;l++) {
  261 + cbkl = cone;
  262 + vt = 4.0*l;
  263 + cr = cone;
  264 + for (k=1;k<=kz;k++) {
  265 + cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
  266 + cbkl += cr;
  267 + }
  268 + ck[l] = ca0*cbkl;
  269 + }
  270 + }
  271 + cg0 = ck[0];
  272 + cg1 = ck[1];
  273 + for (k=2;k<=nm;k++) {
  274 + cg = 2.0*(k-1.0)*cg1/z1+cg0;
  275 + ck[k] = cg;
  276 + cg0 = cg1;
  277 + cg1 = cg;
  278 + }
  279 + if (real(z) < 0.0) {
  280 + fac = 1.0;
  281 + for (k=0;k<=nm;k++) {
  282 + if (imag(z) < 0.0) {
  283 + ck[k] = fac*ck[k]+cii*M_PI*ci[k];
  284 + }
  285 + else {
  286 + ck[k] = fac*ck[k]-cii*M_PI*ci[k];
  287 + }
  288 + ci[k] *= fac;
  289 + fac = -fac;
  290 + }
  291 + }
  292 + cip[0] = ci[1];
  293 + ckp[0] = -ck[1];
  294 + for (k=1;k<=nm;k++) {
  295 + cip[k] = ci[k-1]-(double)k*ci[k]/z;
  296 + ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
  297 + }
  298 + return 0;
  299 +}
  300 +int cbessikv(double v,complex<double>z,double &vm,complex<double> *civ,
  301 + complex<double> *ckv,complex<double> *civp,complex<double> *ckvp)
  302 +{
  303 + complex<double> z1,z2,ca1,ca,cs,cr,ci0,cbi0,cf,cf1,cf2;
  304 + complex<double> ct,cp,cbk0,ca2,cr1,cr2,csu,cws,cb;
  305 + complex<double> cg0,cg1,cgk,cbk1,cvk;
  306 + double a0,v0,v0p,v0n,vt,w0,piv,gap,gan;
  307 + int m,n,k,kz;
  308 +
  309 + a0 = abs(z);
  310 + z1 = z;
  311 + z2 = z*z;
  312 + n = (int)v;
  313 + v0 = v-n;
  314 + piv = M_PI*v0;
  315 + vt = 4.0*v0*v0;
  316 + if (n == 0) n = 1;
  317 + if (a0 < 1e-100) {
  318 + for (k=0;k<=n;k++) {
  319 + civ[k] = czero;
  320 + ckv[k] = complex<double>(-1e308,0);
  321 + civp[k] = czero;
  322 + ckvp[k] = complex<double>(1e308,0);
  323 + }
  324 + if (v0 == 0.0) {
  325 + civ[0] = cone;
  326 + civp[1] = complex<double> (0.5,0.0);
  327 + }
  328 + vm = v;
  329 + return 0;
  330 + }
  331 + if (a0 >= 50.0) kz = 8;
  332 + else if (a0 >= 35.0) kz = 10;
  333 + else kz = 14;
  334 + if (real(z) <= 0.0) z1 = -z;
  335 + if (a0 < 18.0) {
  336 + if (v0 == 0.0) {
  337 + ca1 = cone;
  338 + }
  339 + else {
  340 + v0p = 1.0+v0;
  341 + gap = gamma(v0p);
  342 + ca1 = pow(0.5*z1,v0)/gap;
  343 + }
  344 + ci0 = cone;
  345 + cr = cone;
  346 + for (k=1;k<=50;k++) {
  347 + cr *= 0.25*z2/(k*(k+v0));
  348 + ci0 += cr;
  349 + if (abs(cr/ci0) < eps) break;
  350 + }
  351 + cbi0 = ci0*ca1;
  352 + }
  353 + else {
  354 + ca = exp(z1)/sqrt(2.0*M_PI*z1);
  355 + cs = cone;
  356 + cr = cone;
  357 + for (k=1;k<=kz;k++) {
  358 + cr *= -0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
  359 + cs += cr;
  360 + }
  361 + cbi0 = ca*cs;
  362 + }
  363 + m = msta1(a0,200);
  364 + if (m < n) n = m;
  365 + else m = msta2(a0,n,15);
  366 + cf2 = czero;
  367 + cf1 = complex<double>(1.0e-100,0.0);
  368 + for (k=m;k>=0;k--) {
  369 + cf = 2.0*(v0+k+1.0)*cf1/z1+cf2;
  370 + if (k <= n) civ[k] = cf;
  371 + cf2 = cf1;
  372 + cf1 = cf;
  373 + }
  374 + cs = cbi0/cf;
  375 + for (k=0;k<=n;k++) {
  376 + civ[k] *= cs;
  377 + }
  378 + if (a0 <= 9.0) {
  379 + if (v0 == 0.0) {
  380 + ct = -log(0.5*z1)-el;
  381 + cs = czero;
  382 + w0 = 0.0;
  383 + cr = cone;
  384 + for (k=1;k<=50;k++) {
  385 + w0 += 1.0/k;
  386 + cr *= 0.25*z2/(double)(k*k);
  387 + cp = cr*(w0+ct);
  388 + cs += cp;
  389 + if ((k >= 10) && (abs(cp/cs) < eps)) break;
  390 + }
  391 + cbk0 = ct+cs;
  392 + }
  393 + else {
  394 + v0n = 1.0-v0;
  395 + gan = gamma(v0n);
  396 + ca2 = 1.0/(gan*pow(0.5*z1,v0));
  397 + ca1 = pow(0.5*z1,v0)/gap;
  398 + csu = ca2-ca1;
  399 + cr1 = cone;
  400 + cr2 = cone;
  401 + cws = czero;
  402 + for (k=1;k<=50;k++) {
  403 + cr1 *= 0.25*z2/(k*(k-v0));
  404 + cr2 *= 0.25*z2/(k*(k+v0));
  405 + csu += ca2*cr1-ca1*cr2;
  406 + if ((k >= 10) && (abs((cws-csu)/csu) < eps)) break;
  407 + cws = csu;
  408 + }
  409 + cbk0 = csu*M_PI_2/sin(piv);
  410 + }
  411 + }
  412 + else {
  413 + cb = exp(-z1)*sqrt(M_PI_2/z1);
  414 + cs = cone;
  415 + cr = cone;
  416 + for (k=1;k<=kz;k++) {
  417 + cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
  418 + cs += cr;
  419 + }
  420 + cbk0 = cb*cs;
  421 + }
  422 + cbk1 = (1.0/z1-civ[1]*cbk0)/civ[0];
  423 + ckv[0] = cbk0;
  424 + ckv[1] = cbk1;
  425 + cg0 = cbk0;
  426 + cg1 = cbk1;
  427 + for (k=2;k<=n;k++) {
  428 + cgk = 2.0*(v0+k-1.0)*cg1/z1+cg0;
  429 + ckv[k] = cgk;
  430 + cg0 = cg1;
  431 + cg1 = cgk;
  432 + }
  433 + if (real(z) < 0.0) {
  434 + for (k=0;k<=n;k++) {
  435 + cvk = exp((k+v0)*M_PI*cii);
  436 + if (imag(z) < 0.0) {
  437 + ckv[k] = cvk*ckv[k]+M_PI*cii*civ[k];
  438 + civ[k] /= cvk;
  439 + }
  440 + else if (imag(z) > 0.0) {
  441 + ckv[k] = ckv[k]/cvk-M_PI*cii*civ[k];
  442 + civ[k] *= cvk;
  443 + }
  444 + }
  445 + }
  446 + civp[0] = v0*civ[0]/z+civ[1];
  447 + ckvp[0] = v0*ckv[0]/z-ckv[1];
  448 + for (k=1;k<=n;k++) {
  449 + civp[k] = -(k+v0)*civ[k]/z+civ[k-1];
  450 + ckvp[k] = -(k+v0)*ckv[k]/z-ckv[k-1];
  451 + }
  452 + vm = n+v0;
  453 + return 0;
  454 +}
... ...
CBESSJY.CPP 0 → 100644
  Show/Hide comments   View file @da3d4e0
  1 +++ a/CBESSJY.CPP
  1 +// cbessjy.cpp -- complex Bessel functions.
  2 +// Algorithms and coefficient values from "Computation of Special
  3 +// Functions", Zhang and Jin, John Wiley and Sons, 1996.
  4 +//
  5 +// (C) 2003, C. Bond. All rights reserved.
  6 +//
  7 +#include <complex>
  8 +using namespace std;
  9 +#include "bessel.h"
  10 +double gamma(double);
  11 +
  12 +static complex<double> cii(0.0,1.0);
  13 +static complex<double> cone(1.0,0.0);
  14 +static complex<double> czero(0.0,0.0);
  15 +
  16 +int cbessjy01(complex<double> z,complex<double> &cj0,complex<double> &cj1,
  17 + complex<double> &cy0,complex<double> &cy1,complex<double> &cj0p,
  18 + complex<double> &cj1p,complex<double> &cy0p,complex<double> &cy1p)
  19 +{
  20 + complex<double> z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
  21 + double a0,w0,w1;
  22 + int k,kz;
  23 +
  24 + static double a[] = {
  25 + -7.03125e-2,
  26 + 0.112152099609375,
  27 + -0.5725014209747314,
  28 + 6.074042001273483,
  29 + -1.100171402692467e2,
  30 + 3.038090510922384e3,
  31 + -1.188384262567832e5,
  32 + 6.252951493434797e6,
  33 + -4.259392165047669e8,
  34 + 3.646840080706556e10,
  35 + -3.833534661393944e12,
  36 + 4.854014686852901e14,
  37 + -7.286857349377656e16,
  38 + 1.279721941975975e19};
  39 + static double b[] = {
  40 + 7.32421875e-2,
  41 + -0.2271080017089844,
  42 + 1.727727502584457,
  43 + -2.438052969955606e1,
  44 + 5.513358961220206e2,
  45 + -1.825775547429318e4,
  46 + 8.328593040162893e5,
  47 + -5.006958953198893e7,
  48 + 3.836255180230433e9,
  49 + -3.649010818849833e11,
  50 + 4.218971570284096e13,
  51 + -5.827244631566907e15,
  52 + 9.476288099260110e17,
  53 + -1.792162323051699e20};
  54 + static double a1[] = {
  55 + 0.1171875,
  56 + -0.1441955566406250,
  57 + 0.6765925884246826,
  58 + -6.883914268109947,
  59 + 1.215978918765359e2,
  60 + -3.302272294480852e3,
  61 + 1.276412726461746e5,
  62 + -6.656367718817688e6,
  63 + 4.502786003050393e8,
  64