optics / imie

Commit f095146f8a8052705b1edd70ac7380a8a18aea96

Authored by dmayerich
1 parent bfe3f56b

Case-sensitive errors with the bessel functions.

Showing 30 changed files with 10232 additions and 12646 deletions   Show diff stats
Inline Side-by-side
BESSEL.h deleted
View file @bfe3f56
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 deleted
View file @bfe3f56
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,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 deleted
View file @bfe3f56
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 = (int)(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 = (int)(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,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,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 deleted
View file @bfe3f56
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 deleted
View file @bfe3f56
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   - -3.833857520742790e10,
65   - 4.011838599133198e12,
66   - -5.060568503314727e14,
67   - 7.572616461117958e16,
68   - -1.326257285320556e19};
69   - static double b1[] = {
70   - -0.1025390625,
71   - 0.2775764465332031,
72   - -1.993531733751297,
73   - 2.724882731126854e1,
74   - -6.038440767050702e2,
75   - 1.971837591223663e4,
76   - -8.902978767070678e5,
77   - 5.310411010968522e7,