14#include <gsl/gsl_multifit_nlin.h>
15#include <gsl/gsl_sf_erf.h>
16#include <gsl/gsl_sf_lambert.h>
20using namespace CurveFitting;
22using namespace Kernel;
30 declareParameter(
"I", 0.0,
"integrated intensity of the peak");
31 declareParameter(
"A", 1.0,
32 "exponential constant of rising part of neutron pulse");
33 declareParameter(
"B", 0.05,
"exponential constant of decaying part of neutron pulse");
34 declareParameter(
"X0", 0.0,
"peak position");
35 declareParameter(
"S", 1.0,
36 "standard deviation of gaussian part of peakshape function");
46 const auto width =
fwhm();
47 const auto h = I * M_LN2 / (width - M_LN2 * s);
66 if (!std::isfinite(area)) {
67 area = std::numeric_limits<double>::max() / 2;
82 return w0 * exp(-0.5 * M_LN2 * s / w0) + 2 * sqrt(2 * M_LN2) * s;
94 const auto a = 0.5 * M_LN2;
95 const auto b = 2 * sqrt(2 * M_LN2);
97 setParameter(
"S", w0 * (gsl_sf_lambert_W0(-(a / b) * exp(-(a / b) * (w / w0))) / a + (w / w0) / b));
128 double normFactor = a * b / (a + b) / 2;
130 if (normFactor == 0.0)
132 for (
size_t i = 0; i < nData; i++) {
133 double diff = xValues[i] - x0;
134 if (
fabs(diff) < extent) {
136 double arg1 = a / 2 * (a * s2 + 2 * diff);
137 val += exp(arg1 + gsl_sf_log_erfc((a * s2 + diff) / sqrt(2 * s2)));
138 double arg2 = b / 2 * (b * s2 - 2 * diff);
139 val += exp(arg2 + gsl_sf_log_erfc((b * s2 - diff) / sqrt(2 * s2)));
140 out[i] = I * val * normFactor;
163 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.
1D domain - a wrapper around an array of doubles.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
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.
double getParameter(size_t i) const override
Get i-th parameter.
Provide BackToBackExponential peak shape function interface to IPeakFunction.
void setFwhm(const double w) override
Set new peak width approximately.
void function1D(double *out, const double *xValues, const size_t nData) const override
General implementation of the method for all peaks.
double height() const override
Get approximate height (maximum value) of the peak (not at X0)
void setHeight(const double h) override
Set new height of the peak.
double fwhm() const override
Get approximate peak width.
void functionDeriv1D(API::Jacobian *jacobian, const double *xValues, const size_t nData) override
Evaluate function derivatives numerically.
double expWidth() const
Calculate contribution to the width by the exponentials.