LifeV
bessjy.cpp
Go to the documentation of this file.
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 #include <lifev/navier_stokes/function/bessel/bessel.hpp>
15 namespace bessel
16 {
17 double gamma (double x);
18 //
19 // INPUT:
20 // double x -- argument of Bessel function
21 //
22 // OUTPUT (via address pointers):
23 // double j0 -- Bessel function of 1st kind, 0th order
24 // double j1 -- Bessel function of 1st kind, 1st order
25 // double y0 -- Bessel function of 2nd kind, 0th order
26 // double y1 -- Bessel function of 2nd kind, 1st order
27 // double j0p -- derivative of Bessel function of 1st kind, 0th order
28 // double j1p -- derivative of Bessel function of 1st kind, 1st order
29 // double y0p -- derivative of Bessel function of 2nd kind, 0th order
30 // double y1p -- derivative of Bessel function of 2nd kind, 1st order
31 //
32 // RETURN:
33 // int error code: 0 = OK, 1 = error
34 //
35 // This algorithm computes the above functions using series expansions.
36 //
37 int bessjy01a (double x, double& j0, double& j1, double& y0, double& y1,
38  double& j0p, double& j1p, double& y0p, double& y1p)
39 {
40  double x2, r, ec, w0, w1, r0, r1, cs0, cs1;
41  double cu, p0, q0, p1, q1, t1, t2;
42  int k, kz;
43  static double a[] =
44  {
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  };
60  static double b[] =
61  {
62  7.32421875e-2,
63  -0.2271080017089844,
64  1.727727502584457,
65  -2.438052969955606e1,
66  5.513358961220206e2,
67  -1.825775547429318e4,
68  8.328593040162893e5,
69  -5.006958953198893e7,
70  3.836255180230433e9,
71  -3.649010818849833e11,
72  4.218971570284096e13,
73  -5.827244631566907e15,
74  9.476288099260110e17,
75  -1.792162323051699e20
76  };
77  static double a1[] =
78  {
79  0.1171875,
80  -0.1441955566406250,
81  0.6765925884246826,
82  -6.883914268109947,
83  1.215978918765359e2,
84  -3.302272294480852e3,
85  1.276412726461746e5,
86  -6.656367718817688e6,
87  4.502786003050393e8,
88  -3.833857520742790e10,
89  4.011838599133198e12,
90  -5.060568503314727e14,
91  7.572616461117958e16,
92  -1.326257285320556e19
93  };
94  static double b1[] =
95  {
96  -0.1025390625,
97  0.2775764465332031,
98  -1.993531733751297,
99  2.724882731126854e1,
100  -6.038440767050702e2,
101  1.971837591223663e4,
102  -8.902978767070678e5,
103  5.310411010968522e7,
104  -4.043620325107754e9,
105  3.827011346598605e11,
106  -4.406481417852278e13,
107  6.065091351222699e15,
108  -9.833883876590679e17,
109  1.855045211579828e20
110  };
111 
112  if (x < 0.0)
113  {
114  return 1;
115  }
116  if (x == 0.0)
117  {
118  j0 = 1.0;
119  j1 = 0.0;
120  y0 = -1e308;
121  y1 = -1e308;
122  j0p = 0.0;
123  j1p = 0.5;
124  y0p = 1e308;
125  y1p = 1e308;
126  return 0;
127  }
128  x2 = x * x;
129  if (x <= 12.0)
130  {
131  j0 = 1.0;
132  r = 1.0;
133  for (k = 1; k <= 30; k++)
134  {
135  r *= -0.25 * x2 / (k * k);
136  j0 += r;
137  if (std::fabs (r) < std::fabs (j0) * 1e-15)
138  {
139  break;
140  }
141  }
142  j1 = 1.0;
143  r = 1.0;
144  for (k = 1; k <= 30; k++)
145  {
146  r *= -0.25 * x2 / (k * (k + 1) );
147  j1 += r;
148  if (std::fabs (r) < std::fabs (j1) * 1e-15)
149  {
150  break;
151  }
152  }
153  j1 *= 0.5 * x;
154  ec = std::log (0.5 * x) + el;
155  cs0 = 0.0;
156  w0 = 0.0;
157  r0 = 1.0;
158  for (k = 1; k <= 30; k++)
159  {
160  w0 += 1.0 / k;
161  r0 *= -0.25 * x2 / (k * k);
162  r = r0 * w0;
163  cs0 += r;
164  if (std::fabs (r) < std::fabs (cs0) * 1e-15)
165  {
166  break;
167  }
168  }
169  y0 = M_2_PI * (ec * j0 - cs0);
170  cs1 = 1.0;
171  w1 = 0.0;
172  r1 = 1.0;
173  for (k = 1; k <= 30; k++)
174  {
175  w1 += 1.0 / k;
176  r1 *= -0.25 * x2 / (k * (k + 1) );
177  r = r1 * (2.0 * w1 + 1.0 / (k + 1) );
178  cs1 += r;
179  if (std::fabs (r) < std::fabs (cs1) * 1e-15)
180  {
181  break;
182  }
183  }
184  y1 = M_2_PI * (ec * j1 - 1.0 / x - 0.25 * x * cs1);
185  }
186  else
187  {
188  if (x >= 50.0)
189  {
190  kz = 8; // Can be changed to 10
191  }
192  else if (x >= 35.0)
193  {
194  kz = 10; // " " 12
195  }
196  else
197  {
198  kz = 12; // " " 14
199  }
200  t1 = x - M_PI_4;
201  p0 = 1.0;
202  q0 = -0.125 / x;
203  for (k = 0; k < kz; k++)
204  {
205  p0 += a[k] * std::pow (x, -2 * k - 2);
206  q0 += b[k] * std::pow (x, -2 * k - 3);
207  }
208  cu = std::sqrt (M_2_PI / x);
209  j0 = cu * (p0 * std::cos (t1) - q0 * std::sin (t1) );
210  y0 = cu * (p0 * std::sin (t1) + q0 * std::cos (t1) );
211  t2 = x - 0.75 * M_PI;
212  p1 = 1.0;
213  q1 = 0.375 / x;
214  for (k = 0; k < kz; k++)
215  {
216  p1 += a1[k] * std::pow (x, -2 * k - 2);
217  q1 += b1[k] * std::pow (x, -2 * k - 3);
218  }
219  j1 = cu * (p1 * std::cos (t2) - q1 * std::sin (t2) );
220  y1 = cu * (p1 * std::sin (t2) + q1 * std::cos (t2) );
221  }
222  j0p = -j1;
223  j1p = j0 - j1 / x;
224  y0p = -y1;
225  y1p = y0 - y1 / x;
226  return 0;
227 }
228 //
229 // INPUT:
230 // double x -- argument of Bessel function
231 //
232 // OUTPUT:
233 // double j0 -- Bessel function of 1st kind, 0th order
234 // double j1 -- Bessel function of 1st kind, 1st order
235 // double y0 -- Bessel function of 2nd kind, 0th order
236 // double y1 -- Bessel function of 2nd kind, 1st order
237 // double j0p -- derivative of Bessel function of 1st kind, 0th order
238 // double j1p -- derivative of Bessel function of 1st kind, 1st order
239 // double y0p -- derivative of Bessel function of 2nd kind, 0th order
240 // double y1p -- derivative of Bessel function of 2nd kind, 1st order
241 //
242 // RETURN:
243 // int error code: 0 = OK, 1 = error
244 //
245 // This algorithm computes the functions using polynomial approximations.
246 //
247 int bessjy01b (double x, double& j0, double& j1, double& y0, double& y1,
248  double& j0p, double& j1p, double& y0p, double& y1p)
249 {
250  double t, t2, dtmp, a0, p0, q0, p1, q1, ta0, ta1;
251  if (x < 0.0)
252  {
253  return 1;
254  }
255  if (x == 0.0)
256  {
257  j0 = 1.0;
258  j1 = 0.0;
259  y0 = -1e308;
260  y1 = -1e308;
261  j0p = 0.0;
262  j1p = 0.5;
263  y0p = 1e308;
264  y1p = 1e308;
265  return 0;
266  }
267  if (x <= 4.0)
268  {
269  t = x / 4.0;
270  t2 = t * t;
271  j0 = ( ( ( ( ( (-0.5014415e-3 * t2 + 0.76771853e-2) * t2 - 0.0709253492) * t2 +
272  0.4443584263) * t2 - 1.7777560599) * t2 + 3.9999973021) * t2
273  - 3.9999998721) * t2 + 1.0;
274  j1 = t * ( ( ( ( ( ( (-0.1289769e-3 * t2 + 0.22069155e-2) * t2 - 0.0236616773) * t2 +
275  0.1777582922) * t2 - 0.8888839649) * t2 + 2.6666660544) * t2 -
276  3.999999971) * t2 + 1.9999999998);
277  dtmp = ( ( ( ( ( ( (-0.567433e-4 * t2 + 0.859977e-3) * t2 - 0.94855882e-2) * t2 +
278  0.0772975809) * t2 - 0.4261737419) * t2 + 1.4216421221) * t2 -
279  2.3498519931) * t2 + 1.0766115157) * t2 + 0.3674669052;
280  y0 = M_2_PI * std::log (0.5 * x) * j0 + dtmp;
281  dtmp = ( ( ( ( ( ( (0.6535773e-3 * t2 - 0.0108175626) * t2 + 0.107657607) * t2 -
282  0.7268945577) * t2 + 3.1261399273) * t2 - 7.3980241381) * t2 +
283  6.8529236342) * t2 + 0.3932562018) * t2 - 0.6366197726;
284  y1 = M_2_PI * std::log (0.5 * x) * j1 + dtmp / x;
285  }
286  else
287  {
288  t = 4.0 / x;
289  t2 = t * t;
290  a0 = std::sqrt (M_2_PI / x);
291  p0 = ( ( ( (-0.9285e-5 * t2 + 0.43506e-4) * t2 - 0.122226e-3) * t2 +
292  0.434725e-3) * t2 - 0.4394275e-2) * t2 + 0.999999997;
293  q0 = t * ( ( ( ( (0.8099e-5 * t2 - 0.35614e-4) * t2 + 0.85844e-4) * t2 -
294  0.218024e-3) * t2 + 0.1144106e-2) * t2 - 0.031249995);
295  ta0 = x - M_PI_4;
296  j0 = a0 * (p0 * std::cos (ta0) - q0 * std::sin (ta0) );
297  y0 = a0 * (p0 * std::sin (ta0) + q0 * std::cos (ta0) );
298  p1 = ( ( ( (0.10632e-4 * t2 - 0.50363e-4) * t2 + 0.145575e-3) * t2
299  - 0.559487e-3) * t2 + 0.7323931e-2) * t2 + 1.000000004;
300  q1 = t * ( ( ( ( (-0.9173e-5 * t2 + 0.40658e-4) * t2 - 0.99941e-4) * t2
301  + 0.266891e-3) * t2 - 0.1601836e-2) * t2 + 0.093749994);
302  ta1 = x - 0.75 * M_PI;
303  j1 = a0 * (p1 * std::cos (ta1) - q1 * std::sin (ta1) );
304  y1 = a0 * (p1 * std::sin (ta1) + q1 * std::cos (ta1) );
305  }
306  j0p = -j1;
307  j1p = j0 - j1 / x;
308  y0p = -y1;
309  y1p = y0 - y1 / x;
310  return 0;
311 }
312 int msta1 (double x, int mp)
313 {
314  double a0, f0, f1, f;
315  int i, n0, n1, nn;
316 
317  a0 = std::fabs (x);
318  n0 = (int) (1.1 * a0) + 1;
319  f0 = 0.5 * std::log10 (6.28 * n0) - n0 * std::log10 (1.36 * a0 / n0) - mp;
320  n1 = n0 + 5;
321  f1 = 0.5 * std::log10 (6.28 * n1) - n1 * std::log10 (1.36 * a0 / n1) - mp;
322  for (i = 0; i < 20; i++)
323  {
324  nn = (int) (n1 - (n1 - n0) / (1.0 - f0 / f1) );
325  f = 0.5 * std::log10 (6.28 * nn) - nn * std::log10 (1.36 * a0 / nn) - mp;
326  if (std::abs (nn - n1) < 1)
327  {
328  break;
329  }
330  n0 = n1;
331  f0 = f1;
332  n1 = nn;
333  f1 = f;
334  }
335  return nn;
336 }
337 int msta2 (double x, int n, int mp)
338 {
339  double a0, ejn, hmp, f0, f1, f, obj;
340  int i, n0, n1, nn;
341 
342  a0 = std::fabs (x);
343  hmp = 0.5 * mp;
344  ejn = 0.5 * std::log10 (6.28 * n) - n * std::log10 (1.36 * a0 / n);
345  if (ejn <= hmp)
346  {
347  obj = mp;
348  n0 = (int) (1.1 * a0);
349  if (n0 < 1)
350  {
351  n0 = 1;
352  }
353  }
354  else
355  {
356  obj = hmp + ejn;
357  n0 = n;
358  }
359  f0 = 0.5 * std::log10 (6.28 * n0) - n0 * std::log10 (1.36 * a0 / n0) - obj;
360  n1 = n0 + 5;
361  f1 = 0.5 * std::log10 (6.28 * n1) - n1 * std::log10 (1.36 * a0 / n1) - obj;
362  for (i = 0; i < 20; i++)
363  {
364  nn = (int) (n1 - (n1 - n0) / (1.0 - f0 / f1) );
365  f = 0.5 * std::log10 (6.28 * nn) - nn * std::log10 (1.36 * a0 / nn) - obj;
366  if (std::abs (nn - n1) < 1)
367  {
368  break;
369  }
370  n0 = n1;
371  f0 = f1;
372  n1 = nn;
373  f1 = f;
374  }
375  return nn + 10;
376 }
377 //
378 // INPUT:
379 // double x -- argument of Bessel function of 1st and 2nd kind.
380 // int n -- order
381 //
382 // OUPUT:
383 //
384 // int nm -- highest order actually computed (nm <= n)
385 // double jn[] -- Bessel function of 1st kind, orders from 0 to nm
386 // double yn[] -- Bessel function of 2nd kind, orders from 0 to nm
387 // double j'n[]-- derivative of Bessel function of 1st kind,
388 // orders from 0 to nm
389 // double y'n[]-- derivative of Bessel function of 2nd kind,
390 // orders from 0 to nm
391 //
392 // Computes Bessel functions of all order up to 'n' using recurrence
393 // relations. If 'nm' < 'n' only 'nm' orders are returned.
394 //
395 int bessjyna (int n, double x, int& nm, double* jn, double* yn,
396  double* jnp, double* ynp)
397 {
398  double bj0, bj1, f (0), f0, f1, f2, cs;
399  int i, k, m;
400 
401  nm = n;
402  if ( (x < 0.0) || (n < 0) )
403  {
404  return 1;
405  }
406  if (x < 1e-15)
407  {
408  for (i = 0; i <= n; i++)
409  {
410  jn[i] = 0.0;
411  yn[i] = -1e308;
412  jnp[i] = 0.0;
413  ynp[i] = 1e308;
414  }
415  jn[0] = 1.0;
416  jnp[1] = 0.5;
417  return 0;
418  }
419  bessjy01a (x, jn[0], jn[1], yn[0], yn[1], jnp[0], jnp[1], ynp[0], ynp[1]);
420  if (n < 2)
421  {
422  return 0;
423  }
424  bj0 = jn[0];
425  bj1 = jn[1];
426  if (n < (int) 0.9 * x)
427  {
428  for (k = 2; k <= n; k++)
429  {
430  jn[k] = 2.0 * (k - 1.0) * bj1 / x - bj0;
431  bj0 = bj1;
432  bj1 = jn[k];
433  }
434  }
435  else
436  {
437  m = msta1 (x, 200);
438  if (m < n)
439  {
440  nm = m;
441  }
442  else
443  {
444  m = msta2 (x, n, 15);
445  }
446  f2 = 0.0;
447  f1 = 1.0e-100;
448  for (k = m; k >= 0; k--)
449  {
450  f = 2.0 * (k + 1.0) / x * f1 - f2;
451  if (k <= nm)
452  {
453  jn[k] = f;
454  }
455  f2 = f1;
456  f1 = f;
457  }
458  if (std::fabs (bj0) > std::fabs (bj1) )
459  {
460  cs = bj0 / f;
461  }
462  else
463  {
464  cs = bj1 / f2;
465  }
466  for (k = 0; k <= nm; k++)
467  {
468  jn[k] *= cs;
469  }
470  }
471  for (k = 2; k <= nm; k++)
472  {
473  jnp[k] = jn[k - 1] - k * jn[k] / x;
474  }
475  f0 = yn[0];
476  f1 = yn[1];
477  for (k = 2; k <= nm; k++)
478  {
479  f = 2.0 * (k - 1.0) * f1 / x - f0;
480  yn[k] = f;
481  f0 = f1;
482  f1 = f;
483  }
484  for (k = 2; k <= nm; k++)
485  {
486  ynp[k] = yn[k - 1] - k * yn[k] / x;
487  }
488  return 0;
489 }
490 //
491 // Same input and output conventions as above. Different recurrence
492 // relations used for 'x' < 300.
493 //
494 int bessjynb (int n, double x, int& nm, double* jn, double* yn,
495  double* jnp, double* ynp)
496 {
497  double t1, t2, f (0), f1, f2, bj0, bj1, bjk, by0, by1, cu, s0, su, sv;
498  double ec, bs, byk, p0, p1, q0, q1;
499  static double a[] =
500  {
501  -0.7031250000000000e-1,
502  0.1121520996093750,
503  -0.5725014209747314,
504  6.074042001273483
505  };
506  static double b[] =
507  {
508  0.7324218750000000e-1,
509  -0.2271080017089844,
510  1.727727502584457,
511  -2.438052969955606e1
512  };
513  static double a1[] =
514  {
515  0.1171875,
516  -0.1441955566406250,
517  0.6765925884246826,
518  -6.883914268109947
519  };
520  static double b1[] =
521  {
522  -0.1025390625,
523  0.2775764465332031,
524  -1.993531733751297,
525  2.724882731126854e1
526  };
527 
528  int i, k, m;
529  nm = n;
530  if ( (x < 0.0) || (n < 0) )
531  {
532  return 1;
533  }
534  if (x < 1e-15)
535  {
536  for (i = 0; i <= n; i++)
537  {
538  jn[i] = 0.0;
539  yn[i] = -1e308;
540  jnp[i] = 0.0;
541  ynp[i] = 1e308;
542  }
543  jn[0] = 1.0;
544  jnp[1] = 0.5;
545  return 0;
546  }
547  if (x <= 300.0 || n > (int) (0.9 * x) )
548  {
549  if (n == 0)
550  {
551  nm = 1;
552  }
553  m = msta1 (x, 200);
554  if (m < nm)
555  {
556  nm = m;
557  }
558  else
559  {
560  m = msta2 (x, nm, 15);
561  }
562  bs = 0.0;
563  su = 0.0;
564  sv = 0.0;
565  f2 = 0.0;
566  f1 = 1.0e-100;
567  for (k = m; k >= 0; k--)
568  {
569  f = 2.0 * (k + 1.0) / x * f1 - f2;
570  if (k <= nm)
571  {
572  jn[k] = f;
573  }
574  if ( (k == 2 * (int) (k / 2) ) && (k != 0) )
575  {
576  bs += 2.0 * f;
577  // su += std::pow(-1,k>>1)*f/(double)k;
578  su += (-1) * ( (k & 2) - 1) * f / (double) k;
579  }
580  else if (k > 1)
581  {
582  // sv += std::pow(-1,k>>1)*k*f/(k*k-1.0);
583  sv += (-1) * ( (k & 2) - 1) * (double) k * f / (k * k - 1.0);
584  }
585  f2 = f1;
586  f1 = f;
587  }
588  s0 = bs + f;
589  for (k = 0; k <= nm; k++)
590  {
591  jn[k] /= s0;
592  }
593  ec = std::log (0.5 * x) + 0.5772156649015329;
594  by0 = M_2_PI * (ec * jn[0] - 4.0 * su / s0);
595  yn[0] = by0;
596  by1 = M_2_PI * ( (ec - 1.0) * jn[1] - jn[0] / x - 4.0 * sv / s0);
597  yn[1] = by1;
598  }
599  else
600  {
601  t1 = x - M_PI_4;
602  p0 = 1.0;
603  q0 = -0.125 / x;
604  for (k = 0; k < 4; k++)
605  {
606  p0 += a[k] * std::pow (x, -2 * k - 2);
607  q0 += b[k] * std::pow (x, -2 * k - 3);
608  }
609  cu = std::sqrt (M_2_PI / x);
610  bj0 = cu * (p0 * std::cos (t1) - q0 * std::sin (t1) );
611  by0 = cu * (p0 * std::sin (t1) + q0 * std::cos (t1) );
612  jn[0] = bj0;
613  yn[0] = by0;
614  t2 = x - 0.75 * M_PI;
615  p1 = 1.0;
616  q1 = 0.375 / x;
617  for (k = 0; k < 4; k++)
618  {
619  p1 += a1[k] * std::pow (x, -2 * k - 2);
620  q1 += b1[k] * std::pow (x, -2 * k - 3);
621  }
622  bj1 = cu * (p1 * std::cos (t2) - q1 * std::sin (t2) );
623  by1 = cu * (p1 * std::sin (t2) + q1 * std::cos (t2) );
624  jn[1] = bj1;
625  yn[1] = by1;
626  for (k = 2; k <= nm; k++)
627  {
628  bjk = 2.0 * (k - 1.0) * bj1 / x - bj0;
629  jn[k] = bjk;
630  bj0 = bj1;
631  bj1 = bjk;
632  }
633  }
634  jnp[0] = -jn[1];
635  for (k = 1; k <= nm; k++)
636  {
637  jnp[k] = jn[k - 1] - k * jn[k] / x;
638  }
639  for (k = 2; k <= nm; k++)
640  {
641  byk = 2.0 * (k - 1.0) * by1 / x - by0;
642  yn[k] = byk;
643  by0 = by1;
644  by1 = byk;
645  }
646  ynp[0] = -yn[1];
647  for (k = 1; k <= nm; k++)
648  {
649  ynp[k] = yn[k - 1] - k * yn[k] / x;
650  }
651  return 0;
652 
653 }
654 
655 // The following routine computes Bessel Jv(x) and Yv(x) for
656 // arbitrary positive order (v). For negative order, use:
657 //
658 // J-v(x) = Jv(x)cos(v pi) - Yv(x)sin(v pi)
659 // Y-v(x) = Jv(x)sin(v pi) + Yv(x)cos(v pi)
660 //
661 int bessjyv (double v, double x, double& vm, double* jv, double* yv,
662  double* djv, double* dyv)
663 {
664  double v0, vl, vg, vv, a, a0, r, x2, bjv0 (0), bjv1, bjvl, f (0), f0, f1, f2;
665  double r0, r1, ck, cs, cs0, cs1, sk, qx, px, byv0 (0), byv1 (0), rp, xk, rq;
666  double b, ec, w0, w1, bju0 (0), bju1, pv0, pv1, byvk;
667  int j, k, l, m, n, kz;
668 
669  x2 = x * x;
670  n = (int) v;
671  v0 = v - n;
672  if ( (x < 0.0) || (v < 0.0) )
673  {
674  return 1;
675  }
676  if (x < 1e-15)
677  {
678  for (k = 0; k <= n; k++)
679  {
680  jv[k] = 0.0;
681  yv[k] = -1e308;
682  djv[k] = 0.0;
683  dyv[k] = 1e308;
684  if (v0 == 0.0)
685  {
686  jv[0] = 1.0;
687  djv[1] = 0.5;
688  }
689  else
690  {
691  djv[0] = 1e308;
692  }
693  }
694  vm = v;
695  return 0;
696  }
697  if (x <= 12.0)
698  {
699  for (l = 0; l < 2; l++)
700  {
701  vl = v0 + l;
702  bjvl = 1.0;
703  r = 1.0;
704  for (k = 1; k <= 40; k++)
705  {
706  r *= -0.25 * x2 / (k * (k + vl) );
707  bjvl += r;
708  if (std::fabs (r) < std::fabs (bjvl) * 1e-15)
709  {
710  break;
711  }
712  }
713  vg = 1.0 + vl;
714  a = std::pow (0.5 * x, vl) / gamma (vg);
715  if (l == 0)
716  {
717  bjv0 = bjvl * a;
718  }
719  else
720  {
721  bjv1 = bjvl * a;
722  }
723  }
724  }
725  else
726  {
727  if (x >= 50.0)
728  {
729  kz = 8;
730  }
731  else if (x >= 35.0)
732  {
733  kz = 10;
734  }
735  else
736  {
737  kz = 11;
738  }
739  for (j = 0; j < 2; j++)
740  {
741  vv = 4.0 * (j + v0) * (j + v0);
742  px = 1.0;
743  rp = 1.0;
744  for (k = 1; k <= kz; k++)
745  {
746  rp *= (-0.78125e-2) * (vv - std::pow (4.0 * k - 3.0, 2.0) ) *
747  (vv - std::pow (4.0 * k - 1.0, 2.0) ) / (k * (2.0 * k - 1.0) * x2);
748  px += rp;
749  }
750  qx = 1.0;
751  rq = 1.0;
752  for (k = 1; k <= kz; k++)
753  {
754  rq *= (-0.78125e-2) * (vv - std::pow (4.0 * k - 1.0, 2.0) ) *
755  (vv - std::pow (4.0 * k + 1.0, 2.0) ) / (k * (2.0 * k + 1.0) * x2);
756  qx += rq;
757  }
758  qx *= 0.125 * (vv - 1.0) / x;
759  xk = x - (0.5 * (j + v0) + 0.25) * M_PI;
760  a0 = std::sqrt (M_2_PI / x);
761  ck = std::cos (xk);
762  sk = std::sin (xk);
763 
764  if (j == 0)
765  {
766  bjv0 = a0 * (px * ck - qx * sk);
767  byv0 = a0 * (px * sk + qx * ck);
768  }
769  else
770  {
771  bjv1 = a0 * (px * ck - qx * sk);
772  byv1 = a0 * (px * sk + qx * ck);
773  }
774  }
775  }
776  jv[0] = bjv0;
777  jv[1] = bjv1;
778  djv[0] = v0 * jv[0] / x - jv[1];
779  djv[1] = - (1.0 + v0) * jv[1] / x + jv[0];
780  if ( (n >= 2) && (n <= (int) (0.9 * x) ) )
781  {
782  f0 = bjv0;
783  f1 = bjv1;
784  for (k = 2; k <= n; k++)
785  {
786  f = 2.0 * (k + v0 - 1.0) * f1 / x - f0;
787  jv[k] = f;
788  f0 = f1;
789  f1 = f;
790  }
791  }
792  else if (n >= 2)
793  {
794  m = msta1 (x, 200);
795  if (m < n)
796  {
797  n = m;
798  }
799  else
800  {
801  m = msta2 (x, n, 15);
802  }
803  f2 = 0.0;
804  f1 = 1.0e-100;
805  for (k = m; k >= 0; k--)
806  {
807  f = 2.0 * (v0 + k + 1.0) / x * f1 - f2;
808  if (k <= n)
809  {
810  jv[k] = f;
811  }
812  f2 = f1;
813  f1 = f;
814  }
815  if (std::fabs (bjv0) > std::fabs (bjv1) )
816  {
817  cs = bjv0 / f;
818  }
819  else
820  {
821  cs = bjv1 / f2;
822  }
823  for (k = 0; k <= n; k++)
824  {
825  jv[k] *= cs;
826  }
827  }
828  for (k = 2; k <= n; k++)
829  {
830  djv[k] = - (k + v0) * jv[k] / x + jv[k - 1];
831  }
832  if (x <= 12.0)
833  {
834  if (v0 != 0.0)
835  {
836  for (l = 0; l < 2; l++)
837  {
838  vl = v0 + l;
839  bjvl = 1.0;
840  r = 1.0;
841  for (k = 1; k <= 40; k++)
842  {
843  r *= -0.25 * x2 / (k * (k - vl) );
844  bjvl += r;
845  if (std::fabs (r) < std::fabs (bjvl) * 1e-15)
846  {
847  break;
848  }
849  }
850  vg = 1.0 - vl;
851  b = std::pow (2.0 / x, vl) / gamma (vg);
852  if (l == 0)
853  {
854  bju0 = bjvl * b;
855  }
856  else
857  {
858  bju1 = bjvl * b;
859  }
860  }
861  pv0 = M_PI * v0;
862  pv1 = M_PI * (1.0 + v0);
863  byv0 = (bjv0 * std::cos (pv0) - bju0) / std::sin (pv0);
864  byv1 = (bjv1 * std::cos (pv1) - bju1) / std::sin (pv1);
865  }
866  else
867  {
868  ec = std::log (0.5 * x) + el;
869  cs0 = 0.0;
870  w0 = 0.0;
871  r0 = 1.0;
872  for (k = 1; k <= 30; k++)
873  {
874  w0 += 1.0 / k;
875  r0 *= -0.25 * x2 / (k * k);
876  cs0 += r0 * w0;
877  }
878  byv0 = M_2_PI * (ec * bjv0 - cs0);
879  cs1 = 1.0;
880  w1 = 0.0;
881  r1 = 1.0;
882  for (k = 1; k <= 30; k++)
883  {
884  w1 += 1.0 / k;
885  r1 *= -0.25 * x2 / (k * (k + 1) );
886  cs1 += r1 * (2.0 * w1 + 1.0 / (k + 1.0) );
887  }
888  byv1 = M_2_PI * (ec * bjv1 - 1.0 / x - 0.25 * x * cs1);
889  }
890  }
891  yv[0] = byv0;
892  yv[1] = byv1;
893  for (k = 2; k <= n; k++)
894  {
895  byvk = 2.0 * (v0 + k - 1.0) * byv1 / x - byv0;
896  yv[k] = byvk;
897  byv0 = byv1;
898  byv1 = byvk;
899  }
900  dyv[0] = v0 * yv[0] / x - yv[1];
901  for (k = 1; k <= n; k++)
902  {
903  dyv[k] = - (k + v0) * yv[k] / x + yv[k - 1];
904  }
905  vm = n + v0;
906  return 0;
907 }
908 }
bessel::msta1
int msta1(double x, int mp)
Definition: bessjy.cpp:312
bessel::msta2
int msta2(double x, int n, int mp)
Definition: bessjy.cpp:337
bessel
Definition: bessik.cpp:9
bessel::bessjynb
int bessjynb(int n, double x, int &nm, double *jn, double *yn, double *jnp, double *ynp)
Definition: bessjy.cpp:494
M_PI
#define M_PI
Definition: winmath.h:20
bessel::gamma
double gamma(double x)
Definition: gamma.cpp:15
bessel::bessjyv
int bessjyv(double v, double x, double &vm, double *jv, double *yv, double *jvp, double *yvp)
Definition: bessjy.cpp:661
bessel::bessjy01b
int bessjy01b(double x, double &j0, double &j1, double &y0, double &y1, double &j0p, double &j1p, double &y0p, double &y1p)
Definition: bessjy.cpp:247
bessel::bessjyna
int bessjyna(int n, double x, int &nm, double *jn, double *yn, double *jnp, double *ynp)
Definition: bessjy.cpp:395
el
#define el
Definition: bessel.hpp:15
bessel::bessjy01a
int bessjy01a(double x, double &j0, double &j1, double &y0, double &y1, double &j0p, double &j1p, double &y0p, double &y1p)
Definition: bessjy.cpp:37