68 lines
1.4 KiB
C++
68 lines
1.4 KiB
C++
/*
|
||
* © Copyright 2014 – Urheberrechtshinweis
|
||
* Alle Rechte vorbehalten / All Rights Reserved
|
||
*
|
||
* Programmcode ist urheberrechtlich geschuetzt.
|
||
* Das Urheberrecht liegt, soweit nicht ausdruecklich anders gekennzeichnet, bei Frank Ebner.
|
||
* Keine Verwendung ohne explizite Genehmigung.
|
||
* (vgl. § 106 ff UrhG / § 97 UrhG)
|
||
*/
|
||
|
||
#ifndef BESSEL_H
|
||
#define BESSEL_H
|
||
|
||
namespace Distribution {
|
||
|
||
/** helper class */
|
||
template <typename T> class Bessel {
|
||
|
||
public:
|
||
|
||
/**
|
||
* calculation for I_v(z)
|
||
* https://en.wikipedia.org/wiki/Bessel_function
|
||
* http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html
|
||
*/
|
||
T getModified(const int v, const T z) const {
|
||
|
||
// running factorials (updated within for-loops)
|
||
uint64_t runFac1 = 1;
|
||
uint64_t runFac2 = factorial(v); // Delta(k+v+1) = (k+v+1-1)! = (k+v)! [k starts at 0] -> (v)!
|
||
|
||
// running potential
|
||
T runPot = 1;
|
||
const T pot = T(0.25) * z * z;
|
||
|
||
// start for-loop at k=1 instead of k=0. allows for running factorial without using if (k==0)
|
||
T sum = 0;
|
||
for (int k = 0; k < 16; ) {
|
||
sum += runPot / runFac1 / runFac2;
|
||
++k;
|
||
runFac1 *= k;
|
||
runFac2 *= k;
|
||
runPot *= pot;
|
||
}
|
||
|
||
// done
|
||
return std::pow( (T(0.5) * z), v) * sum;
|
||
|
||
}
|
||
|
||
private:
|
||
|
||
|
||
|
||
|
||
|
||
static uint64_t factorial(const int n) {
|
||
uint64_t res = 1;
|
||
for (int i = 2; i <= n; ++i) {res *= i;}
|
||
return res;
|
||
}
|
||
|
||
};
|
||
|
||
}
|
||
|
||
#endif // BESSEL_H
|