LifeV
cbessik.cpp
Go to the documentation of this file.
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 #include <lifev/navier_stokes/function/bessel/bessel.hpp>
9 
10 namespace bessel
11 {
12 
13 static std::complex<double> cii (0.0, 1.0);
14 static std::complex<double> czero (0.0, 0.0);
15 static std::complex<double> cone (1.0, 0.0);
16 
17 double gamma (double x);
18 
19 int cbessik01 (std::complex<double>z, std::complex<double>& ci0, std::complex<double>& ci1,
20  std::complex<double>& ck0, std::complex<double>& ck1, std::complex<double>& ci0p,
21  std::complex<double>& ci1p, std::complex<double>& ck0p, std::complex<double>& ck1p)
22 {
23  std::complex<double> z1, z2, zr, zr2, cr, ca, cb, cs, ct, cw;
24  double a0, w0;
25  int k, kz;
26  static double a[] =
27  {
28  0.125,
29  7.03125e-2,
30  7.32421875e-2,
31  1.1215209960938e-1,
32  2.2710800170898e-1,
33  5.7250142097473e-1,
34  1.7277275025845,
35  6.0740420012735,
36  2.4380529699556e1,
37  1.1001714026925e2,
38  5.5133589612202e2,
39  3.0380905109224e3
40  };
41  static double b[] =
42  {
43  -0.375,
44  -1.171875e-1,
45  -1.025390625e-1,
46  -1.4419555664063e-1,
47  -2.7757644653320e-1,
48  -6.7659258842468e-1,
49  -1.9935317337513,
50  -6.8839142681099,
51  -2.7248827311269e1,
52  -1.2159789187654e2,
53  -6.0384407670507e2,
54  -3.3022722944809e3
55  };
56  static double a1[] =
57  {
58  0.125,
59  0.2109375,
60  1.0986328125,
61  1.1775970458984e1,
62  2.1461706161499e2,
63  5.9511522710323e3,
64  2.3347645606175e5,
65  1.2312234987631e7,
66  8.401390346421e08,
67  7.2031420482627e10
68  };
69 
70  a0 = std::abs (z);
71  z2 = z * z;
72  z1 = z;
73  if (a0 == 0.0)
74  {
75  ci0 = cone;
76  ci1 = czero;
77  ck0 = std::complex<double> (1e308, 0);
78  ck1 = std::complex<double> (1e308, 0);
79  ci0p = czero;
80  ci1p = std::complex<double> (0.5, 0.0);
81  ck0p = std::complex<double> (-1e308, 0);
82  ck1p = std::complex<double> (-1e308, 0);
83  return 0;
84  }
85  if (real (z) < 0.0)
86  {
87  z1 = -z;
88  }
89  if (a0 <= 18.0)
90  {
91  ci0 = cone;
92  cr = cone;
93  for (k = 1; k <= 50; k++)
94  {
95  cr *= 0.25 * z2 / (double) (k * k);
96  ci0 += cr;
97  if (std::abs (cr / ci0) < eps)
98  {
99  break;
100  }
101  }
102  ci1 = cone;
103  cr = cone;
104  for (k = 1; k <= 50; k++)
105  {
106  cr *= 0.25 * z2 / (double) (k * (k + 1.0) );
107  ci1 += cr;
108  if (std::abs (cr / ci1) < eps)
109  {
110  break;
111  }
112  }
113  ci1 *= 0.5 * z1;
114  }
115  else
116  {
117  if (a0 >= 50.0)
118  {
119  kz = 7;
120  }
121  else if (a0 >= 35.0)
122  {
123  kz = 9;
124  }
125  else
126  {
127  kz = 12;
128  }
129  ca = std::exp (z1) / std::sqrt (2.0 * M_PI * z1);
130  ci0 = cone;
131  zr = 1.0 / z1;
132  for (k = 0; k < kz; k++)
133  {
134  ci0 += a[k] * std::pow (zr, k + 1.0);
135  }
136  ci0 *= ca;
137  ci1 = cone;
138  for (k = 0; k < kz; k++)
139  {
140  ci1 += b[k] * std::pow (zr, k + 1.0);
141  }
142  ci1 *= ca;
143  }
144  if (a0 <= 9.0)
145  {
146  cs = czero;
147  ct = -log (0.5 * z1) - el;
148  w0 = 0.0;
149  cr = cone;
150  for (k = 1; k <= 50; k++)
151  {
152  w0 += 1.0 / k;
153  cr *= 0.25 * z2 / (double) (k * k);
154  cs += cr * (w0 + ct);
155  if (std::abs ( (cs - cw) / cs) < eps)
156  {
157  break;
158  }
159  cw = cs;
160  }
161  ck0 = ct + cs;
162  }
163  else
164  {
165  cb = 0.5 / z1;
166  zr2 = 1.0 / z2;
167  ck0 = cone;
168  for (k = 0; k < 10; k++)
169  {
170  ck0 += a1[k] * std::pow (zr2, k + 1.0);
171  }
172  ck0 *= cb / ci0;
173  }
174  ck1 = (1.0 / z1 - ci1 * ck0) / ci0;
175  if (real (z) < 0.0)
176  {
177  if (imag (z) < 0.0)
178  {
179  ck0 += cii * M_PI * ci0;
180  ck1 = -ck1 + cii * M_PI * ci1;
181  }
182  else if (imag (z) > 0.0)
183  {
184  ck0 -= cii * M_PI * ci0;
185  ck1 = -ck1 - cii * M_PI * ci1;
186  }
187  ci1 = -ci1;
188  }
189  ci0p = ci1;
190  ci1p = ci0 - 1.0 * ci1 / z;
191  ck0p = -ck1;
192  ck1p = -ck0 - 1.0 * ck1 / z;
193  return 0;
194 }
195 int cbessikna (int n, std::complex<double> z, int& nm, std::complex<double>* ci,
196  std::complex<double>* ck, std::complex<double>* cip, std::complex<double>* ckp)
197 {
198  std::complex<double> ci0, ci1, ck0, ck1, ckk, cf, cf1, cf2, cs;
199  double a0;
200  int k, m;
201  a0 = std::abs (z);
202  nm = n;
203  if (a0 < 1.0e-100)
204  {
205  for (k = 0; k <= n; k++)
206  {
207  ci[k] = czero;
208  ck[k] = std::complex<double> (-1e308, 0);
209  cip[k] = czero;
210  ckp[k] = std::complex<double> (1e308, 0);
211  }
212  ci[0] = cone;
213  cip[1] = std::complex<double> (0.5, 0.0);
214  return 0;
215  }
216  cbessik01 (z, ci[0], ci[1], ck[0], ck[1], cip[0], cip[1], ckp[0], ckp[1]);
217  if (n < 2)
218  {
219  return 0;
220  }
221  ci0 = ci[0];
222  ci1 = ci[1];
223  ck0 = ck[0];
224  ck1 = ck[1];
225  m = msta1 (a0, 200);
226  if (m < n)
227  {
228  nm = m;
229  }
230  else
231  {
232  m = msta2 (a0, n, 15);
233  }
234  cf2 = czero;
235  cf1 = std::complex<double> (1.0e-100, 0.0);
236  for (k = m; k >= 0; k--)
237  {
238  cf = 2.0 * (k + 1.0) * cf1 / z + cf2;
239  if (k <= nm)
240  {
241  ci[k] = cf;
242  }
243  cf2 = cf1;
244  cf1 = cf;
245  }
246  cs = ci0 / cf;
247  for (k = 0; k <= nm; k++)
248  {
249  ci[k] *= cs;
250  }
251  for (k = 2; k <= nm; k++)
252  {
253  if (std::abs (ci[k - 1]) > std::abs (ci[k - 2]) )
254  {
255  ckk = (1.0 / z - ci[k] * ck[k - 1]) / ci[k - 1];
256  }
257  else
258  {
259  ckk = (ci[k] * ck[k - 2] + 2.0 * (k - 1.0) / (z * z) ) / ci[k - 2];
260  }
261  ck[k] = ckk;
262  }
263  for (k = 2; k <= nm; k++)
264  {
265  cip[k] = ci[k - 1] - (double) k * ci[k] / z;
266  ckp[k] = -ck[k - 1] - (double) k * ck[k] / z;
267  }
268  return 0;
269 }
270 int cbessiknb (int n, std::complex<double> z, int& nm, std::complex<double>* ci,
271  std::complex<double>* ck, std::complex<double>* cip, std::complex<double>* ckp)
272 {
273  std::complex<double> z1, cbs, csk0, cf, cf0, cf1, ca0, cbkl;
274  std::complex<double> cg, cg0, cg1, cs0, cs, cr;
275  double a0, vt, fac;
276  int k, kz, l, m;
277 
278  a0 = std::abs (z);
279  nm = n;
280  if (a0 < 1.0e-100)
281  {
282  for (k = 0; k <= n; k++)
283  {
284  ci[k] = czero;
285  ck[k] = std::complex<double> (1e308, 0);
286  cip[k] = czero;
287  ckp[k] = std::complex<double> (-1e308, 0);
288  }
289  ci[0] = std::complex<double> (1.0, 0.0);
290  cip[1] = std::complex<double> (0.5, 0.0);
291  return 0;
292  }
293  z1 = z;
294  if (real (z) < 0.0)
295  {
296  z1 = -z;
297  }
298  if (n == 0)
299  {
300  nm = 1;
301  }
302  m = msta1 (a0, 200);
303  if (m < nm)
304  {
305  nm = m;
306  }
307  else
308  {
309  m = msta2 (a0, nm, 15);
310  }
311  cbs = czero;
312  csk0 = czero;
313  cf0 = czero;
314  cf1 = std::complex<double> (1.0e-100, 0.0);
315  for (k = m; k >= 0; k--)
316  {
317  cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
318  if (k <= nm)
319  {
320  ci[k] = cf;
321  }
322  if ( (k != 0) && (k == 2 * (k >> 1) ) )
323  {
324  csk0 += 4.0 * cf / (double) k;
325  }
326  cbs += 2.0 * cf;
327  cf0 = cf1;
328  cf1 = cf;
329  }
330  cs0 = std::exp (z1) / (cbs - cf);
331  for (k = 0; k <= nm; k++)
332  {
333  ci[k] *= cs0;
334  }
335  if (a0 <= 9.0)
336  {
337  ck[0] = - (log (0.5 * z1) + el) * ci[0] + cs0 * csk0;
338  ck[1] = (1.0 / z1 - ci[1] * ck[0]) / ci[0];
339  }
340  else
341  {
342  ca0 = std::sqrt (M_PI_2 / z1) * std::exp (-z1);
343  if (a0 >= 200.0)
344  {
345  kz = 6;
346  }
347  else if (a0 >= 80.0)
348  {
349  kz = 8;
350  }
351  else if (a0 >= 25.0)
352  {
353  kz = 10;
354  }
355  else
356  {
357  kz = 16;
358  }
359  for (l = 0; l < 2; l++)
360  {
361  cbkl = cone;
362  vt = 4.0 * l;
363  cr = cone;
364  for (k = 1; k <= kz; k++)
365  {
366  cr *= 0.125 * (vt - std::pow (2.0 * k - 1.0, 2.0) ) / ( (double) k * z1);
367  cbkl += cr;
368  }
369  ck[l] = ca0 * cbkl;
370  }
371  }
372  cg0 = ck[0];
373  cg1 = ck[1];
374  for (k = 2; k <= nm; k++)
375  {
376  cg = 2.0 * (k - 1.0) * cg1 / z1 + cg0;
377  ck[k] = cg;
378  cg0 = cg1;
379  cg1 = cg;
380  }
381  if (real (z) < 0.0)
382  {
383  fac = 1.0;
384  for (k = 0; k <= nm; k++)
385  {
386  if (imag (z) < 0.0)
387  {
388  ck[k] = fac * ck[k] + cii * M_PI * ci[k];
389  }
390  else
391  {
392  ck[k] = fac * ck[k] - cii * M_PI * ci[k];
393  }
394  ci[k] *= fac;
395  fac = -fac;
396  }
397  }
398  cip[0] = ci[1];
399  ckp[0] = -ck[1];
400  for (k = 1; k <= nm; k++)
401  {
402  cip[k] = ci[k - 1] - (double) k * ci[k] / z;
403  ckp[k] = -ck[k - 1] - (double) k * ck[k] / z;
404  }
405  return 0;
406 }
407 int cbessikv (double v, std::complex<double>z, double& vm, std::complex<double>* civ,
408  std::complex<double>* ckv, std::complex<double>* civp, std::complex<double>* ckvp)
409 {
410  std::complex<double> z1, z2, ca1, ca, cs, cr, ci0, cbi0, cf, cf1, cf2;
411  std::complex<double> ct, cp, cbk0, ca2, cr1, cr2, csu, cws, cb;
412  std::complex<double> cg0, cg1, cgk, cbk1, cvk;
413  double a0, v0, v0p, v0n, vt, w0, piv, gap (0), gan;
414  int m, n, k, kz;
415 
416  a0 = std::abs (z);
417  z1 = z;
418  z2 = z * z;
419  n = (int) v;
420  v0 = v - n;
421  piv = M_PI * v0;
422  vt = 4.0 * v0 * v0;
423  if (n == 0)
424  {
425  n = 1;
426  }
427  if (a0 < 1e-100)
428  {
429  for (k = 0; k <= n; k++)
430  {
431  civ[k] = czero;
432  ckv[k] = std::complex<double> (-1e308, 0);
433  civp[k] = czero;
434  ckvp[k] = std::complex<double> (1e308, 0);
435  }
436  if (v0 == 0.0)
437  {
438  civ[0] = cone;
439  civp[1] = std::complex<double> (0.5, 0.0);
440  }
441  vm = v;
442  return 0;
443  }
444  if (a0 >= 50.0)
445  {
446  kz = 8;
447  }
448  else if (a0 >= 35.0)
449  {
450  kz = 10;
451  }
452  else
453  {
454  kz = 14;
455  }
456  if (real (z) <= 0.0)
457  {
458  z1 = -z;
459  }
460  if (a0 < 18.0)
461  {
462  if (v0 == 0.0)
463  {
464  ca1 = cone;
465  }
466  else
467  {
468  v0p = 1.0 + v0;
469  gap = gamma (v0p);
470  ca1 = std::pow (0.5 * z1, v0) / gap;
471  }
472  ci0 = cone;
473  cr = cone;
474  for (k = 1; k <= 50; k++)
475  {
476  cr *= 0.25 * z2 / (k * (k + v0) );
477  ci0 += cr;
478  if (std::abs (cr / ci0) < eps)
479  {
480  break;
481  }
482  }
483  cbi0 = ci0 * ca1;
484  }
485  else
486  {
487  ca = std::exp (z1) / std::sqrt (2.0 * M_PI * z1);
488  cs = cone;
489  cr = cone;
490  for (k = 1; k <= kz; k++)
491  {
492  cr *= -0.125 * (vt - std::pow (2.0 * k - 1.0, 2.0) ) / ( (double) k * z1);
493  cs += cr;
494  }
495  cbi0 = ca * cs;
496  }
497  m = msta1 (a0, 200);
498  if (m < n)
499  {
500  n = m;
501  }
502  else
503  {
504  m = msta2 (a0, n, 15);
505  }
506  cf2 = czero;
507  cf1 = std::complex<double> (1.0e-100, 0.0);
508  for (k = m; k >= 0; k--)
509  {
510  cf = 2.0 * (v0 + k + 1.0) * cf1 / z1 + cf2;
511  if (k <= n)
512  {
513  civ[k] = cf;
514  }
515  cf2 = cf1;
516  cf1 = cf;
517  }
518  cs = cbi0 / cf;
519  for (k = 0; k <= n; k++)
520  {
521  civ[k] *= cs;
522  }
523  if (a0 <= 9.0)
524  {
525  if (v0 == 0.0)
526  {
527  ct = -log (0.5 * z1) - el;
528  cs = czero;
529  w0 = 0.0;
530  cr = cone;
531  for (k = 1; k <= 50; k++)
532  {
533  w0 += 1.0 / k;
534  cr *= 0.25 * z2 / (double) (k * k);
535  cp = cr * (w0 + ct);
536  cs += cp;
537  if ( (k >= 10) && (std::abs (cp / cs) < eps) )
538  {
539  break;
540  }
541  }
542  cbk0 = ct + cs;
543  }
544  else
545  {
546  v0n = 1.0 - v0;
547  gan = gamma (v0n);
548  ca2 = 1.0 / (gan * std::pow (0.5 * z1, v0) );
549  ca1 = std::pow (0.5 * z1, v0) / gap;
550  csu = ca2 - ca1;
551  cr1 = cone;
552  cr2 = cone;
553  cws = czero;
554  for (k = 1; k <= 50; k++)
555  {
556  cr1 *= 0.25 * z2 / (k * (k - v0) );
557  cr2 *= 0.25 * z2 / (k * (k + v0) );
558  csu += ca2 * cr1 - ca1 * cr2;
559  if ( (k >= 10) && (std::abs ( (cws - csu) / csu) < eps) )
560  {
561  break;
562  }
563  cws = csu;
564  }
565  cbk0 = csu * M_PI_2 / std::sin (piv);
566  }
567  }
568  else
569  {
570  cb = std::exp (-z1) * std::sqrt (M_PI_2 / z1);
571  cs = cone;
572  cr = cone;
573  for (k = 1; k <= kz; k++)
574  {
575  cr *= 0.125 * (vt - std::pow (2.0 * k - 1.0, 2.0) ) / ( (double) k * z1);
576  cs += cr;
577  }
578  cbk0 = cb * cs;
579  }
580  cbk1 = (1.0 / z1 - civ[1] * cbk0) / civ[0];
581  ckv[0] = cbk0;
582  ckv[1] = cbk1;
583  cg0 = cbk0;
584  cg1 = cbk1;
585  for (k = 2; k <= n; k++)
586  {
587  cgk = 2.0 * (v0 + k - 1.0) * cg1 / z1 + cg0;
588  ckv[k] = cgk;
589  cg0 = cg1;
590  cg1 = cgk;
591  }
592  if (real (z) < 0.0)
593  {
594  for (k = 0; k <= n; k++)
595  {
596  cvk = std::exp ( (k + v0) * M_PI * cii);
597  if (imag (z) < 0.0)
598  {
599  ckv[k] = cvk * ckv[k] + M_PI * cii * civ[k];
600  civ[k] /= cvk;
601  }
602  else if (imag (z) > 0.0)
603  {
604  ckv[k] = ckv[k] / cvk - M_PI * cii * civ[k];
605  civ[k] *= cvk;
606  }
607  }
608  }
609  civp[0] = v0 * civ[0] / z + civ[1];
610  ckvp[0] = v0 * ckv[0] / z - ckv[1];
611  for (k = 1; k <= n; k++)
612  {
613  civp[k] = - (k + v0) * civ[k] / z + civ[k - 1];
614  ckvp[k] = - (k + v0) * ckv[k] / z - ckv[k - 1];
615  }
616  vm = n + v0;
617  return 0;
618 }
619 }
eps
#define eps
Definition: bessel.hpp:14
bessel::msta1
int msta1(double x, int mp)
Definition: bessjy.cpp:312
bessel::cii
static std::complex< double > cii(0.0, 1.0)
bessel::cbessikna
int cbessikna(int n, std::complex< double > z, int &nm, std::complex< double > *ci, std::complex< double > *ck, std::complex< double > *cip, std::complex< double > *ckp)
Definition: cbessik.cpp:195
bessel::msta2
int msta2(double x, int n, int mp)
Definition: bessjy.cpp:337
bessel
Definition: bessik.cpp:9
M_PI
#define M_PI
Definition: winmath.h:20
bessel::czero
static std::complex< double > czero(0.0, 0.0)
bessel::gamma
double gamma(double x)
Definition: gamma.cpp:15
bessel::cbessiknb
int cbessiknb(int n, std::complex< double > z, int &nm, std::complex< double > *ci, std::complex< double > *ck, std::complex< double > *cip, std::complex< double > *ckp)
Definition: cbessik.cpp:270
bessel::cbessikv
int cbessikv(double v, std::complex< double >z, double &vm, std::complex< double > *civ, std::complex< double > *ckv, std::complex< double > *civp, std::complex< double > *ckvp)
Definition: cbessik.cpp:407
bessel::cbessik01
int cbessik01(std::complex< double >z, std::complex< double > &ci0, std::complex< double > &ci1, std::complex< double > &ck0, std::complex< double > &ck1, std::complex< double > &ci0p, std::complex< double > &ci1p, std::complex< double > &ck0p, std::complex< double > &ck1p)
Definition: cbessik.cpp:19
el
#define el
Definition: bessel.hpp:15
bessel::cone
static std::complex< double > cone(1.0, 0.0)