Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
special.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12 
13 #ifndef HEX_SPECF
14 #define HEX_SPECF
15 
16 #include <cmath>
17 
18 #include <gsl/gsl_math.h>
19 #include <gsl/gsl_sf.h>
20 
22 #define fac(n) gsl_sf_fact(n)
23 
25 #define sqr(x) (gsl_sf_pow_int((x),2))
26 
27 #include "arrays.h"
28 #include "complex.h"
29 #include "misc.h"
30 
32 #define ALPHA_PLUS(n,l) (sqrt(((n)-(l))*((n)+(l)+1)))
33 #define ALPHA_MINUS(n,l) (sqrt(((n)-(l)-1)*((n)+(l))))
34 
36 #define DELTA(a,b) ((a) == (b) ? 1 : 0)
37 
38 namespace special
39 {
40 
41 namespace constant
42 {
43 
44 const double e = M_E; // e
45 
46 const double pi = M_PI; // π
47 const double pi_half = M_PI_2; // π/2
48 const double pi_quart = M_PI_4; // π/4
49 const double pi_sqrt = M_SQRTPI; // √π
50 const double pi_inv = M_1_PI; // 1/π
51 const double pi_two_inv = M_2_PI; // 2/π
52 const double pi_two_invsqrt = M_2_SQRTPI; // 2/√π
53 
54 const double sqrt_two = M_SQRT2; // √2
55 const double sqrt_three = M_SQRT3; // √3
56 const double sqrt_half = M_SQRT1_2; // 1/√2
57 
58 const double log2e = M_LOG2E; // log₂ e
59 const double log10e = M_LOG10E; // log e
60 const double ln2 = M_LN2; // ln 2
61 const double ln10 = M_LN10; // ln 10
62 const double lnpi = M_LNPI; // ln π
63 
64 const double euler = M_EULER; // γ
65 
66 // other special numbers
67 const double Inf = std::numeric_limits<double>::infinity();
68 const double Nan = std::numeric_limits<double>::quiet_NaN();
69 
70 }; // end of namespace "special::constant"
71 
72 namespace integral
73 {
74 
78 template <class T> T trapz (double h, NumberArray<T> const & y)
79 {
80  return h * (sum(y) - 0.5 * (y.front() + y.back()));
81 }
82 
86 template <class T> T trapz (NumberArray<T> const & x, NumberArray<T> const & y)
87 {
88  assert (x.size() == y.size());
89 
90  T sum = 0;
91  for (unsigned i = 1; i < x.size(); i++)
92  sum += 0.5 * (y[i] + y[i-1]) * (x[i] - x[i-1]);
93 
94  return sum;
95 }
96 
100 template <class T> T simpson (double h, NumberArray<T> const & y)
101 {
102  if (y.size() % 2 != 0)
103  throw exception ("You need to use even number of grid points for Simpson integration.");
104 
105  T sum1 = 0, sum2 = 0;
106 
107  for (int i = 1; i < y.size(); i += 2)
108  sum1 += y[i];
109  for (int i = 2; i < y.size() - 1; i += 2)
110  sum2 += y[i];
111 
112  return h * (y.front() + 4. * sum1 + 2. * sum2 + y.back()) / 3.;
113 }
114 
141 template <class T> T pow_exp_hyperg1F1 (T a, T b, T c, T u, T v, T x, double epsrel = 1e-10, unsigned maxiter = 1000)
142 {
143  // last added term and up-to-now sum of terms
144  T term = 1. / (a + 1.);
145  T bfactor = 1.;
146  T suma = term;
147 
148  // auxiliary variables
149  T c_over_b = c / b;
150 
151  // the outer sum
152  for (unsigned n = 1; std::abs(term) > epsrel * std::abs(suma); n++)
153  {
154  // compute the inner sum
155  T iterm = 1, isum = 0;
156  for (unsigned k = 0; k <= n; k++)
157  {
158  // update the inner sum
159  isum += iterm;
160 
161  // update term for the next k
162  iterm *= (u + T(k)) / (v + T(k)) * ((n - k) / (k + 1.)) * c_over_b;
163  }
164 
165  // update the term and the outer sum
166  bfactor *= b * x / T(n);
167  term = bfactor * isum / (a + T(n + 1));
168  suma += term;
169 
170  // check if we run out of allowed iterations
171  if (n == maxiter)
172  {
173  throw exception
174  (
175  "Maximal number of iterations (%d) reached in pow_exp_hyperg1F1.",
176  maxiter
177  );
178  }
179  }
180 
181  // return the result
182  return suma;
183 }
184 
188 template <class T> NumberArray<T> romberg (const ArrayView<T> y)
189 {
190  NumberArray<T> z = y;
191  unsigned N = y.size();
192 
193  T scale = 1;
194  for (unsigned log4scale = 1; log4scale < N; log4scale++)
195  {
196  scale *= 4;
197  for (unsigned j = log4scale; j < N; j++)
198  {
199  z[j] = (scale * z[j] - z[j-1]) / (scale - T(1));
200  }
201  }
202 
203  return z;
204 }
205 
206 }; // end of namespace "special::integral"
207 
219 int coulomb_zeros (double eta, int L, int nzeros, double * zeros, double epsrel = 1e-8);
220 
271 std::vector<std::vector<int>> FdB_partition (int n);
272 
291 template <class T, class OuterFunctionDerivative, class InnerFunctionDerivative> T chain_rule
292 (
293  OuterFunctionDerivative Df,
294  InnerFunctionDerivative Dg,
295  int n,
296  double x
297 )
298 {
299  // no derivative : evaluate the function
300  if (n == 0)
301  return Df(0,Dg(0,x));
302 
303  // first and higher derivative : use the Faa di Bruno formula
304  T suma = 0, term = 0;
305  for (std::vector<int> & counts : FdB_partition(n))
306  {
307  // evaluate derivative of "f"
308  if ((term = Df(std::accumulate(counts.begin(), counts.end(), 0),Dg(0,x))) == 0)
309  continue;
310 
311  // evaluate all derivatives of "g"
312  for (int j = 1; j <= n; j++)
313  {
314  term *= std::pow(Dg(j,x) / gsl_sf_fact(j), counts[j-1]) / gsl_sf_fact(counts[j-1]);
315 
316  if (term == 0)
317  break;
318  }
319 
320  // update the sum
321  suma += term;
322  }
323  return gsl_sf_fact(n) * suma;
324 }
325 
326 }; // end of namespace "special"
327 
334 Complex hydro_P(unsigned n, unsigned l, Complex z);
335 
341 Complex dhydro_P(unsigned n, unsigned l, Complex z);
342 
343 //
344 // Ricatti-Bessel functions.
345 //
346 
367 Complex ric_j(int n, Complex z);
368 
377 cArray ric_jv(int lmax, Complex z);
378 
384 Complex dric_j(int n, Complex z);
385 
394 cArray dric_jv(int lmax, Complex z);
395 
401 inline double ric_j (int n, double x)
402 {
403  return x * gsl_sf_bessel_jl(n,x);
404 }
405 
411 inline double ric_n (int n, double x)
412 {
413  return x * gsl_sf_bessel_yl(n,x);
414 }
415 
416 //
417 // Modified Bessel function of the first kind.
418 //
419 
425 inline double sph_i_scaled (int n, double x)
426 {
427  return gsl_sf_bessel_il_scaled(n,x);
428 }
429 
435 inline double ric_i_scaled (int n, double x)
436 {
437  return x * gsl_sf_bessel_il_scaled(n,x);
438 }
439 
445 inline double sph_i (int n, double x)
446 {
447  return exp(x) * sph_i_scaled(n,x);
448 }
449 
455 inline double ric_i (int n, double x)
456 {
457  return exp(x)*ric_i_scaled(n,x);
458 }
459 
460 
461 //
462 // Modified Bessel function of the second kind.
463 //
464 
470 inline double sph_k_scaled (int n, double x)
471 {
472  return gsl_sf_bessel_kl_scaled(n,x) * M_2_PI;
473 }
474 
480 inline double ric_k_scaled (int n, double x)
481 {
482  return x * sph_k_scaled(n,x);
483 }
484 
490 inline double sph_k (int n, double x)
491 {
492  return exp(-x) * sph_k_scaled(n,x);
493 }
494 
500 inline double ric_k (int n, double x)
501 {
502  return exp(-x) * ric_k_scaled(n,x);
503 }
504 
505 //
506 // Ricatti-Hankel functions.
507 //
508 
514 inline Complex ric_h_plus (int n, double x)
515 {
516  return Complex(ric_j(n,x), ric_n(n,x));
517 }
518 
519 //
520 // Other special functions.
521 //
522 
526 Complex sphY(int l, int m, double theta, double phi);
527 
531 Complex sphBiY(int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2);
532 
536 template <typename T> class RadialFunction
537 {
538 public:
539 
541  virtual T operator() (double x) const = 0;
542 };
543 
547 inline double dric_j(int n, double x)
548 {
549  if (n == 0)
550  return cos(x);
551 
552  return gsl_sf_bessel_jl(n,x) + (n * ric_j(n-1,x) - (n+1) * ric_j(n+1,x)) / (2*n+1);
553 }
554 
558 inline double dric_n(int n, double x)
559 {
560  if (n == 0)
561  return sin(x);
562 
563  return gsl_sf_bessel_yl(n,x) + (n * ric_n(n-1,x) - (n+1) * ric_n(n+1,x)) / (2*n+1);
564 }
565 
569 inline double dric_i(int n, double x)
570 {
571  if (n == 0)
572  return cosh(x);
573 
574  return exp(x) * (sph_i_scaled(n,x) + (n * ric_i_scaled(n-1,x) + (n+1) * ric_i_scaled(n+1,x)) / (2*n+1));
575 }
576 
580 inline double dric_k(int n, double x)
581 {
582  if (n == 0)
583  return -exp(-x);
584 
585  return exp(-x) * (sph_k_scaled(n,x) - (n * ric_k_scaled(n-1,x) + (n+1) * ric_k_scaled(n+1,x)) / (2*n+1));
586 }
587 
591 inline double dric_i_scaled(int n, double x)
592 {
593  if (n == 0)
594  return exp(-2*x);
595 
596  return sph_i_scaled(n,x) - ric_i_scaled(n,x) + (n * ric_i_scaled(n-1,x) + (n+1) * ric_i_scaled(n+1,x)) / (2*n+1);
597 }
598 
602 inline double dric_k_scaled(int n, double x)
603 {
604  if (n == 0)
605  return 0;
606 
607  return sph_k_scaled(n,x) + ric_k_scaled(n,x) - (n * ric_k_scaled(n-1,x) + (n+1) * ric_k_scaled(n+1,x)) / (2*n+1);
608 }
609 
613 inline Complex dric_h_plus (int n, double x)
614 {
615  return Complex(dric_j(n,x),dric_n(n,x));
616 }
617 
632 int coul_F_michel(int l, double k, double r, double& F, double& Fp);
633 
642 int coul_F (int l, double k, double r, double & F, double & Fp);
643 
651 double coul_F_asy(int l, double k, double r, double sigma = special::constant::Nan);
652 
658 double coul_F_sigma(int l, double k);
659 
669 bool makes_triangle (int two_j1, int two_j2, int two_j3);
670 
681 double logdelta (int two_j1, int two_j2, int two_j3);
682 
707 double Wigner3j_2 (int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3);
708 #define Wigner3j(a,b,c,d,e,f) Wigner3j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))
709 
710 
732 double Wigner6j_2 (int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6);
733 #define Wigner6j(a,b,c,d,e,f) Wigner6j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))
734 
735 
739 double computef (int lambda, int l1, int l2, int l1p, int l2p, int L);
740 
748 double ClebschGordan(int l1, int m1, int l2, int m2, int L, int M);
749 
758 double Gaunt(int l1, int m1, int l2, int m2, int l, int m);
759 
764 int triangle_count(int L, int maxell);
765 
780 double Hyper2F1 (double a, double b, double c, double x);
781 
782 #endif
cArray dric_jv(int lmax, Complex z)
Definition: special.cpp:175
Array view.
Definition: arrays.h:186
T simpson(double h, NumberArray< T > const &y)
Uniform Simpson integration.
Definition: special.h:100
double dric_k(int n, double x)
Derivative of the modified Ricatti-Bessel function of the second kind, .
Definition: special.h:580
Complex ric_j(int l, Complex z)
Definition: special.cpp:170
double ric_n(int n, double x)
Ricatti-Bessel function of the second kind, .
Definition: special.h:411
virtual T operator()(double x) const =0
Evaluate the function.
double ric_k(int n, double x)
Modified Ricatti-Bessel function of the second kind, .
Definition: special.h:500
Complex dric_h_plus(int n, double x)
Derivative of Ricatti-Hankel function of the first kind, .
Definition: special.h:613
double computef(int lambda, int l1, int l2, int l1p, int l2p, int L)
Definition: special.cpp:672
bool makes_triangle(int two_j1, int two_j2, int two_j3)
Check triangle inequality.
Definition: special.cpp:518
T trapz(double h, NumberArray< T > const &y)
Uniform trapezoidal integration.
Definition: special.h:78
Complex ric_h_plus(int n, double x)
Ricatti-Hankel function of the first kind, .
Definition: special.h:514
const double sqrt_half
Definition: special.h:56
T & front(int i=0)
Definition: arrays.h:777
const double ln2
Definition: special.h:60
const double lnpi
Definition: special.h:62
Complex sphBiY(int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2)
Bi-polar spherical harmonic function.
Definition: special.cpp:384
const double pi_two_inv
Definition: special.h:51
Complex dric_j(int l, Complex z)
Definition: special.cpp:200
NumberArray< T > romberg(const ArrayView< T > y)
Romberg integration.
Definition: special.h:188
NumberArray< T > cos(NumberArray< T > const &A)
Return per-element cosine.
Definition: arrays.h:1747
double Wigner6j_2(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6)
Wigner 6j coefficient.
Definition: special.cpp:609
double sph_i(int n, double x)
Modified spherical Bessel function of the first kind, .
Definition: special.h:445
A comfortable number array class.
Definition: arrays.h:171
const double ln10
Definition: special.h:61
class RadialFunction exception
T pow_exp_hyperg1F1(T a, T b, T c, T u, T v, T x, double epsrel=1e-10, unsigned maxiter=1000)
Compute integral of the confluent hypergeometric function.
Definition: special.h:141
const double sqrt_two
Definition: special.h:54
const double log2e
Definition: special.h:58
const double pi_inv
Definition: special.h:50
T & back(int i=0)
Definition: arrays.h:779
double sph_k_scaled(int n, double x)
Scaled modified spherical Bessel function of the second kind, .
Definition: special.h:470
const double e
Definition: special.h:44
int triangle_count(int L, int maxl)
Definition: special.cpp:774
const double euler
Definition: special.h:64
double dric_k_scaled(int n, double x)
Scaled derivative of modified Ricatti-Bessel function of the second kind, .
Definition: special.h:602
double dric_i_scaled(int n, double x)
Derivative of scaled modified Ricatti-Bessel function of the second kind, .
Definition: special.h:591
const double pi_sqrt
Definition: special.h:49
const double log10e
Definition: special.h:59
size_t size() const
Length of the array (number of elements).
Definition: arrays.h:276
const double Inf
Definition: special.h:67
double ClebschGordan(int __j1, int __m1, int __j2, int __m2, int __J, int __M)
Clebsch-Gordan coefficient.
Definition: special.cpp:704
int coul_F_michel(int l, double k, double r, double &F, double &Fp)
Uniform approximation to the Coulomb wave function.
Definition: special.cpp:403
T sum(const ArrayView< T > v)
Sum elements in array.
Definition: arrays.h:1770
double dric_n(int n, double x)
Derivative of Ricatti-Bessel function of the second kind, .
Definition: special.h:558
double Hyper2F1(double a, double b, double c, double x)
The hypergeometric function .
Definition: special.cpp:786
const double sqrt_three
Definition: special.h:55
Complex dhydro_P(unsigned n, unsigned l, Complex z)
Definition: special.cpp:338
NumberArray< T > sin(NumberArray< T > const &A)
Return per-element sine.
Definition: arrays.h:1723
double Gaunt(int l1, int m1, int l2, int m2, int l, int m)
Gaunt&#39;s integral.
Definition: special.cpp:755
double coul_F_asy(int l, double k, double r, double sigma)
Asymptotic form of the regular Coulomb wave.
Definition: special.cpp:511
NumberArray< T > pow(NumberArray< T > const &u, double e)
Return per-element power.
Definition: arrays.h:1683
double dric_i(int n, double x)
Derivative of the modified Ricatti-Bessel function of the second kind, .
Definition: special.h:569
std::vector< std::vector< int > > FdB_partition(int n)
Faa di Bruno partitioning.
Definition: special.cpp:823
double coul_F_sigma(int l, double k)
Coulomb phase shift.
Definition: special.cpp:498
int coulomb_zeros(double eta, int L, int nzeros, double *zeros, double epsrel=1e-8)
Get zeros of the Coulomb wave function .
Definition: special.cpp:60
double ric_i_scaled(int n, double x)
Scaled modified Ricatti-Bessel function of the first kind, .
Definition: special.h:435
double sph_i_scaled(int n, double x)
Scaled modified spherical Bessel function of the first kind, .
Definition: special.h:425
const double pi_two_invsqrt
Definition: special.h:52
T chain_rule(OuterFunctionDerivative Df, InnerFunctionDerivative Dg, int n, double x)
Chain rule for n-th derivative.
Definition: special.h:292
const double pi_quart
Definition: special.h:48
Exception class.
Definition: misc.h:39
double logdelta(int two_j1, int two_j2, int two_j3)
Auxiliary function used in coupling coefficient computations.
Definition: special.cpp:525
cArray ric_jv(int lmax, Complex z)
Definition: special.cpp:107
rArray abs(const cArrayView u)
Definition: arrays.cpp:19
std::complex< double > Complex
Definition: complex.h:20
double ric_k_scaled(int n, double x)
Scaled modified Ricatti-Bessel function of the second kind, .
Definition: special.h:480
double ric_i(int n, double x)
Modified Ricatti-Bessel function of the first kind, .
Definition: special.h:455
Base class for radial functions.
Definition: special.h:536
const double pi
Definition: special.h:46
const double Nan
Definition: special.h:68
Complex sphY(int l, int m, double theta, double phi)
Spherical harmonic function.
Definition: special.cpp:352
Complex hydro_P(unsigned n, unsigned l, Complex z)
Definition: special.cpp:319
size_t size() const
Item count.
Definition: arrays.h:673
double sph_k(int n, double x)
Modified spherical Bessel function of the second kind, .
Definition: special.h:490
int coul_F(int l, double k, double r, double &F, double &Fp)
Evaluate Coulomb wave function (and its derivative).
Definition: special.cpp:448
double Wigner3j_2(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
Wigner 3j coefficient.
Definition: special.cpp:535
const double pi_half
Definition: special.h:47