17#include <gsl/gsl_sf_erf.h>
18#include <gsl/gsl_sf_lambert.h>
26using namespace CurveFitting;
27using namespace CurveFitting::SpecialFunctionSupport;
38 declareParameter(
"X0", -0.0,
"Location of the peak");
39 declareParameter(
"Intensity", 0.0,
"Integrated intensity");
40 declareParameter(
"Alpha", 0.04,
"Exponential rise");
41 declareParameter(
"Beta", 0.02,
"Exponential decay");
42 declareParameter(
"Sigma2", 1.0,
"Sigma squared");
43 declareParameter(
"Gamma", 0.0);
53 double minCutOff = 100.0 * std::numeric_limits<double>::min();
54 if (h0 >= 0 && h0 < minCutOff)
56 else if (h0 < 0 && h0 > -minCutOff)
80 const auto s = H / (2 * sqrt(2 * M_LN2));
86 return w0 * exp(-0.5 * M_LN2 * s / w0) + 2 * sqrt(2 * M_LN2) * s;
98 const auto a = 0.5 * M_LN2;
99 const auto b = 2 * sqrt(2 * M_LN2);
102 const auto fwhm_gauss =
103 b * w0 * (gsl_sf_lambert_W0(-(a / b) * exp(-(a / b) * (w / w0))) / a + (w / w0) / b) / 1.6364;
138 double invert_sqrt2sigma = 1.0 / sqrt(2.0 * sigma2);
139 double N = alpha * beta * 0.5 / (alpha + beta);
144 g_log.
debug() <<
"DB1143: nData = " << nData <<
" From " << xValues[0] <<
" To " << xValues[nData - 1]
145 <<
" X0 = " << x0 <<
" Intensity = " <<
intensity <<
" alpha = " << alpha <<
" beta = " << beta
146 <<
" H = " << H <<
" eta = " << eta <<
'\n';
149 double extent = 10 *
fwhm();
150 for (
size_t id = 0;
id < nData; ++id) {
151 double dT = xValues[id] - x0;
152 if (
fabs(dT) < extent) {
153 double omega =
calOmega(dT, eta, N, alpha, beta, H, sigma2, invert_sqrt2sigma);
177 double invert_sqrt2sigma)
const {
179 std::complex<double> p(alpha *
x, alpha * H * 0.5);
180 std::complex<double> q(-beta *
x, beta * H * 0.5);
182 double u = 0.5 * alpha * (alpha * sigma2 + 2 *
x);
183 double y = (alpha * sigma2 +
x) * invert_sqrt2sigma;
185 double v = 0.5 * beta * (beta * sigma2 - 2 *
x);
186 double z = (beta * sigma2 -
x) * invert_sqrt2sigma;
189 double omega1 = (1 - eta) * N * (exp(u + gsl_sf_log_erfc(
y)) + exp(v + gsl_sf_log_erfc(
z)));
196 double omega = omega1 - omega2;
207 double H_G = sqrt(8.0 * sigma2 * M_LN2);
208 double H_L = 2 * gamma;
210 double temp1 = std::pow(H_L, 5) + 0.07842 * H_G * std::pow(H_L, 4) + 4.47163 * std::pow(H_G, 2) * std::pow(H_L, 3) +
211 2.42843 * std::pow(H_G, 3) * std::pow(H_L, 2) + 2.69269 * std::pow(H_G, 4) * H_L + std::pow(H_G, 5);
213 H = std::pow(temp1, 0.2);
216 double gam_pv = H_L / H;
217 eta = 1.36603 * gam_pv - 0.47719 * std::pow(gam_pv, 2) + 0.11116 * std::pow(gam_pv, 3);
219 if (eta > 1 || eta < 0) {
220 g_log.
error() <<
"Bk2BkExpConvPV: Calculated eta = " << eta <<
" is out of range [0, 1].\n";
230 return M_LN2 * (a + b) / (a * b);
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
IPeaksWorkspace_sptr workspace
Base class that represents the domain of a function.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
virtual double getParameter(size_t i) const =0
Get i-th parameter.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
virtual double intensity() const
Returns the integral intensity of the peak.
void setParameter(const std::string &name, const double &value, bool explicitlySet=true) override
Override parent so that we may bust the cache when a parameter is set.
Represents the Jacobian in IFitFunction::functionDeriv.
Bk2BkExpConvPV : Peak profile as tback-to-back exponential convoluted with pseudo-Voigt.
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Numerical derivative.
double height() const override
Get approximate height (maximum value) of the peak (not at X0)
double centre() const override
overwrite IPeakFunction base class methods
void geneatePeak(double *out, const double *xValues, const size_t nData)
Set up the range of peak calculation for higher efficiency.
double fwhm() const override
Get approximate peak width.
void setHeight(const double h) override
Set peak height.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
double calOmega(double x, double eta, double N, double alpha, double beta, double H, double sigma2, double invert_sqrt2sigma) const
Calculate Omega(x) = ... ...
void setFwhm(const double w) override
Set new peak width approximately using same mapping in BackToBackExponential assuming fwhm_gaus = fwh...
void setCentre(const double c) override
Set peak center.
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Implement the peak calculating formula.
double expWidth() const
Calculate contribution to the width by the exponentials.
void calHandEta(double sigma2, double gamma, double &H, double &eta) const
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Local derivative.
Marks code as not implemented yet.
The Logger class is in charge of the publishing messages from the framework through various channels.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error level.
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::complex< double > MANTID_CURVEFITTING_DLL exponentialIntegral(const std::complex< double > &z)
Compute exp(z)*E1(z) where z is complex and E1(z) is the Exponential Integral.