Functions | Variables
CRBond_Bessel Namespace Reference

Functions

double gamma (double x)
 
int bessik01a (double x, double &i0, double &i1, double &k0, double &k1, double &i0p, double &i1p, double &k0p, double &k1p)
 
int bessik01b (double x, double &i0, double &i1, double &k0, double &k1, double &i0p, double &i1p, double &k0p, double &k1p)
 
int bessikna (int n, double x, int &nm, double *in, double *kn, double *inp, double *knp)
 
int bessiknb (int n, double x, int &nm, double *in, double *kn, double *inp, double *knp)
 
int bessikv (double v, double x, double &vm, double *iv, double *kv, double *ivp, double *kvp)
 
int bessjy01a (double x, double &j0, double &j1, double &y0, double &y1, double &j0p, double &j1p, double &y0p, double &y1p)
 
int bessjy01b (double x, double &j0, double &j1, double &y0, double &y1, double &j0p, double &j1p, double &y0p, double &y1p)
 
int msta1 (double x, int mp)
 
int msta2 (double x, int n, int mp)
 
int bessjyna (int n, double x, int &nm, double *jn, double *yn, double *jnp, double *ynp)
 
int bessjynb (int n, double x, int &nm, double *jn, double *yn, double *jnp, double *ynp)
 
int bessjyv (double v, double x, double &vm, double *jv, double *yv, double *djv, double *dyv)
 
static complex< double > cii (0.0, 1.0)
 
static complex< double > czero (0.0, 0.0)
 
static complex< double > cone (1.0, 0.0)
 
int cbessik01 (complex< double >z, complex< double > &ci0, complex< double > &ci1, complex< double > &ck0, complex< double > &ck1, complex< double > &ci0p, complex< double > &ci1p, complex< double > &ck0p, complex< double > &ck1p)
 
int cbessikna (int n, complex< double > z, int &nm, complex< double > *ci, complex< double > *ck, complex< double > *cip, complex< double > *ckp)
 
int cbessiknb (int n, complex< double > z, int &nm, complex< double > *ci, complex< double > *ck, complex< double > *cip, complex< double > *ckp)
 
int cbessikv (double v, complex< double >z, double &vm, complex< double > *civ, complex< double > *ckv, complex< double > *civp, complex< double > *ckvp)
 
int cbessjy01 (complex< double > z, complex< double > &cj0, complex< double > &cj1, complex< double > &cy0, complex< double > &cy1, complex< double > &cj0p, complex< double > &cj1p, complex< double > &cy0p, complex< double > &cy1p)
 
int cbessjyna (int n, complex< double > z, int &nm, complex< double > *cj, complex< double > *cy, complex< double > *cjp, complex< double > *cyp)
 
int cbessjynb (int n, complex< double > z, int &nm, complex< double > *cj, complex< double > *cy, complex< double > *cjp, complex< double > *cyp)
 
int cbessjyva (double v, complex< double > z, double &vm, complex< double > *cjv, complex< double > *cyv, complex< double > *cjvp, complex< double > *cyvp)
 

Variables

double eps =1.0e-15
 
double el =0.5772156649015329
 

Detailed Description

Collection of Bessel functions and the like, implemented in C++ by C. Bond (http://www.crbond.com/), using the algorithms in Zhang and Jin's book "Computation of Special Functions". John Wiley and Sons, 1996. This is based on the version dated 06/04 which appears to have corrected errors in complex Bessel functions.

Function Documentation

◆ bessik01a()

int CRBond_Bessel::bessik01a ( double  x,
double &  i0,
double &  i1,
double &  k0,
double &  k1,
double &  i0p,
double &  i1p,
double &  k0p,
double &  k1p 
)
43 {
44  double r,x2,ca,cb,ct,ww,w0,xr,xr2;
45  int k,kz;
46  static double a[] = {
47  0.125,
48  7.03125e-2,
49  7.32421875e-2,
50  1.1215209960938e-1,
51  2.2710800170898e-1,
52  5.7250142097473e-1,
53  1.7277275025845,
54  6.0740420012735,
55  2.4380529699556e1,
56  1.1001714026925e2,
57  5.5133589612202e2,
58  3.0380905109224e3};
59  static double b[] = {
60  -0.375,
61  -1.171875e-1,
62  -1.025390625e-1,
63  -1.4419555664063e-1,
64  -2.7757644653320e-1,
65  -6.7659258842468e-1,
66  -1.9935317337513,
67  -6.8839142681099,
68  -2.7248827311269e1,
69  -1.2159789187654e2,
70  -6.0384407670507e2,
71  -3.3022722944809e3};
72  static double a1[] = {
73  0.125,
74  0.2109375,
75  1.0986328125,
76  1.1775970458984e1,
77  2.1461706161499e2,
78  5.9511522710323e3,
79  2.3347645606175e5,
80  1.2312234987631e7};
81 
82  if (x < 0.0) return 1;
83  if (x == 0.0) {
84  i0 = 1.0;
85  i1 = 0.0;
86  k0 = 1e308;
87  k1 = 1e308;
88  i0p = 0.0;
89  i1p = 0.5;
90  k0p = -1e308;
91  k1p = -1e308;
92  return 0;
93  }
94  x2 = x*x;
95  if (x <= 18.0) {
96  i0 = 1.0;
97  r = 1.0;
98  for (k=1;k<=50;k++) {
99  r *= 0.25*x2/(k*k);
100  i0 += r;
101  if (fabs(r/i0) < eps) break;
102  }
103  i1 = 1.0;
104  r = 1.0;
105  for (k=1;k<=50;k++) {
106  r *= 0.25*x2/(k*(k+1));
107  i1 += r;
108  if (fabs(r/i1) < eps) break;
109  }
110  i1 *= 0.5*x;
111  }
112  else {
113  if (x >= 50.0) kz = 7;
114  else if (x >= 35.0) kz = 9;
115  else kz = 12;
116  ca = exp(x)/sqrt(2.0*M_PI*x);
117  i0 = 1.0;
118  xr = 1.0/x;
119  for (k=0;k<kz;k++) {
120  i0 += a[k]*pow(xr,k+1);
121  }
122  i0 *= ca;
123  i1 = 1.0;
124  for (k=0;k<kz;k++) {
125  i1 += b[k]*pow(xr,k+1);
126  }
127  i1 *= ca;
128  }
129  if (x <= 9.0) {
130  ct = -(log(0.5*x)+el);
131  k0 = 0.0;
132  w0 = 0.0;
133  r = 1.0;
134  ww = 0.0;
135  for (k=1;k<=50;k++) {
136  w0 += 1.0/k;
137  r *= 0.25*x2/(k*k);
138  k0 += r*(w0+ct);
139  if (fabs((k0-ww)/k0) < eps) break;
140  ww = k0;
141  }
142  k0 += ct;
143  }
144  else {
145  cb = 0.5/x;
146  xr2 = 1.0/x2;
147  k0 = 1.0;
148  for (k=0;k<8;k++) {
149  k0 += a1[k]*pow(xr2,k+1);
150  }
151  k0 *= cb/i0;
152  }
153  k1 = (1.0/x - i1*k0)/i0;
154  i0p = i1;
155  i1p = i0-i1/x;
156  k0p = -k1;
157  k1p = -k0-k1/x;
158  return 0;
159 }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::ArrayBase::pow
EIGEN_DEVICE_FUNC const GlobalUnaryPowReturnType< Derived, ScalarExponent > pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
Definition: GlobalFunctions.h:137
a
const Scalar * a
Definition: level2_cplx_impl.h:32
k
char char char int int * k
Definition: level2_impl.h:374
M_PI
#define M_PI
Definition: main.h:121
CRBond_Bessel::eps
double eps
Definition: crbond_bessel.cc:24
CRBond_Bessel::el
double el
Definition: crbond_bessel.cc:27
Eigen::bfloat16_impl::exp
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 exp(const bfloat16 &a)
Definition: BFloat16.h:615
Eigen::bfloat16_impl::log
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 log(const bfloat16 &a)
Definition: BFloat16.h:618
Global_parameters::x2
Vector< double > x2(const Vector< double > &coord)
Cartesian coordinates centered at the point (1.5,1)
Definition: poisson/poisson_with_singularity/two_d_poisson.cc:102
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
plotDoE.x
list x
Definition: plotDoE.py:28

References a, b, el, eps, Eigen::bfloat16_impl::exp(), boost::multiprecision::fabs(), k, Eigen::bfloat16_impl::log(), M_PI, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, sqrt(), plotDoE::x, and Global_parameters::x2().

Referenced by bessikna().

◆ bessik01b()

int CRBond_Bessel::bessik01b ( double  x,
double &  i0,
double &  i1,
double &  k0,
double &  k1,
double &  i0p,
double &  i1p,
double &  k0p,
double &  k1p 
)
163 {
164  double t,t2,dtmp,dtmp1;
165 
166  if (x < 0.0) return 1;
167  if (x == 0.0) {
168  i0 = 1.0;
169  i1 = 0.0;
170  k0 = 1e308;
171  k1 = 1e308;
172  i0p = 0.0;
173  i1p = 0.5;
174  k0p = -1e308;
175  k1p = -1e308;
176  return 0;
177  }
178  if (x < 3.75) {
179  t = x/3.75;
180  t2 = t*t;
181  i0 = (((((0.0045813*t2+0.0360768)*t2+0.2659732)*t2+
182  1.2067492)*t2+3.0899424)*t2+3.5156229)*t2+1.0;
183  i1 = x*(((((0.00032411*t2+0.00301532)*t2+0.02658733*t2+
184  0.15084934)*t2+0.51498869)*t2+0.87890594)*t2+0.5);
185  }
186  else {
187  t = 3.75/x;
188  dtmp1 = exp(x)/sqrt(x);
189  dtmp = (((((((0.00392377*t-0.01647633)*t+0.026355537)*t-0.02057706)*t+
190  0.00916281)*t-0.00157565)*t+0.00225319)*t+0.01328592)*t+0.39894228;
191  i0 = dtmp*dtmp1;
192  dtmp = (((((((-0.00420059*t+0.01787654)*t-0.02895312)*t+0.02282967)*t-
193  0.01031555)*t+0.00163801)*t-0.00362018)*t-0.03988024)*t+0.39894228;
194  i1 = dtmp*dtmp1;
195  }
196  if (x < 2.0) {
197  t = 0.5*x;
198  t2 = t*t; // already calculated above
199  dtmp = (((((0.0000074*t2+0.0001075)*t2+0.00262698)*t2+0.0348859)*t2+
200  0.23069756)*t2+0.4227842)*t2-0.57721566;
201  k0 = dtmp - i0*log(t);
202  dtmp = (((((-0.00004686*t2-0.00110404)*t2-0.01919402)*t2-
203  0.18156897)*t2-0.67278578)*t2+0.15443144)*t2+1.0;
204  k1 = dtmp/x + i1*log(t);
205  }
206  else {
207  t = 2.0/x;
208  dtmp1 = exp(-x)/sqrt(x);
209  dtmp = (((((0.00053208*t-0.0025154)*t+0.00587872)*t-0.01062446)*t+
210  0.02189568)*t-0.07832358)*t+1.25331414;
211  k0 = dtmp*dtmp1;
212  dtmp = (((((-0.00068245*t+0.00325614)*t-0.00780353)*t+0.01504268)*t-
213  0.0365562)*t+0.23498619)*t+1.25331414;
214  k1 = dtmp*dtmp1;
215  }
216  i0p = i1;
217  i1p = i0 - i1/x;
218  k0p = -k1;
219  k1p = -k0 - k1/x;
220  return 0;
221 }
plotPSD.t
t
Definition: plotPSD.py:36

References Eigen::bfloat16_impl::exp(), Eigen::bfloat16_impl::log(), sqrt(), plotPSD::t, and plotDoE::x.

◆ bessikna()

int CRBond_Bessel::bessikna ( int  n,
double  x,
int &  nm,
double *  in,
double *  kn,
double *  inp,
double *  knp 
)
224 {
225  double bi0,bi1,bk0,bk1,g,g0,g1,f,f0,f1,h,h0,h1,s0;
226  int k,m,ecode;
227 
228  if ((x < 0.0) || (n < 0)) return 1;
229  if (x < eps) {
230  for (k=0;k<=n;k++) {
231  in[k] = 0.0;
232  kn[k] = 1e308;
233  inp[k] = 0.0;
234  knp[k] = -1e308;
235  }
236  in[0] = 1.0;
237  inp[1] = 0.5;
238  return 0;
239  }
240  nm = n;
241  ecode = bessik01a(x,in[0],in[1],kn[0],kn[1],inp[0],inp[1],knp[0],knp[1]);
242  if (n < 2) return 0;
243  bi0 = in[0];
244  bi1 = in[1];
245  bk0 = kn[0];
246  bk1 = kn[1];
247  if ((x > 40.0) && (n < (int)(0.25*x))) {
248  h0 = bi0;
249  h1 = bi1;
250  for (k=2;k<=n;k++) {
251  h = -2.0*(k-1.0)*h1/x+h0;
252  in[k] = h;
253  h0 = h1;
254  h1 = h;
255  }
256  }
257  else {
258  m = msta1(x,200);
259  if (m < n) nm = m;
260  else m = msta2(x,n,15);
261  f0 = 0.0;
262  f1 = 1.0e-100;
263  for (k=m;k>=0;k--) {
264  f = 2.0*(k+1.0)*f1/x+f0;
265  if (x <= nm) in[k] = f;
266  f0 = f1;
267  f1 = f;
268  }
269  s0 = bi0/f;
270  for (k=0;k<=m;k++) {
271  in[k] *= s0;
272  }
273  }
274  g0 = bk0;
275  g1 = bk1;
276  for (k=2;k<=nm;k++) {
277  g = 2.0*(k-1.0)*g1/x+g0;
278  kn[k] = g;
279  g0 = g1;
280  g1 = g;
281  }
282  for (k=2;k<=nm;k++) {
283  inp[k] = in[k-1]-k*in[k]/x;
284  knp[k] = -kn[k-1]-k*kn[k]/x;
285  }
286  return 0;
287 }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
m
int * m
Definition: level2_cplx_impl.h:294
CRBond_Bessel::msta1
int msta1(double x, int mp)
Definition: crbond_bessel.cc:779
CRBond_Bessel::bessik01a
int bessik01a(double x, double &i0, double &i1, double &k0, double &k1, double &i0p, double &i1p, double &k0p, double &k1p)
Definition: crbond_bessel.cc:41
CRBond_Bessel::msta2
int msta2(double x, int n, int mp)
Definition: crbond_bessel.cc:802
Global_parameters::f1
double f1(const Vector< double > &coord)
f1 function, in front of the C1 unknown
Definition: poisson/poisson_with_singularity/two_d_poisson.cc:147

References bessik01a(), eps, f(), Global_parameters::f1(), k, m, msta1(), msta2(), n, and plotDoE::x.

◆ bessiknb()

int CRBond_Bessel::bessiknb ( int  n,
double  x,
int &  nm,
double *  in,
double *  kn,
double *  inp,
double *  knp 
)
290 {
291  double s0,bs,f,f0,f1,sk0,a0,bkl,vt,r,g,g0,g1;
292  int k,kz,m,l;
293 
294  if ((x < 0.0) || (n < 0)) return 1;
295  if (x < eps) {
296  for (k=0;k<=n;k++) {
297  in[k] = 0.0;
298  kn[k] = 1e308;
299  inp[k] = 0.0;
300  knp[k] = -1e308;
301  }
302  in[0] = 1.0;
303  inp[1] = 0.5;
304  return 0;
305  }
306  nm = n;
307  if (n == 0) nm = 1;
308  m = msta1(x,200);
309  if (m < nm) nm = m;
310  else m = msta2(x,nm,15);
311  bs = 0.0;
312  sk0 = 0.0;
313  f0 = 0.0;
314  f1 = 1.0e-100;
315  for (k=m;k>=0;k--) {
316  f = 2.0*(k+1.0)*f1/x+f0;
317  if (k <= nm) in[k] = f;
318  if ((k != 0) && (k == 2*(int)(k/2))) {
319  sk0 += 4.0*f/k;
320  }
321  bs += 2.0*f;
322  f0 = f1;
323  f1 = f;
324  }
325  s0 = exp(x)/(bs-f);
326  for (k=0;k<=nm;k++) {
327  in[k] *= s0;
328  }
329  if (x <= 8.0) {
330  kn[0] = -(log(0.5*x)+el)*in[0]+s0*sk0;
331  kn[1] = (1.0/x-in[1]*kn[0])/in[0];
332  }
333  else {
334  a0 = sqrt(M_PI_2/x)*exp(-x);
335  if (x >= 200.0) kz = 6;
336  else if (x >= 80.0) kz = 8;
337  else if (x >= 25.0) kz = 10;
338  else kz = 16;
339  for (l=0;l<2;l++) {
340  bkl = 1.0;
341  vt = 4.0*l;
342  r = 1.0;
343  for (k=1;k<=kz;k++) {
344  r *= 0.125*(vt-pow(2.0*k-1.0,2))/(k*x);
345  bkl += r;
346  }
347  kn[l] = a0*bkl;
348  }
349  }
350  g0 = kn[0];
351  g1 = kn[1];
352  for (k=2;k<=nm;k++) {
353  g = 2.0*(k-1.0)*g1/x+g0;
354  kn[k] = g;
355  g0 = g1;
356  g1 = g;
357  }
358  inp[0] = in[1];
359  knp[0] = -kn[1];
360  for (k=1;k<=nm;k++) {
361  inp[k] = in[k-1]-k*in[k]/x;
362  knp[k] = -kn[k-1]-k*kn[k]/x;
363  }
364  return 0;
365 }
int
return int(ret)+1

References el, eps, Eigen::bfloat16_impl::exp(), f(), Global_parameters::f1(), int(), k, Eigen::bfloat16_impl::log(), m, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, sqrt(), and plotDoE::x.

◆ bessikv()

int CRBond_Bessel::bessikv ( double  v,
double  x,
double &  vm,
double *  iv,
double *  kv,
double *  ivp,
double *  kvp 
)
376 {
377  double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum;
378  double f,f1,f2,ct,cs,wa,gan,ww,w0,v0n;
379  double r1,r2,bk0,bk1,bk2,a2,cb;
380  int n,k,kz,m;
381 
382  if ((v < 0.0) || (x < 0.0)) return 1;
383  x2 = x*x;
384  n = (int)v;
385  v0 = v-n;
386  if (n == 0) n = 1;
387  if (x == 0.0) {
388  for (k=0;k<=n;k++) {
389  iv[k] = 0.0;
390  kv[k] = -1e308;
391  ivp[k] = 0.0;
392  kvp[k] = 1e308;
393  }
394  if (v0 == 0.0) {
395  iv[0] = 1.0;
396  ivp[1] = 0.5;
397  }
398  vm = v;
399  return 0;
400  }
401  piv = M_PI*v0;
402  vt = 4.0*v0*v0;
403  if (v0 == 0.0) {
404  a1 = 1.0;
405  }
406  else {
407  v0p = 1.0+v0;
408  gap = gamma(v0p);
409  a1 = pow(0.5*x,v0)/gap;
410  }
411  if (x >= 50.0) kz = 8;
412  else if (x >= 35.0) kz = 10;
413  else kz = 14;
414  if (x <= 18.0) {
415  bi0 = 1.0;
416  r = 1.0;
417  for (k=1;k<=30;k++) {
418  r *= 0.25*x2/(k*(k+v0));
419  bi0 += r;
420  if (fabs(r/bi0) < eps) break;
421  }
422  bi0 *= a1;
423  }
424  else {
425  ca = exp(x)/sqrt(2.0*M_PI*x);
426  sum = 1.0;
427  r = 1.0;
428  for (k=1;k<=kz;k++) {
429  r *= -0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
430  sum += r;
431  }
432  bi0 = ca*sum;
433  }
434  m = msta1(x,200);
435  if (m < n) n = m;
436  else m = msta2(x,n,15);
437  f2 = 0.0;
438  f1 = 1.0e-100;
439  for (k=m;k>=0;k--) {
440  f = 2.0*(v0+k+1.0)*f1/x+f2;
441  if (k <= n) iv[k] = f;
442  f2 = f1;
443  f1 = f;
444  }
445  cs = bi0/f;
446  for (k=0;k<=n;k++) {
447  iv[k] *= cs;
448  }
449  ivp[0] = v0*iv[0]/x+iv[1];
450  for (k=1;k<=n;k++) {
451  ivp[k] = -(k+v0)*iv[k]/x+iv[k-1];
452  }
453  ww = 0.0;
454  if (x <= 9.0) {
455  if (v0 == 0.0) {
456  ct = -log(0.5*x)-el;
457  cs = 0.0;
458  w0 = 0.0;
459  r = 1.0;
460  for (k=1;k<=50;k++) {
461  w0 += 1.0/k;
462  r *= 0.25*x2/(k*k);
463  cs += r*(w0+ct);
464  wa = fabs(cs);
465  if (fabs((wa-ww)/wa) < eps) break;
466  ww = wa;
467  }
468  bk0 = ct+cs;
469  }
470  else {
471  v0n = 1.0-v0;
472  gan = gamma(v0n);
473  a2 = 1.0/(gan*pow(0.5*x,v0));
474  a1 = pow(0.5*x,v0)/gap;
475  sum = a2-a1;
476  r1 = 1.0;
477  r2 = 1.0;
478  for (k=1;k<=120;k++) {
479  r1 *= 0.25*x2/(k*(k-v0));
480  r2 *= 0.25*x2/(k*(k+v0));
481  sum += a2*r1-a1*r2;
482  wa = fabs(sum);
483  if (fabs((wa-ww)/wa) < eps) break;
484  ww = wa;
485  }
486  bk0 = M_PI_2*sum/sin(piv);
487  }
488  }
489  else {
490  cb = exp(-x)*sqrt(M_PI_2/x);
491  sum = 1.0;
492  r = 1.0;
493  for (k=1;k<=kz;k++) {
494  r *= 0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x);
495  sum += r;
496  }
497  bk0 = cb*sum;
498  }
499  bk1 = (1.0/x-iv[1]*bk0)/iv[0];
500  kv[0] = bk0;
501  kv[1] = bk1;
502  for (k=2;k<=n;k++) {
503  bk2 = 2.0*(v0+k-1.0)*bk1/x+bk0;
504  kv[k] = bk2;
505  bk0 = bk1;
506  bk1 = bk2;
507  }
508  kvp[0] = v0*kv[0]/x-kv[1];
509  for (k=1;k<=n;k++) {
510  kvp[k] = -(k+v0)*kv[k]/x-kv[k-1];
511  }
512  vm = n+v0;
513  return 0;
514 }
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
CRBond_Bessel::gamma
double gamma(double x)
Definition: crbond_bessel.cc:2402
Global_parameters::f2
double f2(const Vector< double > &coord)
f2 function, in front of the C2 unknown
Definition: poisson/poisson_with_singularity/two_d_poisson.cc:233

References el, eps, Eigen::bfloat16_impl::exp(), f(), Global_parameters::f1(), Global_parameters::f2(), boost::multiprecision::fabs(), gamma(), int(), k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, sin(), sqrt(), v, plotDoE::x, and Global_parameters::x2().

◆ bessjy01a()

int CRBond_Bessel::bessjy01a ( double  x,
double &  j0,
double &  j1,
double &  y0,
double &  y1,
double &  j0p,
double &  j1p,
double &  y0p,
double &  y1p 
)
553 {
554  double x2,r,ec,w0,w1,r0,r1,cs0,cs1;
555  double cu,p0,q0,p1,q1,t1,t2;
556  int k,kz;
557  static double a[] = {
558  -7.03125e-2,
559  0.112152099609375,
560  -0.5725014209747314,
561  6.074042001273483,
562  -1.100171402692467e2,
563  3.038090510922384e3,
564  -1.188384262567832e5,
565  6.252951493434797e6,
566  -4.259392165047669e8,
567  3.646840080706556e10,
568  -3.833534661393944e12,
569  4.854014686852901e14,
570  -7.286857349377656e16,
571  1.279721941975975e19};
572  static double b[] = {
573  7.32421875e-2,
574  -0.2271080017089844,
575  1.727727502584457,
576  -2.438052969955606e1,
577  5.513358961220206e2,
578  -1.825775547429318e4,
579  8.328593040162893e5,
580  -5.006958953198893e7,
581  3.836255180230433e9,
582  -3.649010818849833e11,
583  4.218971570284096e13,
584  -5.827244631566907e15,
585  9.476288099260110e17,
586  -1.792162323051699e20};
587  static double a1[] = {
588  0.1171875,
589  -0.1441955566406250,
590  0.6765925884246826,
591  -6.883914268109947,
592  1.215978918765359e2,
593  -3.302272294480852e3,
594  1.276412726461746e5,
595  -6.656367718817688e6,
596  4.502786003050393e8,
597  -3.833857520742790e10,
598  4.011838599133198e12,
599  -5.060568503314727e14,
600  7.572616461117958e16,
601  -1.326257285320556e19};
602  static double b1[] = {
603  -0.1025390625,
604  0.2775764465332031,
605  -1.993531733751297,
606  2.724882731126854e1,
607  -6.038440767050702e2,
608  1.971837591223663e4,
609  -8.902978767070678e5,
610  5.310411010968522e7,
611  -4.043620325107754e9,
612  3.827011346598605e11,
613  -4.406481417852278e13,
614  6.065091351222699e15,
615  -9.833883876590679e17,
616  1.855045211579828e20};
617 
618  if (x < 0.0) return 1;
619  if (x == 0.0) {
620  j0 = 1.0;
621  j1 = 0.0;
622  y0 = -1e308;
623  y1 = -1e308;
624  j0p = 0.0;
625  j1p = 0.5;
626  y0p = 1e308;
627  y1p = 1e308;
628  return 0;
629  }
630  x2 = x*x;
631  if (x <= 12.0) {
632  j0 = 1.0;
633  r = 1.0;
634  for (k=1;k<=30;k++) {
635  r *= -0.25*x2/(k*k);
636  j0 += r;
637  if (fabs(r) < fabs(j0)*1e-15) break;
638  }
639  j1 = 1.0;
640  r = 1.0;
641  for (k=1;k<=30;k++) {
642  r *= -0.25*x2/(k*(k+1));
643  j1 += r;
644  if (fabs(r) < fabs(j1)*1e-15) break;
645  }
646  j1 *= 0.5*x;
647  ec = log(0.5*x)+el;
648  cs0 = 0.0;
649  w0 = 0.0;
650  r0 = 1.0;
651  for (k=1;k<=30;k++) {
652  w0 += 1.0/k;
653  r0 *= -0.25*x2/(k*k);
654  r = r0 * w0;
655  cs0 += r;
656  if (fabs(r) < fabs(cs0)*1e-15) break;
657  }
658  y0 = M_2_PI*(ec*j0-cs0);
659  cs1 = 1.0;
660  w1 = 0.0;
661  r1 = 1.0;
662  for (k=1;k<=30;k++) {
663  w1 += 1.0/k;
664  r1 *= -0.25*x2/(k*(k+1));
665  r = r1*(2.0*w1+1.0/(k+1));
666  cs1 += r;
667  if (fabs(r) < fabs(cs1)*1e-15) break;
668  }
669  y1 = M_2_PI * (ec*j1-1.0/x-0.25*x*cs1);
670  }
671  else {
672  if (x >= 50.0) kz = 8; // Can be changed to 10
673  else if (x >= 35.0) kz = 10; // " " 12
674  else kz = 12; // " " 14
675  t1 = x-M_PI_4;
676  p0 = 1.0;
677  q0 = -0.125/x;
678  for (k=0;k<kz;k++) {
679  p0 += a[k]*pow(x,-2*k-2);
680  q0 += b[k]*pow(x,-2*k-3);
681  }
682  cu = sqrt(M_2_PI/x);
683  j0 = cu*(p0*cos(t1)-q0*sin(t1));
684  y0 = cu*(p0*sin(t1)+q0*cos(t1));
685  t2 = x-0.75*M_PI;
686  p1 = 1.0;
687  q1 = 0.375/x;
688  for (k=0;k<kz;k++) {
689  p1 += a1[k]*pow(x,-2*k-2);
690  q1 += b1[k]*pow(x,-2*k-3);
691  }
692  j1 = cu*(p1*cos(t2)-q1*sin(t2));
693  y1 = cu*(p1*sin(t2)+q1*cos(t2));
694  }
695  j0p = -j1;
696  j1p = j0-j1/x;
697  y0p = -y1;
698  y1p = y0-y1/x;
699  return 0;
700 }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
p0
Vector3f p0
Definition: MatrixBase_all.cpp:2
p1
Vector3f p1
Definition: MatrixBase_all.cpp:2

References a, b, cos(), e(), el, boost::multiprecision::fabs(), k, Eigen::bfloat16_impl::log(), M_PI, p0, p1, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, sin(), sqrt(), plotDoE::x, and Global_parameters::x2().

Referenced by AxisymOscillatingDisk::accel(), Global_Parameters::BesselJ(), Global_Parameters::BesselY(), bessjyna(), InterfaceProblem< ELEMENT, TIMESTEPPER >::deform_free_surface(), AxisymOscillatingDisk::position(), AxisymOscillatingDisk::residual_for_dispersion(), and AxisymOscillatingDisk::veloc().

◆ bessjy01b()

int CRBond_Bessel::bessjy01b ( double  x,
double &  j0,
double &  j1,
double &  y0,
double &  y1,
double &  j0p,
double &  j1p,
double &  y0p,
double &  y1p 
)
722 {
723  double t,t2,dtmp,a0,p0,q0,p1,q1,ta0,ta1;
724  if (x < 0.0) return 1;
725  if (x == 0.0) {
726  j0 = 1.0;
727  j1 = 0.0;
728  y0 = -1e308;
729  y1 = -1e308;
730  j0p = 0.0;
731  j1p = 0.5;
732  y0p = 1e308;
733  y1p = 1e308;
734  return 0;
735  }
736  if(x <= 4.0) {
737  t = x/4.0;
738  t2 = t*t;
739  j0 = ((((((-0.5014415e-3*t2+0.76771853e-2)*t2-0.0709253492)*t2+
740  0.4443584263)*t2-1.7777560599)*t2+3.9999973021)*t2
741  -3.9999998721)*t2+1.0;
742  j1 = t*(((((((-0.1289769e-3*t2+0.22069155e-2)*t2-0.0236616773)*t2+
743  0.1777582922)*t2-0.8888839649)*t2+2.6666660544)*t2-
744  3.999999971)*t2+1.9999999998);
745  dtmp = (((((((-0.567433e-4*t2+0.859977e-3)*t2-0.94855882e-2)*t2+
746  0.0772975809)*t2-0.4261737419)*t2+1.4216421221)*t2-
747  2.3498519931)*t2+1.0766115157)*t2+0.3674669052;
748  y0 = M_2_PI*log(0.5*x)*j0+dtmp;
749  dtmp = (((((((0.6535773e-3*t2-0.0108175626)*t2+0.107657607)*t2-
750  0.7268945577)*t2+3.1261399273)*t2-7.3980241381)*t2+
751  6.8529236342)*t2+0.3932562018)*t2-0.6366197726;
752  y1 = M_2_PI*log(0.5*x)*j1+dtmp/x;
753  }
754  else {
755  t = 4.0/x;
756  t2 = t*t;
757  a0 = sqrt(M_2_PI/x);
758  p0 = ((((-0.9285e-5*t2+0.43506e-4)*t2-0.122226e-3)*t2+
759  0.434725e-3)*t2-0.4394275e-2)*t2+0.999999997;
760  q0 = t*(((((0.8099e-5*t2-0.35614e-4)*t2+0.85844e-4)*t2-
761  0.218024e-3)*t2+0.1144106e-2)*t2-0.031249995);
762  ta0 = x-M_PI_4;
763  j0 = a0*(p0*cos(ta0)-q0*sin(ta0));
764  y0 = a0*(p0*sin(ta0)+q0*cos(ta0));
765  p1 = ((((0.10632e-4*t2-0.50363e-4)*t2+0.145575e-3)*t2
766  -0.559487e-3)*t2+0.7323931e-2)*t2+1.000000004;
767  q1 = t*(((((-0.9173e-5*t2+0.40658e-4)*t2-0.99941e-4)*t2
768  +0.266891e-3)*t2-0.1601836e-2)*t2+0.093749994);
769  ta1 = x-0.75*M_PI;
770  j1 = a0*(p1*cos(ta1)-q1*sin(ta1));
771  y1 = a0*(p1*sin(ta1)+q1*cos(ta1));
772  }
773  j0p = -j1;
774  j1p = j0-j1/x;
775  y0p = -y1;
776  y1p = y0-y1/x;
777  return 0;
778 }

References cos(), Eigen::bfloat16_impl::log(), M_PI, p0, p1, sin(), sqrt(), plotPSD::t, and plotDoE::x.

◆ bessjyna()

int CRBond_Bessel::bessjyna ( int  n,
double  x,
int &  nm,
double *  jn,
double *  yn,
double *  jnp,
double *  ynp 
)
854 {
855  double bj0,bj1,f,f0,f1,f2,cs;
856  int i,k,m,ecode;
857 
858  nm = n;
859  if ((x < 0.0) || (n < 0)) return 1;
860  if (x < 1e-15) {
861  for (i=0;i<=n;i++) {
862  jn[i] = 0.0;
863  yn[i] = -1e308;
864  jnp[i] = 0.0;
865  ynp[i] = 1e308;
866  }
867  jn[0] = 1.0;
868  jnp[1] = 0.5;
869  return 0;
870  }
871  ecode = bessjy01a(x,jn[0],jn[1],yn[0],yn[1],jnp[0],jnp[1],ynp[0],ynp[1]);
872  if (n < 2) return 0;
873  bj0 = jn[0];
874  bj1 = jn[1];
875  if (n < (int)0.9*x) {
876  for (k=2;k<=n;k++) {
877  jn[k] = 2.0*(k-1.0)*bj1/x-bj0;
878  bj0 = bj1;
879  bj1 = jn[k];
880  }
881  }
882  else {
883  m = msta1(x,200);
884  if (m < n) nm = m;
885  else m = msta2(x,n,15);
886  f2 = 0.0;
887  f1 = 1.0e-100;
888  for (k=m;k>=0;k--) {
889  f = 2.0*(k+1.0)/x*f1-f2;
890  if (k <= nm) jn[k] = f;
891  f2 = f1;
892  f1 = f;
893  }
894  if (fabs(bj0) > fabs(bj1)) cs = bj0/f;
895  else cs = bj1/f2;
896  for (k=0;k<=nm;k++) {
897  jn[k] *= cs;
898  }
899  }
900  for (k=2;k<=nm;k++) {
901  jnp[k] = jn[k-1]-k*jn[k]/x;
902  }
903  f0 = yn[0];
904  f1 = yn[1];
905  for (k=2;k<=nm;k++) {
906  f = 2.0*(k-1.0)*f1/x-f0;
907  yn[k] = f;
908  f0 = f1;
909  f1 = f;
910  }
911  for (k=2;k<=nm;k++) {
912  ynp[k] = yn[k-1]-k*yn[k]/x;
913  }
914  return 0;
915 }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
CRBond_Bessel::bessjy01a
int bessjy01a(double x, double &j0, double &j1, double &y0, double &y1, double &j0p, double &j1p, double &y0p, double &y1p)
Definition: crbond_bessel.cc:551

References bessjy01a(), e(), f(), Global_parameters::f1(), Global_parameters::f2(), boost::multiprecision::fabs(), i, k, m, msta1(), msta2(), n, and plotDoE::x.

Referenced by GlobalParameters::get_exact_u(), Global_Parameters::Hankel_first(), oomph::Hankel_functions_for_helmholtz_problem::Hankel_first(), and GlobalParameters::prescribed_incoming_flux().

◆ bessjynb()

int CRBond_Bessel::bessjynb ( int  n,
double  x,
int &  nm,
double *  jn,
double *  yn,
double *  jnp,
double *  ynp 
)
922 {
923  double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
924  double ec,bs,byk,p0,p1,q0,q1;
925  static double a[] = {
926  -0.7031250000000000e-1,
927  0.1121520996093750,
928  -0.5725014209747314,
929  6.074042001273483};
930  static double b[] = {
931  0.7324218750000000e-1,
932  -0.2271080017089844,
933  1.727727502584457,
934  -2.438052969955606e1};
935  static double a1[] = {
936  0.1171875,
937  -0.1441955566406250,
938  0.6765925884246826,
939  -6.883914268109947};
940  static double b1[] = {
941  -0.1025390625,
942  0.2775764465332031,
943  -1.993531733751297,
944  2.724882731126854e1};
945 
946  int i,k,m;
947  nm = n;
948  if ((x < 0.0) || (n < 0)) return 1;
949  if (x < 1e-15) {
950  for (i=0;i<=n;i++) {
951  jn[i] = 0.0;
952  yn[i] = -1e308;
953  jnp[i] = 0.0;
954  ynp[i] = 1e308;
955  }
956  jn[0] = 1.0;
957  jnp[1] = 0.5;
958  return 0;
959  }
960  if (x <= 300.0 || n > (int)(0.9*x)) {
961  if (n == 0) nm = 1;
962  m = msta1(x,200);
963  if (m < nm) nm = m;
964  else m = msta2(x,nm,15);
965  bs = 0.0;
966  su = 0.0;
967  sv = 0.0;
968  f2 = 0.0;
969  f1 = 1.0e-100;
970  for (k = m;k>=0;k--) {
971  f = 2.0*(k+1.0)/x*f1 - f2;
972  if (k <= nm) jn[k] = f;
973  if ((k == 2*(int)(k/2)) && (k != 0)) {
974  bs += 2.0*f;
975 // su += pow(-1,k>>1)*f/(double)k;
976  su += (-1)*((k & 2)-1)*f/(double)k;
977  }
978  else if (k > 1) {
979 // sv += pow(-1,k>>1)*k*f/(k*k-1.0);
980  sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0);
981  }
982  f2 = f1;
983  f1 = f;
984  }
985  s0 = bs+f;
986  for (k=0;k<=nm;k++) {
987  jn[k] /= s0;
988  }
989  ec = log(0.5*x) +0.5772156649015329;
990  by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
991  yn[0] = by0;
992  by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
993  yn[1] = by1;
994  }
995  else {
996  t1 = x-M_PI_4;
997  p0 = 1.0;
998  q0 = -0.125/x;
999  for (k=0;k<4;k++) {
1000  p0 += a[k]*pow(x,-2*k-2);
1001  q0 += b[k]*pow(x,-2*k-3);
1002  }
1003  cu = sqrt(M_2_PI/x);
1004  bj0 = cu*(p0*cos(t1)-q0*sin(t1));
1005  by0 = cu*(p0*sin(t1)+q0*cos(t1));
1006  jn[0] = bj0;
1007  yn[0] = by0;
1008  t2 = x-0.75*M_PI;
1009  p1 = 1.0;
1010  q1 = 0.375/x;
1011  for (k=0;k<4;k++) {
1012  p1 += a1[k]*pow(x,-2*k-2);
1013  q1 += b1[k]*pow(x,-2*k-3);
1014  }
1015  bj1 = cu*(p1*cos(t2)-q1*sin(t2));
1016  by1 = cu*(p1*sin(t2)+q1*cos(t2));
1017  jn[1] = bj1;
1018  yn[1] = by1;
1019  for (k=2;k<=nm;k++) {
1020  bjk = 2.0*(k-1.0)*bj1/x-bj0;
1021  jn[k] = bjk;
1022  bj0 = bj1;
1023  bj1 = bjk;
1024  }
1025  }
1026  jnp[0] = -jn[1];
1027  for (k=1;k<=nm;k++) {
1028  jnp[k] = jn[k-1]-k*jn[k]/x;
1029  }
1030  for (k=2;k<=nm;k++) {
1031  byk = 2.0*(k-1.0)*by1/x-by0;
1032  yn[k] = byk;
1033  by0 = by1;
1034  by1 = byk;
1035  }
1036  ynp[0] = -yn[1];
1037  for (k=1;k<=nm;k++) {
1038  ynp[k] = yn[k-1]-k*yn[k]/x;
1039  }
1040  return 0;
1041 
1042 }

References a, b, cos(), e(), f(), Global_parameters::f1(), Global_parameters::f2(), i, k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, p0, p1, Eigen::bfloat16_impl::pow(), sin(), sqrt(), and plotDoE::x.

◆ bessjyv()

int CRBond_Bessel::bessjyv ( double  v,
double  x,
double &  vm,
double *  jv,
double *  yv,
double *  djv,
double *  dyv 
)
1052 {
1053  double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2;
1054  double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq;
1055  double b,ec,w0,w1,bju0,bju1,pv0,pv1,byvk;
1056  int j,k,l,m,n,kz;
1057 
1058  x2 = x*x;
1059  n = (int)v;
1060  v0 = v-n;
1061  if ((x < 0.0) || (v < 0.0)) return 1;
1062  if (x < 1e-15) {
1063  for (k=0;k<=n;k++) {
1064  jv[k] = 0.0;
1065  yv[k] = -1e308;
1066  djv[k] = 0.0;
1067  dyv[k] = 1e308;
1068  if (v0 == 0.0) {
1069  jv[0] = 1.0;
1070  djv[1] = 0.5;
1071  }
1072  else djv[0] = 1e308;
1073  }
1074  vm = v;
1075  return 0;
1076  }
1077  if (x <= 12.0) {
1078  for (l=0;l<2;l++) {
1079  vl = v0 + l;
1080  bjvl = 1.0;
1081  r = 1.0;
1082  for (k=1;k<=40;k++) {
1083  r *= -0.25*x2/(k*(k+vl));
1084  bjvl += r;
1085  if (fabs(r) < fabs(bjvl)*1e-15) break;
1086  }
1087  vg = 1.0 + vl;
1088  a = pow(0.5*x,vl)/gamma(vg);
1089  if (l == 0) bjv0 = bjvl*a;
1090  else bjv1 = bjvl*a;
1091  }
1092  }
1093  else {
1094  if (x >= 50.0) kz = 8;
1095  else if (x >= 35.0) kz = 10;
1096  else kz = 11;
1097  for (j=0;j<2;j++) {
1098  vv = 4.0*(j+v0)*(j+v0);
1099  px = 1.0;
1100  rp = 1.0;
1101  for (k=1;k<=kz;k++) {
1102  rp *= (-0.78125e-2)*(vv-pow(4.0*k-3.0,2.0))*
1103  (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*x2);
1104  px += rp;
1105  }
1106  qx = 1.0;
1107  rq = 1.0;
1108  for (k=1;k<=kz;k++) {
1109  rq *= (-0.78125e-2)*(vv-pow(4.0*k-1.0,2.0))*
1110  (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*x2);
1111  qx += rq;
1112  }
1113  qx *= 0.125*(vv-1.0)/x;
1114  xk = x-(0.5*(j+v0)+0.25)*M_PI;
1115  a0 = sqrt(M_2_PI/x);
1116  ck = cos(xk);
1117  sk = sin(xk);
1118 
1119  if (j == 0) {
1120  bjv0 = a0*(px*ck-qx*sk);
1121  byv0 = a0*(px*sk+qx*ck);
1122  }
1123  else {
1124  bjv1 = a0*(px*ck-qx*sk);
1125  byv1 = a0*(px*sk+qx*ck);
1126  }
1127  }
1128  }
1129  jv[0] = bjv0;
1130  jv[1] = bjv1;
1131  djv[0] = v0*jv[0]/x-jv[1];
1132  djv[1] = -(1.0+v0)*jv[1]/x+jv[0];
1133  if ((n >= 2) && (n <= (int)(0.9*x))) {
1134  f0 = bjv0;
1135  f1 = bjv1;
1136  for (k=2;k<=n;k++) {
1137  f = 2.0*(k+v0-1.0)*f1/x-f0;
1138  jv[k] = f;
1139  f0 = f1;
1140  f1 = f;
1141  }
1142  }
1143  else if (n >= 2) {
1144  m = msta1(x,200);
1145  if (m < n) n = m;
1146  else m = msta2(x,n,15);
1147  f2 = 0.0;
1148  f1 = 1.0e-100;
1149  for (k=m;k>=0;k--) {
1150  f = 2.0*(v0+k+1.0)/x*f1-f2;
1151  if (k <= n) jv[k] = f;
1152  f2 = f1;
1153  f1 = f;
1154  }
1155  if (fabs(bjv0) > fabs(bjv1)) cs = bjv0/f;
1156  else cs = bjv1/f2;
1157  for (k=0;k<=n;k++) {
1158  jv[k] *= cs;
1159  }
1160  }
1161  for (k=2;k<=n;k++) {
1162  djv[k] = -(k+v0)*jv[k]/x+jv[k-1];
1163  }
1164  if (x <= 12.0) {
1165  if (v0 != 0.0) {
1166  for (l=0;l<2;l++) {
1167  vl = v0 +l;
1168  bjvl = 1.0;
1169  r = 1.0;
1170  for (k=1;k<=40;k++) {
1171  r *= -0.25*x2/(k*(k-vl));
1172  bjvl += r;
1173  if (fabs(r) < fabs(bjvl)*1e-15) break;
1174  }
1175  vg = 1.0-vl;
1176  b = pow(2.0/x,vl)/gamma(vg);
1177  if (l == 0) bju0 = bjvl*b;
1178  else bju1 = bjvl*b;
1179  }
1180  pv0 = M_PI*v0;
1181  pv1 = M_PI*(1.0+v0);
1182  byv0 = (bjv0*cos(pv0)-bju0)/sin(pv0);
1183  byv1 = (bjv1*cos(pv1)-bju1)/sin(pv1);
1184  }
1185  else {
1186  ec = log(0.5*x)+el;
1187  cs0 = 0.0;
1188  w0 = 0.0;
1189  r0 = 1.0;
1190  for (k=1;k<=30;k++) {
1191  w0 += 1.0/k;
1192  r0 *= -0.25*x2/(k*k);
1193  cs0 += r0*w0;
1194  }
1195  byv0 = M_2_PI*(ec*bjv0-cs0);
1196  cs1 = 1.0;
1197  w1 = 0.0;
1198  r1 = 1.0;
1199  for (k=1;k<=30;k++) {
1200  w1 += 1.0/k;
1201  r1 *= -0.25*x2/(k*(k+1));
1202  cs1 += r1*(2.0*w1+1.0/(k+1.0));
1203  }
1204  byv1 = M_2_PI*(ec*bjv1-1.0/x-0.25*x*cs1);
1205  }
1206  }
1207  yv[0] = byv0;
1208  yv[1] = byv1;
1209  for (k=2;k<=n;k++) {
1210  byvk = 2.0*(v0+k-1.0)*byv1/x-byv0;
1211  yv[k] = byvk;
1212  byv0 = byv1;
1213  byv1 = byvk;
1214  }
1215  dyv[0] = v0*yv[0]/x-yv[1];
1216  for (k=1;k<=n;k++) {
1217  dyv[k] = -(k+v0)*yv[k]/x+yv[k-1];
1218  }
1219  vm = n + v0;
1220  return 0;
1221 }
px
RealScalar RealScalar * px
Definition: level1_cplx_impl.h:27
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References a, b, cos(), e(), el, f(), Global_parameters::f1(), Global_parameters::f2(), boost::multiprecision::fabs(), gamma(), int(), j, k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), px, UniformPSDSelfTest::r, sin(), sqrt(), v, plotDoE::x, and Global_parameters::x2().

Referenced by ProblemParameters::exact_minus_dudr(), PlanarWave::get_exact_u(), ProblemParameters::get_exact_u(), main(), and oomph::FourierDecomposedHelmholtzDtNMesh< ELEMENT >::setup_gamma().

◆ cbessik01()

int CRBond_Bessel::cbessik01 ( complex< double >  z,
complex< double > &  ci0,
complex< double > &  ci1,
complex< double > &  ck0,
complex< double > &  ck1,
complex< double > &  ci0p,
complex< double > &  ci1p,
complex< double > &  ck0p,
complex< double > &  ck1p 
)
1240 {
1241  complex<double> z1,z2,zr,zr2,cr,ca,cb,cs,ct,cw;
1242  double a0,w0;
1243  int k,kz;
1244  static double a[] = {
1245  0.125,
1246  7.03125e-2,
1247  7.32421875e-2,
1248  1.1215209960938e-1,
1249  2.2710800170898e-1,
1250  5.7250142097473e-1,
1251  1.7277275025845,
1252  6.0740420012735,
1253  2.4380529699556e1,
1254  1.1001714026925e2,
1255  5.5133589612202e2,
1256  3.0380905109224e3};
1257  static double b[] = {
1258  -0.375,
1259  -1.171875e-1,
1260  -1.025390625e-1,
1261  -1.4419555664063e-1,
1262  -2.7757644653320e-1,
1263  -6.7659258842468e-1,
1264  -1.9935317337513,
1265  -6.8839142681099,
1266  -2.7248827311269e1,
1267  -1.2159789187654e2,
1268  -6.0384407670507e2,
1269  -3.3022722944809e3};
1270  static double a1[] = {
1271  0.125,
1272  0.2109375,
1273  1.0986328125,
1274  1.1775970458984e1,
1275  2.1461706161499e2,
1276  5.9511522710323e3,
1277  2.3347645606175e5,
1278  1.2312234987631e7,
1279  8.401390346421e08,
1280  7.2031420482627e10};
1281 
1282  a0 = abs(z);
1283  z2 = z*z;
1284  z1 = z;
1285  if (a0 == 0.0) {
1286  ci0 = cone;
1287  ci1 = czero;
1288  ck0 = complex<double> (1e308,0);
1289  ck1 = complex<double> (1e308,0);
1290  ci0p = czero;
1291  ci1p = complex<double>(0.5,0.0);
1292  ck0p = complex<double>(-1e308,0);
1293  ck1p = complex<double>(-1e308,0);
1294  return 0;
1295  }
1296  if (real(z) < 0.0) z1 = -z;
1297  if (a0 <= 18.0) {
1298  ci0 = cone;
1299  cr = cone;
1300  for (k=1;k<=50;k++) {
1301  cr *= 0.25*z2/(double)(k*k);
1302  ci0 += cr;
1303  if (abs(cr/ci0) < eps) break;
1304  }
1305  ci1 = cone;
1306  cr = cone;
1307  for (k=1;k<=50;k++) {
1308  cr *= 0.25*z2/(double)(k*(k+1.0));
1309  ci1 += cr;
1310  if (abs(cr/ci1) < eps) break;
1311  }
1312  ci1 *= 0.5*z1;
1313  }
1314  else {
1315  if (a0 >= 50.0) kz = 7;
1316  else if (a0 >= 35.0) kz = 9;
1317  else kz = 12;
1318  ca = exp(z1)/sqrt(2.0*M_PI*z1);
1319  ci0 = cone;
1320  zr = 1.0/z1;
1321  for (k=0;k<kz;k++) {
1322  ci0 += a[k]*pow(zr,k+1.0);
1323  }
1324  ci0 *= ca;
1325  ci1 = cone;
1326  for (k=0;k<kz;k++) {
1327  ci1 += b[k]*pow(zr,k+1.0);
1328  }
1329  ci1 *= ca;
1330  }
1331  if (a0 <= 9.0) {
1332  cs = czero;
1333  ct = -log(0.5*z1)-el;
1334  w0 = 0.0;
1335  cr = cone;
1336  for (k=1;k<=50;k++) {
1337  w0 += 1.0/k;
1338  cr *= 0.25*z2/(double)(k*k);
1339  cs += cr*(w0+ct);
1340  if (abs((cs-cw)/cs) < eps) break;
1341  cw = cs;
1342  }
1343  ck0 = ct+cs;
1344  }
1345  else {
1346  cb = 0.5/z1;
1347  zr2 = 1.0/z2;
1348  ck0 = cone;
1349  for (k=0;k<10;k++) {
1350  ck0 += a1[k]*pow(zr2,k+1.0);
1351  }
1352  ck0 *= cb/ci0;
1353  }
1354  ck1 = (1.0/z1 - ci1*ck0)/ci0;
1355  if (real(z) < 0.0) {
1356  if (imag(z) < 0.0) {
1357  ck0 += cii*M_PI*ci0;
1358  ck1 = -ck1+cii*M_PI*ci1;
1359  }
1360  else if (imag(z) > 0.0) {
1361  ck0 -= cii*M_PI*ci0;
1362  ck1 = -ck1-cii*M_PI*ci1;
1363  }
1364  ci1 = -ci1;
1365  }
1366  ci0p = ci1;
1367  ci1p = ci0-1.0*ci1/z;
1368  ck0p = -ck1;
1369  ck1p = -ck0-1.0*ck1/z;
1370  return 0;
1371 }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
imag
AnnoyingScalar imag(const AnnoyingScalar &)
Definition: AnnoyingScalar.h:132
real
float real
Definition: datatypes.h:10
CRBond_Bessel::cone
static complex< double > cone(1.0, 0.0)
CRBond_Bessel::czero
static complex< double > czero(0.0, 0.0)
CRBond_Bessel::cii
static complex< double > cii(0.0, 1.0)

References a, abs(), b, cii(), cone(), czero(), el, eps, Eigen::bfloat16_impl::exp(), imag(), k, Eigen::bfloat16_impl::log(), M_PI, Eigen::bfloat16_impl::pow(), and sqrt().

Referenced by cbessikna().

◆ cbessikna()

int CRBond_Bessel::cbessikna ( int  n,
complex< double >  z,
int &  nm,
complex< double > *  ci,
complex< double > *  ck,
complex< double > *  cip,
complex< double > *  ckp 
)
1374 {
1375  complex<double> ci0,ci1,ck0,ck1,ckk,cf,cf1,cf2,cs;
1376  double a0;
1377  int k,m,ecode;
1378  a0 = abs(z);
1379  nm = n;
1380  if (a0 < 1.0e-100) {
1381  for (k=0;k<=n;k++) {
1382  ci[k] = czero;
1383  ck[k] = complex<double>(-1e308,0);
1384  cip[k] = czero;
1385  ckp[k] = complex<double>(1e308,0);
1386  }
1387  ci[0] = complex<double>(1e308,0);
1388  cip[1] = complex<double>(0.5,0.0);
1389  return 0;
1390  }
1391  ecode = cbessik01(z,ci[0],ci[1],ck[0],ck[1],cip[0],cip[1],ckp[0],ckp[1]);
1392  if (n < 2) return 0;
1393  ci0 = ci[0];
1394  ci1 = ci[1];
1395  ck0 = ck[0];
1396  ck1 = ck[1];
1397  m = msta1(a0,200);
1398  if (m < n) nm = m;
1399  else m = msta2(a0,n,15);
1400  cf2 = czero;
1401  cf1 = complex<double>(1.0e-100,0.0);
1402  for (k=m;k>=0;k--) {
1403  cf = 2.0*(k+1.0)*cf1/z+cf2;
1404  if (k <= nm) ci[k] = cf;
1405  cf2 = cf1;
1406  cf1 = cf;
1407  }
1408  cs = ci0/cf;
1409  for (k=0;k<=nm;k++) {
1410  ci[k] *= cs;
1411  }
1412  for (k=2;k<=nm;k++) {
1413  if (abs(ci[k-1]) > abs(ci[k-2])) {
1414  ckk = (1.0/z-ci[k]*ck[k-1])/ci[k-1];
1415  }
1416  else {
1417  ckk = (ci[k]*ck[k-2]+2.0*(k-1.0)/(z*z))/ci[k-2];
1418  }
1419  ck[k] = ckk;
1420  }
1421  for (k=2;k<=nm;k++) {
1422  cip[k] = ci[k-1]-(double)k*ci[k]/z;
1423  ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
1424  }
1425  return 0;
1426 }
CRBond_Bessel::cbessik01
int cbessik01(complex< double >z, complex< double > &ci0, complex< double > &ci1, complex< double > &ck0, complex< double > &ck1, complex< double > &ci0p, complex< double > &ci1p, complex< double > &ck0p, complex< double > &ck1p)
Definition: crbond_bessel.cc:1237

References abs(), cbessik01(), czero(), k, m, msta1(), msta2(), and n.

◆ cbessiknb()

int CRBond_Bessel::cbessiknb ( int  n,
complex< double >  z,
int &  nm,
complex< double > *  ci,
complex< double > *  ck,
complex< double > *  cip,
complex< double > *  ckp 
)
1429 {
1430  complex<double> z1,cbs,csk0,cf,cf0,cf1,ca0,cbkl;
1431  complex<double> cg,cg0,cg1,cs0,cs,cr;
1432  double a0,vt,fac;
1433  int k,kz,l,m;
1434 
1435  a0 = abs(z);
1436  nm = n;
1437  if (a0 < 1.0e-100) {
1438  for (k=0;k<=n;k++) {
1439  ci[k] = czero;
1440  ck[k] = complex<double>(1e308,0);
1441  cip[k] = czero;
1442  ckp[k] = complex<double>(-1e308,0);
1443  }
1444  ci[0] = complex<double>(1.0,0.0);
1445  cip[1] = complex<double>(0.5,0.0);
1446  return 0;
1447  }
1448  z1 = z;
1449  if (real(z) < 0.0) z1 = -z;
1450  if (n == 0) nm = 1;
1451  m = msta1(a0,200);
1452  if (m < nm) nm = m;
1453  else m = msta2(a0,nm,15);
1454  cbs = czero;
1455  csk0 = czero;
1456  cf0 = czero;
1457  cf1 = complex<double>(1.0e-100,0.0);
1458  for (k=m;k>=0;k--) {
1459  cf = 2.0*(k+1.0)*cf1/z1+cf0;
1460  if (k <=nm) ci[k] = cf;
1461  if ((k != 0) && (k == 2*(k>>1))) csk0 += 4.0*cf/(double)k;
1462  cbs += 2.0*cf;
1463  cf0 = cf1;
1464  cf1 = cf;
1465  }
1466  cs0 = exp(z1)/(cbs-cf);
1467  for (k=0;k<=nm;k++) {
1468  ci[k] *= cs0;
1469  }
1470  if (a0 <= 9.0) {
1471  ck[0] = -(log(0.5*z1)+el)*ci[0]+cs0*csk0;
1472  ck[1] = (1.0/z1-ci[1]*ck[0])/ci[0];
1473  }
1474  else {
1475  ca0 = sqrt(M_PI_2/z1)*exp(-z1);
1476  if (a0 >= 200.0) kz = 6;
1477  else if (a0 >= 80.0) kz = 8;
1478  else if (a0 >= 25.0) kz = 10;
1479  else kz = 16;
1480  for (l=0;l<2;l++) {
1481  cbkl = cone;
1482  vt = 4.0*l;
1483  cr = cone;
1484  for (k=1;k<=kz;k++) {
1485  cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
1486  cbkl += cr;
1487  }
1488  ck[l] = ca0*cbkl;
1489  }
1490  }
1491  cg0 = ck[0];
1492  cg1 = ck[1];
1493  for (k=2;k<=nm;k++) {
1494  cg = 2.0*(k-1.0)*cg1/z1+cg0;
1495  ck[k] = cg;
1496  cg0 = cg1;
1497  cg1 = cg;
1498  }
1499  if (real(z) < 0.0) {
1500  fac = 1.0;
1501  for (k=0;k<=nm;k++) {
1502  if (imag(z) < 0.0) {
1503  ck[k] = fac*ck[k]+cii*M_PI*ci[k];
1504  }
1505  else {
1506  ck[k] = fac*ck[k]-cii*M_PI*ci[k];
1507  }
1508  ci[k] *= fac;
1509  fac = -fac;
1510  }
1511  }
1512  cip[0] = ci[1];
1513  ckp[0] = -ck[1];
1514  for (k=1;k<=nm;k++) {
1515  cip[k] = ci[k-1]-(double)k*ci[k]/z;
1516  ckp[k] = -ck[k-1]-(double)k*ck[k]/z;
1517  }
1518  return 0;
1519 }

References abs(), cii(), cone(), czero(), el, Eigen::bfloat16_impl::exp(), imag(), k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), and sqrt().

◆ cbessikv()

int CRBond_Bessel::cbessikv ( double  v,
complex< double >  z,
double &  vm,
complex< double > *  civ,
complex< double > *  ckv,
complex< double > *  civp,
complex< double > *  ckvp 
)
1522 {
1523  complex<double> z1,z2,ca1,ca,cs,cr,ci0,cbi0,cf,cf1,cf2;
1524  complex<double> ct,cp,cbk0,ca2,cr1,cr2,csu,cws,cb;
1525  complex<double> cg0,cg1,cgk,cbk1,cvk;
1526  double a0,v0,v0p,v0n,vt,w0,piv,gap,gan;
1527  int m,n,k,kz;
1528 
1529  a0 = abs(z);
1530  z1 = z;
1531  z2 = z*z;
1532  n = (int)v;
1533  v0 = v-n;
1534  piv = M_PI*v0;
1535  vt = 4.0*v0*v0;
1536  if (n == 0) n = 1;
1537  if (a0 < 1e-100) {
1538  for (k=0;k<=n;k++) {
1539  civ[k] = czero;
1540  ckv[k] = complex<double>(-1e308,0);
1541  civp[k] = czero;
1542  ckvp[k] = complex<double>(1e308,0);
1543  }
1544  if (v0 == 0.0) {
1545  civ[0] = cone;
1546  civp[1] = complex<double> (0.5,0.0);
1547  }
1548  vm = v;
1549  return 0;
1550  }
1551  if (a0 >= 50.0) kz = 8;
1552  else if (a0 >= 35.0) kz = 10;
1553  else kz = 14;
1554  if (real(z) <= 0.0) z1 = -z;
1555  if (a0 < 18.0) {
1556  if (v0 == 0.0) {
1557  ca1 = cone;
1558  }
1559  else {
1560  v0p = 1.0+v0;
1561  gap = gamma(v0p);
1562  ca1 = pow(0.5*z1,v0)/gap;
1563  }
1564  ci0 = cone;
1565  cr = cone;
1566  for (k=1;k<=50;k++) {
1567  cr *= 0.25*z2/(k*(k+v0));
1568  ci0 += cr;
1569  if (abs(cr/ci0) < eps) break;
1570  }
1571  cbi0 = ci0*ca1;
1572  }
1573  else {
1574  ca = exp(z1)/sqrt(2.0*M_PI*z1);
1575  cs = cone;
1576  cr = cone;
1577  for (k=1;k<=kz;k++) {
1578  cr *= -0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
1579  cs += cr;
1580  }
1581  cbi0 = ca*cs;
1582  }
1583  m = msta1(a0,200);
1584  if (m < n) n = m;
1585  else m = msta2(a0,n,15);
1586  cf2 = czero;
1587  cf1 = complex<double>(1.0e-100,0.0);
1588  for (k=m;k>=0;k--) {
1589  cf = 2.0*(v0+k+1.0)*cf1/z1+cf2;
1590  if (k <= n) civ[k] = cf;
1591  cf2 = cf1;
1592  cf1 = cf;
1593  }
1594  cs = cbi0/cf;
1595  for (k=0;k<=n;k++) {
1596  civ[k] *= cs;
1597  }
1598  if (a0 <= 9.0) {
1599  if (v0 == 0.0) {
1600  ct = -log(0.5*z1)-el;
1601  cs = czero;
1602  w0 = 0.0;
1603  cr = cone;
1604  for (k=1;k<=50;k++) {
1605  w0 += 1.0/k;
1606  cr *= 0.25*z2/(double)(k*k);
1607  cp = cr*(w0+ct);
1608  cs += cp;
1609  if ((k >= 10) && (abs(cp/cs) < eps)) break;
1610  }
1611  cbk0 = ct+cs;
1612  }
1613  else {
1614  v0n = 1.0-v0;
1615  gan = gamma(v0n);
1616  ca2 = 1.0/(gan*pow(0.5*z1,v0));
1617  ca1 = pow(0.5*z1,v0)/gap;
1618  csu = ca2-ca1;
1619  cr1 = cone;
1620  cr2 = cone;
1621  cws = czero;
1622  for (k=1;k<=50;k++) {
1623  cr1 *= 0.25*z2/(k*(k-v0));
1624  cr2 *= 0.25*z2/(k*(k+v0));
1625  csu += ca2*cr1-ca1*cr2;
1626  if ((k >= 10) && (abs((cws-csu)/csu) < eps)) break;
1627  cws = csu;
1628  }
1629  cbk0 = csu*M_PI_2/sin(piv);
1630  }
1631  }
1632  else {
1633  cb = exp(-z1)*sqrt(M_PI_2/z1);
1634  cs = cone;
1635  cr = cone;
1636  for (k=1;k<=kz;k++) {
1637  cr *= 0.125*(vt-pow(2.0*k-1.0,2.0))/((double)k*z1);
1638  cs += cr;
1639  }
1640  cbk0 = cb*cs;
1641  }
1642  cbk1 = (1.0/z1-civ[1]*cbk0)/civ[0];
1643  ckv[0] = cbk0;
1644  ckv[1] = cbk1;
1645  cg0 = cbk0;
1646  cg1 = cbk1;
1647  for (k=2;k<=n;k++) {
1648  cgk = 2.0*(v0+k-1.0)*cg1/z1+cg0;
1649  ckv[k] = cgk;
1650  cg0 = cg1;
1651  cg1 = cgk;
1652  }
1653  if (real(z) < 0.0) {
1654  for (k=0;k<=n;k++) {
1655  cvk = exp((k+v0)*M_PI*cii);
1656  if (imag(z) < 0.0) {
1657  ckv[k] = cvk*ckv[k]+M_PI*cii*civ[k];
1658  civ[k] /= cvk;
1659  }
1660  else if (imag(z) > 0.0) {
1661  ckv[k] = ckv[k]/cvk-M_PI*cii*civ[k];
1662  civ[k] *= cvk;
1663  }
1664  }
1665  }
1666  civp[0] = v0*civ[0]/z+civ[1];
1667  ckvp[0] = v0*ckv[0]/z-ckv[1];
1668  for (k=1;k<=n;k++) {
1669  civp[k] = -(k+v0)*civ[k]/z+civ[k-1];
1670  ckvp[k] = -(k+v0)*ckv[k]/z-ckv[k-1];
1671  }
1672  vm = n+v0;
1673  return 0;
1674 }

References abs(), cii(), cone(), czero(), e(), el, eps, Eigen::bfloat16_impl::exp(), gamma(), imag(), int(), k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), sin(), sqrt(), and v.

◆ cbessjy01()

int CRBond_Bessel::cbessjy01 ( complex< double >  z,
complex< double > &  cj0,
complex< double > &  cj1,
complex< double > &  cy0,
complex< double > &  cy1,
complex< double > &  cj0p,
complex< double > &  cj1p,
complex< double > &  cy0p,
complex< double > &  cy1p 
)
1692 {
1693  complex<double> z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
1694  double a0,w0,w1;
1695  int k,kz;
1696 
1697  static double a[] = {
1698  -7.03125e-2,
1699  0.112152099609375,
1700  -0.5725014209747314,
1701  6.074042001273483,
1702  -1.100171402692467e2,
1703  3.038090510922384e3,
1704  -1.188384262567832e5,
1705  6.252951493434797e6,
1706  -4.259392165047669e8,
1707  3.646840080706556e10,
1708  -3.833534661393944e12,
1709  4.854014686852901e14,
1710  -7.286857349377656e16,
1711  1.279721941975975e19};
1712  static double b[] = {
1713  7.32421875e-2,
1714  -0.2271080017089844,
1715  1.727727502584457,
1716  -2.438052969955606e1,
1717  5.513358961220206e2,
1718  -1.825775547429318e4,
1719  8.328593040162893e5,
1720  -5.006958953198893e7,
1721  3.836255180230433e9,
1722  -3.649010818849833e11,
1723  4.218971570284096e13,
1724  -5.827244631566907e15,
1725  9.476288099260110e17,
1726  -1.792162323051699e20};
1727  static double a1[] = {
1728  0.1171875,
1729  -0.1441955566406250,
1730  0.6765925884246826,
1731  -6.883914268109947,
1732  1.215978918765359e2,
1733  -3.302272294480852e3,
1734  1.276412726461746e5,
1735  -6.656367718817688e6,
1736  4.502786003050393e8,
1737  -3.833857520742790e10,
1738  4.011838599133198e12,
1739  -5.060568503314727e14,
1740  7.572616461117958e16,
1741  -1.326257285320556e19};
1742  static double b1[] = {
1743  -0.1025390625,
1744  0.2775764465332031,
1745  -1.993531733751297,
1746  2.724882731126854e1,
1747  -6.038440767050702e2,
1748  1.971837591223663e4,
1749  -8.902978767070678e5,
1750  5.310411010968522e7,
1751  -4.043620325107754e9,
1752  3.827011346598605e11,
1753  -4.406481417852278e13,
1754  6.065091351222699e15,
1755  -9.833883876590679e17,
1756  1.855045211579828e20};
1757 
1758  a0 = abs(z);
1759  z2 = z*z;
1760  z1 = z;
1761  if (a0 == 0.0) {
1762  cj0 = cone;
1763  cj1 = czero;
1764  cy0 = complex<double>(-1e308,0);
1765  cy1 = complex<double>(-1e308,0);
1766  cj0p = czero;
1767  cj1p = complex<double>(0.5,0.0);
1768  cy0p = complex<double>(1e308,0);
1769  cy1p = complex<double>(1e308,0);
1770  return 0;
1771  }
1772  if (real(z) < 0.0) z1 = -z;
1773  if (a0 <= 12.0) {
1774  cj0 = cone;
1775  cr = cone;
1776  for (k=1;k<=40;k++) {
1777  cr *= -0.25*z2/(double)(k*k);
1778  cj0 += cr;
1779  if (abs(cr) < abs(cj0)*eps) break;
1780  }
1781  cj1 = cone;
1782  cr = cone;
1783  for (k=1;k<=40;k++) {
1784  cr *= -0.25*z2/(k*(k+1.0));
1785  cj1 += cr;
1786  if (abs(cr) < abs(cj1)*eps) break;
1787  }
1788  cj1 *= 0.5*z1;
1789  w0 = 0.0;
1790  cr = cone;
1791  cs = czero;
1792  for (k=1;k<=40;k++) {
1793  w0 += 1.0/k;
1794  cr *= -0.25*z2/(double)(k*k);
1795  cp = cr*w0;
1796  cs += cp;
1797  if (abs(cp) < abs(cs)*eps) break;
1798  }
1799  cy0 = M_2_PI*((log(0.5*z1)+el)*cj0-cs);
1800  w1 = 0.0;
1801  cr = cone;
1802  cs = cone;
1803  for (k=1;k<=40;k++) {
1804  w1 += 1.0/k;
1805  cr *= -0.25*z2/(k*(k+1.0));
1806  cp = cr*(2.0*w1+1.0/(k+1.0));
1807  cs += cp;
1808  if (abs(cp) < abs(cs)*eps) break;
1809  }
1810  cy1 = M_2_PI*((log(0.5*z1)+el)*cj1-1.0/z1-0.25*z1*cs);
1811  }
1812  else {
1813  if (a0 >= 50.0) kz = 8; // can be changed to 10
1814  else if (a0 >= 35.0) kz = 10; // " " " 12
1815  else kz = 12; // " " " 14
1816  ct1 = z1 - M_PI_4;
1817  cp0 = cone;
1818  for (k=0;k<kz;k++) {
1819  cp0 += a[k]*pow(z1,-2.0*k-2.0);
1820  }
1821  cq0 = -0.125/z1;
1822  for (k=0;k<kz;k++) {
1823  cq0 += b[k]*pow(z1,-2.0*k-3.0);
1824  }
1825  cu = sqrt(M_2_PI/z1);
1826  cj0 = cu*(cp0*cos(ct1)-cq0*sin(ct1));
1827  cy0 = cu*(cp0*sin(ct1)+cq0*cos(ct1));
1828  ct2 = z1 - 0.75*M_PI;
1829  cp1 = cone;
1830  for (k=0;k<kz;k++) {
1831  cp1 += a1[k]*pow(z1,-2.0*k-2.0);
1832  }
1833  cq1 = 0.375/z1;
1834  for (k=0;k<kz;k++) {
1835  cq1 += b1[k]*pow(z1,-2.0*k-3.0);
1836  }
1837  cj1 = cu*(cp1*cos(ct2)-cq1*sin(ct2));
1838  cy1 = cu*(cp1*sin(ct2)+cq1*cos(ct2));
1839  }
1840  if (real(z) < 0.0) {
1841  if (imag(z) < 0.0) {
1842  cy0 -= 2.0*cii*cj0;
1843  cy1 = -(cy1-2.0*cii*cj1);
1844  }
1845  else if (imag(z) > 0.0) {
1846  cy0 += 2.0*cii*cj0;
1847  cy1 = -(cy1+2.0*cii*cj1);
1848  }
1849  cj1 = -cj1;
1850  }
1851  cj0p = -cj1;
1852  cj1p = cj0-cj1/z;
1853  cy0p = -cy1;
1854  cy1p = cy0-cy1/z;
1855  return 0;
1856 }

References a, abs(), b, cii(), cone(), cos(), czero(), el, eps, imag(), k, Eigen::bfloat16_impl::log(), M_PI, Eigen::bfloat16_impl::pow(), sin(), and sqrt().

Referenced by cbessjyna().

◆ cbessjyna()

int CRBond_Bessel::cbessjyna ( int  n,
complex< double >  z,
int &  nm,
complex< double > *  cj,
complex< double > *  cy,
complex< double > *  cjp,
complex< double > *  cyp 
)
1860 {
1861  complex<double> cbj0,cbj1,cby0,cby1,cj0,cjk,cj1,cf,cf1,cf2;
1862  complex<double> cs,cg0,cg1,cyk,cyl1,cyl2,cylk,cp11,cp12,cp21,cp22;
1863  complex<double> ch0,ch1,ch2;
1864  double a0,yak,ya1,ya0,wa;
1865  int m,k,lb,lb0;
1866 
1867  if (n < 0) return 1;
1868  a0 = abs(z);
1869  nm = n;
1870  if (a0 < 1.0e-100) {
1871  for (k=0;k<=n;k++) {
1872  cj[k] = czero;
1873  cy[k] = complex<double> (-1e308,0);
1874  cjp[k] = czero;
1875  cyp[k] = complex<double>(1e308,0);
1876  }
1877  cj[0] = cone;
1878  cjp[1] = complex<double>(0.5,0.0);
1879  return 0;
1880  }
1881  cbessjy01(z,cj[0],cj[1],cy[0],cy[1],cjp[0],cjp[1],cyp[0],cyp[1]);
1882  cbj0 = cj[0];
1883  cbj1 = cj[1];
1884  cby0 = cy[0];
1885  cby1 = cy[1];
1886  if (n <= 1) return 0;
1887  if (n < (int)0.25*a0) {
1888  cj0 = cbj0;
1889  cj1 = cbj1;
1890  for (k=2;k<=n;k++) {
1891  cjk = 2.0*(k-1.0)*cj1/z-cj0;
1892  cj[k] = cjk;
1893  cj0 = cj1;
1894  cj1 = cjk;
1895  }
1896  }
1897  else {
1898  m = msta1(a0,200);
1899  if (m < n) nm = m;
1900  else m = msta2(a0,n,15);
1901  cf2 = czero;
1902  cf1 = complex<double> (1.0e-100,0.0);
1903  for (k=m;k>=0;k--) {
1904  cf = 2.0*(k+1.0)*cf1/z-cf2;
1905  if (k <=nm) cj[k] = cf;
1906  cf2 = cf1;
1907  cf1 = cf;
1908  }
1909  if (abs(cbj0) > abs(cbj1)) cs = cbj0/cf;
1910  else cs = cbj1/cf2;
1911  for (k=0;k<=nm;k++) {
1912  cj[k] *= cs;
1913  }
1914  }
1915  for (k=2;k<=nm;k++) {
1916  cjp[k] = cj[k-1]-(double)k*cj[k]/z;
1917  }
1918  ya0 = abs(cby0);
1919  lb = 0;
1920  cg0 = cby0;
1921  cg1 = cby1;
1922  for (k=2;k<=nm;k++) {
1923  cyk = 2.0*(k-1.0)*cg1/z-cg0;
1924  yak = abs(cyk);
1925  ya1 = abs(cg0);
1926  if ((yak < ya0) && (yak < ya1)) lb = k;
1927  cy[k] = cyk;
1928  cg0 = cg1;
1929  cg1 = cyk;
1930  }
1931  lb0 = 0;
1932  if ((lb > 4) && (imag(z) != 0.0)) {
1933  while (lb != lb0) {
1934  ch2 = cone;
1935  ch1 = czero;
1936  lb0 = lb;
1937  for (k=lb;k>=1;k--) {
1938  ch0 = 2.0*k*ch1/z-ch2;
1939  ch2 = ch1;
1940  ch1 = ch0;
1941  }
1942  cp12 = ch0;
1943  cp22 = ch2;
1944  ch2 = czero;
1945  ch1 = cone;
1946  for (k=lb;k>=1;k--) {
1947  ch0 = 2.0*k*ch1/z-ch2;
1948  ch2 = ch1;
1949  ch1 = ch0;
1950  }
1951  cp11 = ch0;
1952  cp21 = ch2;
1953  if (lb == nm)
1954  cj[lb+1] = 2.0*lb*cj[lb]/z-cj[lb-1];
1955  if (abs(cj[0]) > abs(cj[1])) {
1956  cy[lb+1] = (cj[lb+1]*cby0-2.0*cp11/(M_PI*z))/cj[0];
1957  cy[lb] = (cj[lb]*cby0+2.0*cp12/(M_PI*z))/cj[0];
1958  }
1959  else {
1960  cy[lb+1] = (cj[lb+1]*cby1-2.0*cp21/(M_PI*z))/cj[1];
1961  cy[lb] = (cj[lb]*cby1+2.0*cp22/(M_PI*z))/cj[1];
1962  }
1963  cyl2 = cy[lb+1];
1964  cyl1 = cy[lb];
1965  for (k=lb-1;k>=0;k--) {
1966  cylk = 2.0*(k+1.0)*cyl1/z-cyl2;
1967  cy[k] = cylk;
1968  cyl2 = cyl1;
1969  cyl1 = cylk;
1970  }
1971  cyl1 = cy[lb];
1972  cyl2 = cy[lb+1];
1973  for (k=lb+1;k<n;k++) {
1974  cylk = 2.0*k*cyl2/z-cyl1;
1975  cy[k+1] = cylk;
1976  cyl1 = cyl2;
1977  cyl2 = cylk;
1978  }
1979  for (k=2;k<=nm;k++) {
1980  wa = abs(cy[k]);
1981  if (wa < abs(cy[k-1])) lb = k;
1982  }
1983  }
1984  }
1985  for (k=2;k<=nm;k++) {
1986  cyp[k] = cy[k-1]-(double)k*cy[k]/z;
1987  }
1988  return 0;
1989 }
CRBond_Bessel::cbessjy01
int cbessjy01(complex< double > z, complex< double > &cj0, complex< double > &cj1, complex< double > &cy0, complex< double > &cy1, complex< double > &cj0p, complex< double > &cj1p, complex< double > &cy0p, complex< double > &cy1p)
Definition: crbond_bessel.cc:1689

References abs(), cbessjy01(), cone(), czero(), imag(), k, m, M_PI, msta1(), msta2(), and n.

Referenced by oomph::Hankel_functions_for_helmholtz_problem::CHankel_first().

◆ cbessjynb()

int CRBond_Bessel::cbessjynb ( int  n,
complex< double >  z,
int &  nm,
complex< double > *  cj,
complex< double > *  cy,
complex< double > *  cjp,
complex< double > *  cyp 
)
1993 {
1994  complex<double> cf,cf0,cf1,cf2,cbs,csu,csv,cs0,ce;
1995  complex<double> ct1,cp0,cq0,cp1,cq1,cu,cbj0,cby0,cbj1,cby1;
1996  complex<double> cyy,cbjk,ct2;
1997  double a0,y0;
1998  int k,m;
1999  static double a[] = {
2000  -0.7031250000000000e-1,
2001  0.1121520996093750,
2002  -0.5725014209747314,
2003  6.074042001273483};
2004  static double b[] = {
2005  0.7324218750000000e-1,
2006  -0.2271080017089844,
2007  1.727727502584457,
2008  -2.438052969955606e1};
2009  static double a1[] = {
2010  0.1171875,
2011  -0.1441955566406250,
2012  0.6765925884246826,
2013  -6.883914268109947};
2014  static double b1[] = {
2015  -0.1025390625,
2016  0.2775764465332031,
2017  -1.993531733751297,
2018  2.724882731126854e1};
2019 
2020  y0 = abs(imag(z));
2021  a0 = abs(z);
2022  nm = n;
2023  if (a0 < 1.0e-100) {
2024  for (k=0;k<=n;k++) {
2025  cj[k] = czero;
2026  cy[k] = complex<double> (-1e308,0);
2027  cjp[k] = czero;
2028  cyp[k] = complex<double>(1e308,0);
2029  }
2030  cj[0] = cone;
2031  cjp[1] = complex<double>(0.5,0.0);
2032  return 0;
2033  }
2034  if ((a0 <= 300.0) || (n > (int)(0.25*a0))) {
2035  if (n == 0) nm = 1;
2036  m = msta1(a0,200);
2037  if (m < nm) nm = m;
2038  else m = msta2(a0,nm,15);
2039  cbs = czero;
2040  csu = czero;
2041  csv = czero;
2042  cf2 = czero;
2043  cf1 = complex<double> (1.0e-100,0.0);
2044  for (k=m;k>=0;k--) {
2045  cf = 2.0*(k+1.0)*cf1/z-cf2;
2046  if (k <= nm) cj[k] = cf;
2047  if (((k & 1) == 0) && (k != 0)) {
2048  if (y0 <= 1.0) {
2049  cbs += 2.0*cf;
2050  }
2051  else {
2052  cbs += (-1)*((k & 2)-1)*2.0*cf;
2053  }
2054  csu += (double)((-1)*((k & 2)-1))*cf/(double)k;
2055  }
2056  else if (k > 1) {
2057  csv += (double)((-1)*((k & 2)-1)*k)*cf/(double)(k*k-1.0);
2058  }
2059  cf2 = cf1;
2060  cf1 = cf;
2061  }
2062  if (y0 <= 1.0) cs0 = cbs+cf;
2063  else cs0 = (cbs+cf)/cos(z);
2064  for (k=0;k<=nm;k++) {
2065  cj[k] /= cs0;
2066  }
2067  ce = log(0.5*z)+el;
2068  cy[0] = M_2_PI*(ce*cj[0]-4.0*csu/cs0);
2069  cy[1] = M_2_PI*(-cj[0]/z+(ce-1.0)*cj[1]-4.0*csv/cs0);
2070  }
2071  else {
2072  ct1 = z-M_PI_4;
2073  cp0 = cone;
2074  for (k=0;k<4;k++) {
2075  cp0 += a[k]*pow(z,-2.0*k-2.0);
2076  }
2077  cq0 = -0.125/z;
2078  for (k=0;k<4;k++) {
2079  cq0 += b[k] *pow(z,-2.0*k-3.0);
2080  }
2081  cu = sqrt(M_2_PI/z);
2082  cbj0 = cu*(cp0*cos(ct1)-cq0*sin(ct1));
2083  cby0 = cu*(cp0*sin(ct1)+cq0*cos(ct1));
2084  cj[0] = cbj0;
2085  cy[0] = cby0;
2086  ct2 = z-0.75*M_PI;
2087  cp1 = cone;
2088  for (k=0;k<4;k++) {
2089  cp1 += a1[k]*pow(z,-2.0*k-2.0);
2090  }
2091  cq1 = 0.375/z;
2092  for (k=0;k<4;k++) {
2093  cq1 += b1[k]*pow(z,-2.0*k-3.0);
2094  }
2095  cbj1 = cu*(cp1*cos(ct2)-cq1*sin(ct2));
2096  cby1 = cu*(cp1*sin(ct2)+cq1*sin(ct2));
2097  cj[1] = cbj1;
2098  cy[1] = cby1;
2099  for (k=2;k<=n;k++) {
2100  cbjk = 2.0*(k-1.0)*cbj1/z-cbj0;
2101  cj[k] = cbjk;
2102  cbj0 = cbj1;
2103  cbj1 = cbjk;
2104  }
2105  }
2106  cjp[0] = -cj[1];
2107  for (k=1;k<=nm;k++) {
2108  cjp[k] = cj[k-1]-(double)k*cj[k]/z;
2109  }
2110  if (abs(cj[0]) > 1.0)
2111  cy[1] = (cj[1]*cy[0]-2.0/(M_PI*z))/cj[0];
2112  for (k=2;k<=nm;k++) {
2113  if (abs(cj[k-1]) >= abs(cj[k-2]))
2114  cyy = (cj[k]*cy[k-1]-2.0/(M_PI*z))/cj[k-1];
2115  else
2116  cyy = (cj[k]*cy[k-2]-4.0*(k-1.0)/(M_PI*z*z))/cj[k-2];
2117  cy[k] = cyy;
2118  }
2119  cyp[0] = -cy[1];
2120  for (k=1;k<=nm;k++) {
2121  cyp[k] = cy[k-1]-(double)k*cy[k]/z;
2122  }
2123 
2124  return 0;
2125 }

References a, abs(), b, cone(), cos(), czero(), el, imag(), k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), sin(), and sqrt().

◆ cbessjyva()

int CRBond_Bessel::cbessjyva ( double  v,
complex< double >  z,
double &  vm,
complex< double > *  cjv,
complex< double > *  cyv,
complex< double > *  cjvp,
complex< double > *  cyvp 
)
2129 {
2130  complex<double> z1,z2,zk,cjvl,cr,ca,cjv0,cjv1,cpz,crp;
2131  complex<double> cqz,crq,ca0,cck,csk,cyv0,cyv1,cju0,cju1,cb;
2132  complex<double> cs,cs0,cr0,cs1,cr1,cec,cf,cf0,cf1,cf2;
2133  complex<double> cfac0,cfac1,cg0,cg1,cyk,cp11,cp12,cp21,cp22;
2134  complex<double> ch0,ch1,ch2,cyl1,cyl2,cylk;
2135 
2136  double a0,v0,pv0,pv1,vl,ga,gb,vg,vv,w0,w1,ya0,yak,ya1,wa;
2137  int j,n,k,kz,l,lb,lb0,m;
2138 
2139  a0 = abs(z);
2140  z1 = z;
2141  z2 = z*z;
2142  n = (int)v;
2143 
2144 
2145  v0 = v-n;
2146 
2147  pv0 = M_PI*v0;
2148  pv1 = M_PI*(1.0+v0);
2149  if (a0 < 1.0e-100) {
2150  for (k=0;k<=n;k++) {
2151  cjv[k] = czero;
2152  cyv[k] = complex<double> (-1e308,0);
2153  cjvp[k] = czero;
2154  cyvp[k] = complex<double> (1e308,0);
2155 
2156  }
2157  if (v0 == 0.0) {
2158  cjv[0] = cone;
2159  cjvp[1] = complex<double> (0.5,0.0);
2160  }
2161  else {
2162  cjv[0] = complex<double> (1e308,0);
2163  }
2164  vm = v;
2165  return 0;
2166  }
2167  if (real(z1) < 0.0) z1 = -z;
2168  if (a0 <= 12.0) {
2169  for (l=0;l<2;l++) {
2170  vl = v0+l;
2171  cjvl = cone;
2172  cr = cone;
2173  for (k=1;k<=40;k++) {
2174  cr *= -0.25*z2/(k*(k+vl));
2175  cjvl += cr;
2176  if (abs(cr) < abs(cjvl)*eps) break;
2177  }
2178  vg = 1.0 + vl;
2179  ga = gamma(vg);
2180  ca = pow(0.5*z1,vl)/ga;
2181  if (l == 0) cjv0 = cjvl*ca;
2182  else cjv1 = cjvl*ca;
2183  }
2184  }
2185  else {
2186  if (a0 >= 50.0) kz = 8;
2187  else if (a0 >= 35.0) kz = 10;
2188  else kz = 11;
2189  for (j=0;j<2;j++) {
2190  vv = 4.0*(j+v0)*(j+v0);
2191  cpz = cone;
2192  crp = cone;
2193  for (k=1;k<=kz;k++) {
2194  crp = -0.78125e-2*crp*(vv-pow(4.0*k-3.0,2.0))*
2195  (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*z2);
2196  cpz += crp;
2197  }
2198  cqz = cone;
2199  crq = cone;
2200  for (k=1;k<=kz;k++) {
2201  crq = -0.78125e-2*crq*(vv-pow(4.0*k-1.0,2.0))*
2202  (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*z2);
2203  cqz += crq;
2204  }
2205  cqz *= 0.125*(vv-1.0)/z1;
2206  zk = z1-(0.5*(j+v0)+0.25)*M_PI;
2207  ca0 = sqrt(M_2_PI/z1);
2208  cck = cos(zk);
2209  csk = sin(zk);
2210  if (j == 0) {
2211  cjv0 = ca0*(cpz*cck-cqz*csk);
2212  cyv0 = ca0*(cpz*csk+cqz+cck);
2213  }
2214  else {
2215  cjv1 = ca0*(cpz*cck-cqz*csk);
2216  cyv1 = ca0*(cpz*csk+cqz*cck);
2217  }
2218  }
2219  }
2220  if (a0 <= 12.0) {
2221  if (v0 != 0.0) {
2222  for (l=0;l<2;l++) {
2223  vl = v0+l;
2224  cjvl = cone;
2225  cr = cone;
2226  for (k=1;k<=40;k++) {
2227  cr *= -0.25*z2/(k*(k-vl));
2228  cjvl += cr;
2229  if (abs(cr) < abs(cjvl)*eps) break;
2230  }
2231  vg = 1.0-vl;
2232  gb = gamma(vg);
2233  cb = pow(2.0/z1,vl)/gb;
2234  if (l == 0) cju0 = cjvl*cb;
2235  else cju1 = cjvl*cb;
2236  }
2237  cyv0 = (cjv0*cos(pv0)-cju0)/sin(pv0);
2238  cyv1 = (cjv1*cos(pv1)-cju1)/sin(pv1);
2239  }
2240  else {
2241  cec = log(0.5*z1)+el;
2242  cs0 = czero;
2243  w0 = 0.0;
2244  cr0 = cone;
2245  for (k=1;k<=30;k++) {
2246  w0 += 1.0/k;
2247  cr0 *= -0.25*z2/(double)(k*k);
2248  cs0 += cr0*w0;
2249  }
2250  cyv0 = M_2_PI*(cec*cjv0-cs0);
2251  cs1 = cone;
2252  w1 = 0.0;
2253  cr1 = cone;
2254  for (k=1;k<=30;k++) {
2255  w1 += 1.0/k;
2256  cr1 *= -0.25*z2/(k*(k+1.0));
2257  cs1 += cr1*(2.0*w1+1.0/(k+1.0));
2258  }
2259  cyv1 = M_2_PI*(cec*cjv1-1.0/z1-0.25*z1*cs1);
2260  }
2261  }
2262  if (real(z) < 0.0) {
2263  cfac0 = exp(pv0*cii);
2264  cfac1 = exp(pv1*cii);
2265  if (imag(z) < 0.0) {
2266  cyv0 = cfac0*cyv0-2.0*cii*cos(pv0)*cjv0;
2267  cyv1 = cfac1*cyv1-2.0*cii*cos(pv1)*cjv1;
2268  cjv0 /= cfac0;
2269  cjv1 /= cfac1;
2270  }
2271  else if (imag(z) > 0.0) {
2272  cyv0 = cyv0/cfac0+2.0*cii*cos(pv0)*cjv0;
2273  cyv1 = cyv1/cfac1+2.0*cii*cos(pv1)*cjv1;
2274  cjv0 *= cfac0;
2275  cjv1 *= cfac1;
2276  }
2277  }
2278  cjv[0] = cjv0;
2279  cjv[1] = cjv1;
2280  if ((n >= 2) && (n <= (int)(0.25*a0))) {
2281  cf0 = cjv0;
2282  cf1 = cjv1;
2283  for (k=2;k<= n;k++) {
2284  cf = 2.0*(k+v0-1.0)*cf1/z-cf0;
2285  cjv[k] = cf;
2286  cf0 = cf1;
2287  cf1 = cf;
2288  }
2289  }
2290  else if (n >= 2) {
2291  m = msta1(a0,200);
2292  if (m < n) n = m;
2293  else m = msta2(a0,n,15);
2294  cf2 = czero;
2295  cf1 = complex<double>(1.0e-100,0.0);
2296  for (k=m;k>=0;k--) {
2297  cf = 2.0*(v0+k+1.0)*cf1/z-cf2;
2298  if (k <= n) cjv[k] = cf;
2299  cf2 = cf1;
2300  cf1 = cf;
2301  }
2302  if (abs(cjv0) > abs(cjv1)) cs = cjv0/cf;
2303  else cs = cjv1/cf2;
2304  for (k=0;k<=n;k++) {
2305  cjv[k] *= cs;
2306  }
2307  }
2308  cjvp[0] = v0*cjv[0]/z-cjv[1];
2309  for (k=1;k<=n;k++) {
2310  cjvp[k] = -(k+v0)*cjv[k]/z+cjv[k-1];
2311  }
2312  cyv[0] = cyv0;
2313  cyv[1] = cyv1;
2314  ya0 = abs(cyv0);
2315  lb = 0;
2316  cg0 = cyv0;
2317  cg1 = cyv1;
2318  for (k=2;k<=n;k++) {
2319  cyk = 2.0*(v0+k-1.0)*cg1/z-cg0;
2320  yak = abs(cyk);
2321  ya1 = abs(cg0);
2322  if ((yak < ya0) && (yak< ya1)) lb = k;
2323  cyv[k] = cyk;
2324  cg0 = cg1;
2325  cg1 = cyk;
2326  }
2327  lb0 = 0;
2328  if ((lb > 4) && (imag(z) != 0.0)) {
2329  while(lb != lb0) {
2330  ch2 = cone;
2331  ch1 = czero;
2332  lb0 = lb;
2333  for (k=lb;k>=1;k--) {
2334  ch0 = 2.0*(k+v0)*ch1/z-ch2;
2335  ch2 = ch1;
2336  ch1 = ch0;
2337  }
2338  cp12 = ch0;
2339  cp22 = ch2;
2340  ch2 = czero;
2341  ch1 = cone;
2342  for (k=lb;k>=1;k--) {
2343  ch0 = 2.0*(k+v0)*ch1/z-ch2;
2344  ch2 = ch1;
2345  ch1 = ch0;
2346  }
2347  cp11 = ch0;
2348  cp21 = ch2;
2349  if (lb == n)
2350  cjv[lb+1] = 2.0*(lb+v0)*cjv[lb]/z-cjv[lb-1];
2351  if (abs(cjv[0]) > abs(cjv[1])) {
2352  cyv[lb+1] = (cjv[lb+1]*cyv0-2.0*cp11/(M_PI*z))/cjv[0];
2353  cyv[lb] = (cjv[lb]*cyv0+2.0*cp12/(M_PI*z))/cjv[0];
2354  }
2355  else {
2356  cyv[lb+1] = (cjv[lb+1]*cyv1-2.0*cp21/(M_PI*z))/cjv[1];
2357  cyv[lb] = (cjv[lb]*cyv1+2.0*cp22/(M_PI*z))/cjv[1];
2358  }
2359  cyl2 = cyv[lb+1];
2360  cyl1 = cyv[lb];
2361  for (k=lb-1;k>=0;k--) {
2362  cylk = 2.0*(k+v0+1.0)*cyl1/z-cyl2;
2363  cyv[k] = cylk;
2364  cyl2 = cyl1;
2365  cyl1 = cylk;
2366  }
2367  cyl1 = cyv[lb];
2368  cyl2 = cyv[lb+1];
2369  for (k=lb+1;k<n;k++) {
2370  cylk = 2.0*(k+v0)*cyl2/z-cyl1;
2371  cyv[k+1] = cylk;
2372  cyl1 = cyl2;
2373  cyl2 = cylk;
2374  }
2375  for (k=2;k<=n;k++) {
2376  wa = abs(cyv[k]);
2377  if (wa < abs(cyv[k-1])) lb = k;
2378  }
2379  }
2380  }
2381  cyvp[0] = v0*cyv[0]/z-cyv[1];
2382  for (k=1;k<=n;k++) {
2383  cyvp[k] = cyv[k-1]-(k+v0)*cyv[k]/z;
2384  }
2385  vm = n+v0;
2386  return 0;
2387 }

References abs(), cii(), cone(), cos(), czero(), el, eps, Eigen::bfloat16_impl::exp(), gamma(), imag(), int(), j, k, Eigen::bfloat16_impl::log(), m, M_PI, msta1(), msta2(), n, Eigen::bfloat16_impl::pow(), sin(), sqrt(), and v.

◆ cii()

static complex<double> CRBond_Bessel::cii ( 0.  0,
1.  0 
)
static

Referenced by cbessik01(), cbessiknb(), cbessikv(), cbessjy01(), and cbessjyva().

◆ cone()

static complex<double> CRBond_Bessel::cone ( 1.  0,
0.  0 
)
static

Referenced by cbessik01(), cbessiknb(), cbessikv(), cbessjy01(), cbessjyna(), cbessjynb(), cbessjyva(), and Silo::setupInitialConditions().

◆ czero()

static complex<double> CRBond_Bessel::czero ( 0.  0,
0.  0 
)
static

Referenced by cbessik01(), cbessikna(), cbessiknb(), cbessikv(), cbessjy01(), cbessjyna(), cbessjynb(), and cbessjyva().

◆ gamma()

double CRBond_Bessel::gamma ( double  x)
2403 {
2404  int i,k,m;
2405  double ga,gr,r,z;
2406 
2407  static double g[] = {
2408  1.0,
2409  0.5772156649015329,
2410  -0.6558780715202538,
2411  -0.420026350340952e-1,
2412  0.1665386113822915,
2413  -0.421977345555443e-1,
2414  -0.9621971527877e-2,
2415  0.7218943246663e-2,
2416  -0.11651675918591e-2,
2417  -0.2152416741149e-3,
2418  0.1280502823882e-3,
2419  -0.201348547807e-4,
2420  -0.12504934821e-5,
2421  0.1133027232e-5,
2422  -0.2056338417e-6,
2423  0.6116095e-8,
2424  0.50020075e-8,
2425  -0.11812746e-8,
2426  0.1043427e-9,
2427  0.77823e-11,
2428  -0.36968e-11,
2429  0.51e-12,
2430  -0.206e-13,
2431  -0.54e-14,
2432  0.14e-14};
2433 
2434  if (x > 171.0) return 1e308; // This value is an overflow flag.
2435  if (x == (int)x) {
2436  if (x > 0.0) {
2437  ga = 1.0; // use factorial
2438  for (i=2;i<x;i++) {
2439  ga *= i;
2440  }
2441  }
2442  else
2443  ga = 1e308;
2444  }
2445  else {
2446  if (fabs(x) > 1.0) {
2447  z = fabs(x);
2448  m = (int)z;
2449  r = 1.0;
2450  for (k=1;k<=m;k++) {
2451  r *= (z-k);
2452  }
2453  z -= m;
2454  }
2455  else
2456  z = x;
2457  gr = g[24];
2458  for (k=23;k>=0;k--) {
2459  gr = gr*z+g[k];
2460  }
2461  ga = 1.0/(gr*z);
2462  if (fabs(x) > 1.0) {
2463  ga *= r;
2464  if (x < 0.0) {
2465  ga = -M_PI/(x*ga*sin(M_PI*x));
2466  }
2467  }
2468  }
2469  return ga;
2470 }

References boost::multiprecision::fabs(), i, int(), k, m, M_PI, UniformPSDSelfTest::r, sin(), and plotDoE::x.

Referenced by bessikv(), bessjyv(), cbessikv(), and cbessjyva().

◆ msta1()

int CRBond_Bessel::msta1 ( double  x,
int  mp 
)
780 {
781  double a0,f0,f1,f;
782  int i,n0,n1,nn;
783 
784  a0 = fabs(x);
785  n0 = (int)(1.1*a0)+1;
786  f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-mp;
787  n1 = n0+5;
788  f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp;
789  for (i=0;i<20;i++) {
790  nn = n1-(n1-n0)/(1.0-f0/f1);
791  f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp;
792  int aux=int(std::abs(float(nn-n1)));
793  //if (int(std::abs(nn-n1)) < 1) break;
794  if (aux < 1) break;
795  n0 = n1;
796  f0 = f1;
797  n1 = nn;
798  f1 = f;
799  }
800  return nn;
801 }
Eigen::bfloat16_impl::log10
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 log10(const bfloat16 &a)
Definition: BFloat16.h:620

References abs(), f(), Global_parameters::f1(), boost::multiprecision::fabs(), i, int(), Eigen::bfloat16_impl::log10(), and plotDoE::x.

Referenced by bessikna(), bessiknb(), bessikv(), bessjyna(), bessjynb(), bessjyv(), cbessikna(), cbessiknb(), cbessikv(), cbessjyna(), cbessjynb(), and cbessjyva().

◆ msta2()

int CRBond_Bessel::msta2 ( double  x,
int  n,
int  mp 
)
803 {
804  double a0,ejn,hmp,f0,f1,f,obj;
805  int i,n0,n1,nn;
806 
807  a0 = fabs(x);
808  hmp = 0.5*mp;
809  ejn = 0.5*log10(6.28*n)-n*log10(1.36*a0/n);
810  if (ejn <= hmp) {
811  obj = mp;
812  n0 = (int)(1.1*a0);
813  if (n0 < 1) n0 = 1;
814  }
815  else {
816  obj = hmp+ejn;
817  n0 = n;
818  }
819  f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-obj;
820  n1 = n0+5;
821  f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj;
822  for (i=0;i<20;i++) {
823  nn = n1-(n1-n0)/(1.0-f0/f1);
824  f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj;
825  int aux=int(std::abs(float(nn-n1)));
826  if (aux < 1) break;
827  n0 = n1;
828  f0 = f1;
829  n1 = nn;
830  f1 = f;
831  }
832  return nn+10;
833 }

References abs(), f(), Global_parameters::f1(), boost::multiprecision::fabs(), i, int(), Eigen::bfloat16_impl::log10(), n, and plotDoE::x.

Referenced by bessikna(), bessiknb(), bessikv(), bessjyna(), bessjynb(), bessjyv(), cbessikna(), cbessiknb(), cbessikv(), cbessjyna(), cbessjynb(), and cbessjyva().

Variable Documentation

◆ el

double CRBond_Bessel::el =0.5772156649015329

Referenced by bessik01a(), bessiknb(), bessikv(), bessjy01a(), bessjyv(), cbessik01(), cbessiknb(), cbessikv(), cbessjy01(), cbessjynb(), cbessjyva(), and Eigen::internal::AmbiVector< Scalar_, StorageIndex_ >::coeffRef().

◆ eps

double CRBond_Bessel::eps =1.0e-15

Referenced by Membrane::addVertex(), array_special_functions(), bessik01a(), bessikna(), bessiknb(), bessikv(), Eigen::internal::blueNorm_impl(), cbessik01(), cbessikv(), cbessjy01(), cbessjyva(), check_product(), Eigen::internal::chkder(), NurbsUtils::close(), Eigen::NumericalDiff< Functor_, mode >::df(), Eigen::EigenSolver< MatrixType_ >::doComputeEigenvectors(), oomph::HerschelBulkleyMenDutRegConstitutiveEquation< DIM >::dviscosity_dinvariant(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::EIGEN_STATIC_ASSERT(), Eigen::internal::fdjac1(), MembraneDemo::fixMembraneEdges(), MembraneSelfTest::fixMembraneEdges(), MeshTriangle::isInsideTriangle(), LawinenBox::LawinenBox(), TwoDDGProblem< ELEMENT >::limit(), Membrane::loadFromSTL(), pblueNorm(), quaternion(), SilbertPeriodic::set_study(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficientNoDispt(), LinearViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SPHNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearPlasticViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterLinNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), SinterNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), BaseCG::setEps(), LinearPlasticViscoelasticNormalSpecies::setRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setRestitutionCoefficient(), SinterLinNormalSpecies::setRestitutionCoefficient(), SPHNormalSpecies::setRestitutionCoefficient(), LinearPlasticViscoelasticNormalSpecies::setStiffnessAndRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setStiffnessAndRestitutionCoefficient(), SinterLinNormalSpecies::setStiffnessAndRestitutionCoefficient(), SinterNormalSpecies::setStiffnessAndRestitutionCoefficient(), SPHNormalSpecies::setStiffnessAndRestitutionCoefficient(), sparse_basic(), sparse_extra(), sparse_vector(), transformations_computed_scaling_continuity(), and oomph::HerschelBulkleyMenDutRegConstitutiveEquation< DIM >::viscosity().