Skip to content
Snippets Groups Projects
Commit 281566a2 authored by Van Herck, Walter's avatar Van Herck, Walter
Browse files

Test code beautyfier

parent 0b695765
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Core/Tools/src/MathFunctions.cpp
+ 122
− 144
View file @ 281566a2
...@@ -24,24 +24,24 @@ ...@@ -24,24 +24,24 @@
double MathFunctions::Gaussian(double value, double average, double std_dev) double MathFunctions::Gaussian(double value, double average, double std_dev)
{ {
return StandardNormal((value-average)/std_dev) / std_dev; return StandardNormal((value - average) / std_dev) / std_dev;
} }
double MathFunctions::IntegratedGaussian(double value, double average, double std_dev) double MathFunctions::IntegratedGaussian(double value, double average, double std_dev)
{ {
double normalized_value = (value - average)/std_dev; double normalized_value = (value - average) / std_dev;
double root2 = std::sqrt(2.0); double root2 = std::sqrt(2.0);
return (gsl_sf_erf(normalized_value/root2) + 1.0)/2.0; return (gsl_sf_erf(normalized_value / root2) + 1.0) / 2.0;
} }
double MathFunctions::StandardNormal(double value) double MathFunctions::StandardNormal(double value)
{ {
return std::exp(- value*value / 2.0) / std::sqrt(2*M_PI); return std::exp(-value * value / 2.0) / std::sqrt(2 * M_PI);
} }
double MathFunctions::GenerateStandardNormalRandom() // using GSL double MathFunctions::GenerateStandardNormalRandom() // using GSL
{ {
gsl_rng * r; gsl_rng *r;
r = gsl_rng_alloc(gsl_rng_ranlxs2); r = gsl_rng_alloc(gsl_rng_ranlxs2);
double result = gsl_ran_ugaussian(r); double result = gsl_ran_ugaussian(r);
gsl_rng_free(r); gsl_rng_free(r);
...@@ -51,42 +51,43 @@ double MathFunctions::GenerateStandardNormalRandom() // using GSL ...@@ -51,42 +51,43 @@ double MathFunctions::GenerateStandardNormalRandom() // using GSL
//! @brief simple (and unoptimized) wrapper function //! @brief simple (and unoptimized) wrapper function
//! for the discrete fast Fourier transformation library (fftw3) //! for the discrete fast Fourier transformation library (fftw3)
std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<complex_t >& data, MathFunctions::EFFTDirection ftCase) std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<complex_t> &data,
MathFunctions::EFFTDirection ftCase)
{ {
double scale(1.); double scale(1.);
size_t npx = data.size(); size_t npx = data.size();
fftw_complex *ftData = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npx); fftw_complex *ftData = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
fftw_complex *ftResult = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npx); fftw_complex *ftResult = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * npx);
memset(ftData, 0, sizeof(fftw_complex) * npx); memset(ftData, 0, sizeof(fftw_complex) * npx);
memset(ftResult, 0, sizeof(fftw_complex) * npx); memset(ftResult, 0, sizeof(fftw_complex) * npx);
for(size_t i=0; i<npx; i++){ for (size_t i = 0; i < npx; i++) {
ftData[i][0] = data[i].real(); ftData[i][0] = data[i].real();
ftData[i][1] = data[i].imag(); ftData[i][1] = data[i].imag();
} }
fftw_plan plan; fftw_plan plan;
switch (ftCase) switch (ftCase) {
{ case MathFunctions::FORWARD_FFT:
case MathFunctions::FORWARD_FFT: plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_FORWARD, FFTW_ESTIMATE);
plan = fftw_plan_dft_1d( (int)npx, ftData, ftResult, FFTW_FORWARD, FFTW_ESTIMATE ); break;
break; case MathFunctions::BACKWARD_FFT:
case MathFunctions::BACKWARD_FFT: plan = fftw_plan_dft_1d((int)npx, ftData, ftResult, FFTW_BACKWARD, FFTW_ESTIMATE);
plan = fftw_plan_dft_1d( (int)npx, ftData, ftResult, FFTW_BACKWARD, FFTW_ESTIMATE ); scale = 1. / double(npx);
scale = 1./double(npx); break;
break; default:
default: throw std::runtime_error(
throw std::runtime_error("MathFunctions::FastFourierTransform -> Panic! Unknown transform case."); "MathFunctions::FastFourierTransform -> Panic! Unknown transform case.");
} }
fftw_execute(plan); fftw_execute(plan);
// saving data for user // saving data for user
std::vector<complex_t > outData; std::vector<complex_t> outData;
outData.resize(npx); outData.resize(npx);
for(size_t i=0; i<npx; i++){ for (size_t i = 0; i < npx; i++) {
outData[i] = scale*complex_t(ftResult[i][0], ftResult[i][1]); outData[i] = scale * complex_t(ftResult[i][0], ftResult[i][1]);
} }
fftw_destroy_plan(plan); fftw_destroy_plan(plan);
...@@ -100,36 +101,42 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<co ...@@ -100,36 +101,42 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<co
//! for the discrete fast Fourier transformation library (fftw3); //! for the discrete fast Fourier transformation library (fftw3);
//! transforms real to complex //! transforms real to complex
std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<double >& data, MathFunctions::EFFTDirection ftCase) std::vector<complex_t> MathFunctions::FastFourierTransform(const std::vector<double> &data,
MathFunctions::EFFTDirection ftCase)
{ {
std::vector<complex_t> cdata; std::vector<complex_t> cdata;
cdata.resize(data.size()); cdata.resize(data.size());
for(size_t i=0; i<data.size(); i++) { cdata[i] = complex_t(data[i], 0); } for (size_t i = 0; i < data.size(); i++) {
cdata[i] = complex_t(data[i], 0);
}
return MathFunctions::FastFourierTransform(cdata, ftCase); return MathFunctions::FastFourierTransform(cdata, ftCase);
} }
//! convolution of two real vectors of equal size //! convolution of two real vectors of equal size
std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double>& signal, const std::vector<double>& resfunc) std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double> &signal,
const std::vector<double> &resfunc)
{ {
if(signal.size() != resfunc.size() ) { if (signal.size() != resfunc.size()) {
throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for two vectors of equal size. Use Convolve class instead."); throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for "
"two vectors of equal size. Use Convolve class instead.");
} }
std::vector<complex_t > fft_signal = MathFunctions::FastFourierTransform(signal, MathFunctions::FORWARD_FFT); std::vector<complex_t> fft_signal
std::vector<complex_t > fft_resfunc = MathFunctions::FastFourierTransform(resfunc, MathFunctions::FORWARD_FFT); = MathFunctions::FastFourierTransform(signal, MathFunctions::FORWARD_FFT);
std::vector<complex_t> fft_resfunc
= MathFunctions::FastFourierTransform(resfunc, MathFunctions::FORWARD_FFT);
std::vector<complex_t > fft_prod; std::vector<complex_t> fft_prod;
fft_prod.resize(fft_signal.size()); fft_prod.resize(fft_signal.size());
for(size_t i=0; i<fft_signal.size(); i++){ for (size_t i = 0; i < fft_signal.size(); i++) {
fft_prod[i] = fft_signal[i] * fft_resfunc[i]; fft_prod[i] = fft_signal[i] * fft_resfunc[i];
} }
std::vector<complex_t > result = MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BACKWARD_FFT); std::vector<complex_t> result
= MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BACKWARD_FFT);
return result; return result;
} }
//! Complex Bessel function of 1st kind, 0st order //! Complex Bessel function of 1st kind, 0st order
//! //!
//! Taken from http://www.crbond.com/math.htm //! Taken from http://www.crbond.com/math.htm
...@@ -140,73 +147,58 @@ std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double>& sig ...@@ -140,73 +147,58 @@ std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double>& sig
complex_t MathFunctions::crbond_bessel_J0(const complex_t &z) complex_t MathFunctions::crbond_bessel_J0(const complex_t &z)
{ {
complex_t cj0; complex_t cj0;
static const complex_t cone(1.0,0.0); static const complex_t cone(1.0, 0.0);
static const complex_t czero(0.0,0.0); static const complex_t czero(0.0, 0.0);
static const double eps=1e-15; static const double eps = 1e-15;
static double a[] = { static double a[] = { -7.03125e-2, 0.112152099609375, -0.5725014209747314,
-7.03125e-2, 6.074042001273483, -1.100171402692467e2, 3.038090510922384e3,
0.112152099609375, -1.188384262567832e5, 6.252951493434797e6, -4.259392165047669e8,
-0.5725014209747314, 3.646840080706556e10, -3.833534661393944e12, 4.854014686852901e14,
6.074042001273483, -7.286857349377656e16, 1.279721941975975e19 };
-1.100171402692467e2, static double b[] = { 7.32421875e-2, -0.2271080017089844, 1.727727502584457,
3.038090510922384e3, -2.438052969955606e1, 5.513358961220206e2, -1.825775547429318e4,
-1.188384262567832e5, 8.328593040162893e5, -5.006958953198893e7, 3.836255180230433e9,
6.252951493434797e6, -3.649010818849833e11, 4.218971570284096e13, -5.827244631566907e15,
-4.259392165047669e8, 9.476288099260110e17, -1.792162323051699e20 };
3.646840080706556e10,
-3.833534661393944e12,
4.854014686852901e14,
-7.286857349377656e16,
1.279721941975975e19};
static double b[] = {
7.32421875e-2,
-0.2271080017089844,
1.727727502584457,
-2.438052969955606e1,
5.513358961220206e2,
-1.825775547429318e4,
8.328593040162893e5,
-5.006958953198893e7,
3.836255180230433e9,
-3.649010818849833e11,
4.218971570284096e13,
-5.827244631566907e15,
9.476288099260110e17,
-1.792162323051699e20};
double a0 = std::abs(z); double a0 = std::abs(z);
complex_t z1 = z; complex_t z1 = z;
if (a0 == 0.0) return cone; if (a0 == 0.0)
if (std::real(z) < 0.0) z1 = -z; return cone;
if (std::real(z) < 0.0)
z1 = -z;
if (a0 <= 12.0) { if (a0 <= 12.0) {
complex_t z2 = z*z; complex_t z2 = z * z;
cj0 = cone; cj0 = cone;
complex_t cr = cone; complex_t cr = cone;
for (size_t k=1; k<=40; ++k) { for (size_t k = 1; k <= 40; ++k) {
cr *= -0.25*z2/(double)(k*k); cr *= -0.25 * z2 / (double)(k * k);
cj0 += cr; cj0 += cr;
if (std::abs(cr) < std::abs(cj0)*eps) break; if (std::abs(cr) < std::abs(cj0) * eps)
break;
} }
} else { } else {
size_t kz; size_t kz;
if (a0 >= 50.0) kz = 8; // can be changed to 10 if (a0 >= 50.0)
else if (a0 >= 35.0) kz = 10; // " " " 12 kz = 8; // can be changed to 10
else kz = 12; // " " " 14 else if (a0 >= 35.0)
kz = 10; // " " " 12
else
kz = 12; // " " " 14
complex_t ct1 = z1 - M_PI_4; complex_t ct1 = z1 - M_PI_4;
complex_t cp0 = cone; complex_t cp0 = cone;
complex_t cq0 = -0.125/z1; complex_t cq0 = -0.125 / z1;
complex_t ptmp=std::pow(z1,-2.0); complex_t ptmp = std::pow(z1, -2.0);
for (size_t k=0; k<kz; ++k) { for (size_t k = 0; k < kz; ++k) {
cp0 += a[k]*ptmp; cp0 += a[k] * ptmp;
cq0 += b[k]*ptmp/z1; cq0 += b[k] * ptmp / z1;
ptmp /= (z1*z1); ptmp /= (z1 * z1);
} }
cj0 = std::sqrt(M_2_PI/z1)*(cp0*std::cos(ct1)-cq0*std::sin(ct1)); cj0 = std::sqrt(M_2_PI / z1) * (cp0 * std::cos(ct1) - cq0 * std::sin(ct1));
} }
return cj0; return cj0;
} }
//! Complex Bessel function of 1st kind, 1st order //! Complex Bessel function of 1st kind, 1st order
//! //!
//! Taken from http://www.crbond.com/math.htm //! Taken from http://www.crbond.com/math.htm
...@@ -217,88 +209,74 @@ complex_t MathFunctions::crbond_bessel_J0(const complex_t &z) ...@@ -217,88 +209,74 @@ complex_t MathFunctions::crbond_bessel_J0(const complex_t &z)
complex_t MathFunctions::crbond_bessel_J1(const complex_t &z) complex_t MathFunctions::crbond_bessel_J1(const complex_t &z)
{ {
complex_t cj1; complex_t cj1;
static const complex_t cone(1.0,0.0); static const complex_t cone(1.0, 0.0);
static const complex_t czero(0.0,0.0); static const complex_t czero(0.0, 0.0);
static const double eps=1e-15; static const double eps = 1e-15;
static double a1[] = { static double a1[] = { 0.1171875, -0.1441955566406250, 0.6765925884246826,
0.1171875, -6.883914268109947, 1.215978918765359e2, -3.302272294480852e3,
-0.1441955566406250, 1.276412726461746e5, -6.656367718817688e6, 4.502786003050393e8,
0.6765925884246826, -3.833857520742790e10, 4.011838599133198e12, -5.060568503314727e14,
-6.883914268109947, 7.572616461117958e16, -1.326257285320556e19 };
1.215978918765359e2, static double b1[] = { -0.1025390625, 0.2775764465332031, -1.993531733751297,
-3.302272294480852e3, 2.724882731126854e1, -6.038440767050702e2, 1.971837591223663e4,
1.276412726461746e5, -8.902978767070678e5, 5.310411010968522e7, -4.043620325107754e9,
-6.656367718817688e6, 3.827011346598605e11, -4.406481417852278e13, 6.065091351222699e15,
4.502786003050393e8, -9.833883876590679e17, 1.855045211579828e20 };
-3.833857520742790e10,
4.011838599133198e12,
-5.060568503314727e14,
7.572616461117958e16,
-1.326257285320556e19};
static double b1[] = {
-0.1025390625,
0.2775764465332031,
-1.993531733751297,
2.724882731126854e1,
-6.038440767050702e2,
1.971837591223663e4,
-8.902978767070678e5,
5.310411010968522e7,
-4.043620325107754e9,
3.827011346598605e11,
-4.406481417852278e13,
6.065091351222699e15,
-9.833883876590679e17,
1.855045211579828e20};
double a0 = std::abs(z); double a0 = std::abs(z);
if (a0 == 0.0) return czero; if (a0 == 0.0)
return czero;
complex_t z1 = z; complex_t z1 = z;
if (std::real(z) < 0.0) z1 = -z; if (std::real(z) < 0.0)
z1 = -z;
if (a0 <= 12.0) { if (a0 <= 12.0) {
complex_t z2 = z*z; complex_t z2 = z * z;
cj1 = cone; cj1 = cone;
complex_t cr = cone; complex_t cr = cone;
for (size_t k=1; k<=40; ++k) { for (size_t k = 1; k <= 40; ++k) {
cr *= -0.25*z2/(k*(k+1.0)); cr *= -0.25 * z2 / (k * (k + 1.0));
cj1 += cr; cj1 += cr;
if (std::abs(cr) < std::abs(cj1)*eps) break; if (std::abs(cr) < std::abs(cj1) * eps)
break;
} }
cj1 *= 0.5*z1; cj1 *= 0.5 * z1;
} else { } else {
size_t kz; size_t kz;
if (a0 >= 50.0) kz = 8; // can be changed to 10 if (a0 >= 50.0)
else if (a0 >= 35.0) kz = 10; // " " " 12 kz = 8; // can be changed to 10
else kz = 12; // " " " 14 else if (a0 >= 35.0)
kz = 10; // " " " 12
else
kz = 12; // " " " 14
complex_t cp1 = cone; complex_t cp1 = cone;
complex_t cq1 = 0.375/z1; complex_t cq1 = 0.375 / z1;
complex_t ptmp=std::pow(z1,-2.0); complex_t ptmp = std::pow(z1, -2.0);
for (size_t k=0; k<kz; ++k) { for (size_t k = 0; k < kz; ++k) {
cp1 += a1[k]*ptmp; cp1 += a1[k] * ptmp;
cq1 += b1[k]*ptmp/z1; cq1 += b1[k] * ptmp / z1;
ptmp /= (z1*z1); ptmp /= (z1 * z1);
} }
complex_t ct2 = z1 - 0.75*M_PI; complex_t ct2 = z1 - 0.75 * M_PI;
cj1 = std::sqrt(M_2_PI/z1)*(cp1*std::cos(ct2)-cq1*std::sin(ct2)); cj1 = std::sqrt(M_2_PI / z1) * (cp1 * std::cos(ct2) - cq1 * std::sin(ct2));
} }
if (std::real(z) < 0.0) cj1 = -cj1; if (std::real(z) < 0.0)
cj1 = -cj1;
return cj1; return cj1;
} }
complex_t MathFunctions::geometricSum(complex_t z, int exponent) complex_t MathFunctions::geometricSum(complex_t z, int exponent)
{ {
if (exponent<1) { if (exponent < 1) {
throw LogicErrorException( throw LogicErrorException("MathFunctions::geometricSeries:"
"MathFunctions::geometricSeries:" " exponent should be > 0");
" exponent should be > 0");
} }
complex_t result(0.0, 0.0); complex_t result(0.0, 0.0);
double nd = (double)exponent; double nd = (double)exponent;
--exponent; --exponent;
while (exponent>0) { while (exponent > 0) {
result += std::pow(z, exponent)*(nd-exponent); result += std::pow(z, exponent) * (nd - exponent);
--exponent; --exponent;
} }
return result; return result;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment