Skip to content
Snippets Groups Projects
Commit 63669d58 authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

minor edits while reading

parent 14e19782
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,8 @@ ...@@ -3,7 +3,8 @@
// BornAgain: simulate and fit scattering at grazing incidence // BornAgain: simulate and fit scattering at grazing incidence
// //
//! @file Tools/inc/MathFunctions.h //! @file Tools/inc/MathFunctions.h
//! @brief Define many functions in namespace MathFunctions. //! @brief Define many functions in namespace MathFunctions,
//! and provide inline implementation for most of them
//! //!
//! @homepage http://apps.jcns.fz-juelich.de/BornAgain //! @homepage http://apps.jcns.fz-juelich.de/BornAgain
//! @license GNU General Public License v3 or higher (see COPYING) //! @license GNU General Public License v3 or higher (see COPYING)
...@@ -28,8 +29,11 @@ ...@@ -28,8 +29,11 @@
#include "gsl/gsl_sf_expint.h" #include "gsl/gsl_sf_expint.h"
#include "gsl/gsl_integration.h" #include "gsl/gsl_integration.h"
//! Various mathematical functions.
namespace MathFunctions namespace MathFunctions
{ {
double Gaussian(double value, double average, double std_dev); double Gaussian(double value, double average, double std_dev);
double IntegratedGaussian(double value, double average, double std_dev); double IntegratedGaussian(double value, double average, double std_dev);
...@@ -116,20 +120,17 @@ inline double MathFunctions::Sinc(double value) // Sin(x)/x ...@@ -116,20 +120,17 @@ inline double MathFunctions::Sinc(double value) // Sin(x)/x
inline complex_t MathFunctions::Sinc(const complex_t &value) // Sin(x)/x inline complex_t MathFunctions::Sinc(const complex_t &value) // Sin(x)/x
{ {
if(std::abs(value)<Numeric::double_epsilon) { if(std::abs(value)<Numeric::double_epsilon)
return complex_t(1.0, 0.0); return complex_t(1.0, 0.0);
}
return std::sin(value)/value; return std::sin(value)/value;
} }
inline complex_t MathFunctions::Laue(const complex_t &value, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x) inline complex_t MathFunctions::Laue(const complex_t &value, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x)
{ {
if (N==0) { if (N==0)
return complex_t(1.0, 0.0); return complex_t(1.0, 0.0);
} if(std::abs(value)<Numeric::double_epsilon)
if(std::abs(value)<Numeric::double_epsilon) {
return complex_t(N+1.0, 0.0); return complex_t(N+1.0, 0.0);
}
return std::exp(complex_t(0.0, 1.0)*value*(double)N/2.0)*std::sin(value*(N+1.0)/2.0)/std::sin(value/2.0); return std::exp(complex_t(0.0, 1.0)*value*(double)N/2.0)*std::sin(value*(N+1.0)/2.0)/std::sin(value/2.0);
} }
......
...@@ -23,7 +23,7 @@ ...@@ -23,7 +23,7 @@
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)
...@@ -45,11 +45,9 @@ double MathFunctions::GenerateStandardNormalRandom() // using Box-Muller transf ...@@ -45,11 +45,9 @@ double MathFunctions::GenerateStandardNormalRandom() // using Box-Muller transf
return std::sqrt(-2.0*std::log(u1))*std::cos(2*M_PI*u2); return std::sqrt(-2.0*std::log(u1))*std::cos(2*M_PI*u2);
} }
//! @brief simple (and unoptimized) wrapper function
//! for the discrete fast fourier transformation library (fftw3)
/* ************************************************************************* */
// simple (and unoptimized) wrapper function for the discrete fast fourier
// transformation library (fftw3)
/* ************************************************************************* */
std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<complex_t > &data, MathFunctions::TransformCase ftCase) std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<complex_t > &data, MathFunctions::TransformCase ftCase)
{ {
double scale(1.); double scale(1.);
...@@ -95,13 +93,10 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<co ...@@ -95,13 +93,10 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<co
return outData; return outData;
} }
//! @brief simple (and unoptimized) wrapper function
//! for the discrete fast fourier transformation library (fftw3);
//! transforms real to complex
/* ************************************************************************* */
// simple (and unoptimized) wrapper function for the discrete fast fourier
// transformation library (fftw3)
//
// transforms real to complex
/* ************************************************************************* */
std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<double > &data, MathFunctions::TransformCase ftCase) std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<double > &data, MathFunctions::TransformCase ftCase)
{ {
std::vector<complex_t> cdata; std::vector<complex_t> cdata;
...@@ -110,85 +105,8 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<do ...@@ -110,85 +105,8 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<do
return MathFunctions::FastFourierTransform(cdata, ftCase); return MathFunctions::FastFourierTransform(cdata, ftCase);
} }
//complex_t MathFunctions::Bessel_J1(complex_t value) //! convolution of two real vectors of equal size
//{
// complex_t z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
// complex_t result;
// double a0;
// int k,kz;
//
// static complex_t czero(0.0, 0.0);
// static complex_t cone(1.0, 0.0);
// static double a1[] = {
// 0.1171875,
// -0.1441955566406250,
// 0.6765925884246826,
// -6.883914268109947,
// 1.215978918765359e2,
// -3.302272294480852e3,
// 1.276412726461746e5,
// -6.656367718817688e6,
// 4.502786003050393e8,
// -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};
//
// a0 = abs(value);
// z2 = value*value;
// z1 = value;
// if (a0 == 0.0) {
// result = czero;
// return result;
// }
// if (real(value) < 0.0) z1 = -value;
// if (a0 <= 12.0) {
// result = cone;
// cr = cone;
// for (k=1;k<=40;k++) {
// cr *= -0.25*z2/(k*(k+1.0));
// result += cr;
// if (abs(cr) < abs(result)*Numeric::double_epsilon) break;
// }
// result *= 0.5*z1;
// return result;
// }
// else {
// if (a0 >= 50.0) kz = 8; // can be changed to 10
// else if (a0 >= 35.0) kz = 10; // " " " 12
// else kz = 12; // " " " 14
// cp1 = cone;
// for (k=0;k<kz;k++) {
// cp1 += a1[k]*pow(z1,-2.0*k-2.0);
// }
// cq1 = 0.375/z1;
// for (k=0;k<kz;k++) {
// cq1 += b1[k]*pow(z1,-2.0*k-3.0);
// }
// result = cu*(cp1*cos(ct2)-cq1*sin(ct2));
// return result;
// }
//}
/* ************************************************************************* */
// 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() ) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment