16#include <gsl/gsl_sf_erf.h>
25using namespace CurveFitting;
28Kernel::Logger
g_log(
"NeutronBk2BkExpConvPV");
54 declareParameter(
"Dtt1", 1.0,
"coefficient 1 for d-spacing calculation for epithermal neutron part");
55 declareParameter(
"Dtt2", 1.0,
"coefficient 2 for d-spacing calculation for epithermal neutron part");
59 declareParameter(
"Alph0", 1.6,
"exponential constant for rising part of epithermal neutron pulse");
60 declareParameter(
"Alph1", 1.5,
"exponential constant for rising part of expithermal neutron pulse");
61 declareParameter(
"Beta0", 1.6,
"exponential constant of decaying part of epithermal neutron pulse");
62 declareParameter(
"Beta1", 1.5,
"exponential constant of decaying part of epithermal neutron pulse");
66 "variance parameter 1 of the Gaussian "
67 "component of the psuedovoigt function");
69 "variance parameter 2 of the Gaussian "
70 "component of the psuedovoigt function");
72 "variance parameter 3 of the Gaussian "
73 "component of the psuedovoigt function");
76 "FWHM parameter 1 of the Lorentzian component "
77 "of the psuedovoigt function");
79 "FWHM parameter 2 of the Lorentzian component "
80 "of the psuedovoigt function");
82 "FWHM parameter 3 of the Lorentzian component "
83 "of the psuedovoigt function");
86 declareParameter(
"LatticeConstant", 10.0,
"lattice constant for the sample");
119 if (paramname ==
"Alpha")
121 else if (paramname ==
"Beta")
123 else if (paramname ==
"Sigma2")
125 else if (paramname ==
"Gamma")
127 else if (paramname ==
"d_h")
129 else if (paramname ==
"Eta")
131 else if (paramname ==
"TOF_h")
133 else if (paramname ==
"FWHM")
137 errss <<
"Parameter " << paramname <<
" does not exist in peak function " << this->
name()
138 <<
"'s calculated parameters. "
139 <<
"Candidates are Alpha, Beta, Sigma2, Gamma, d_h and FWHM. ";
141 throw runtime_error(errss.str());
177 mH =
static_cast<int>(h);
178 mK =
static_cast<int>(k);
179 mL =
static_cast<int>(l);
184 errmsg <<
"H = K = L = 0 is not allowed";
186 throw std::invalid_argument(errmsg.str());
189 g_log.
debug() <<
"Set (HKL) from input parameter = " <<
mH <<
", " <<
mK <<
", " <<
mL <<
"\n";
194 double dh, tof_h, eta, alpha, beta, H, sigma2, gamma, N;
199 m_unitCell.
set(latticeconstant, latticeconstant, latticeconstant, 90.0, 90.0, 90.0);
214 alpha = alph0 + alph1 / dh;
215 beta = beta0 + beta1 / pow(dh, 4.);
216 tof_h = zero + dtt1 * dh + dtt2 * dh * dh;
218 sigma2 = sig0 * sig0 + sig1 * sig1 * std::pow(dh, 2) + sig2 * sig2 * std::pow(dh, 4);
219 gamma = gam0 + gam1 * dh + gam2 * std::pow(dh, 2);
221 g_log.
debug() <<
"[F001] TOF_h = " << tof_h <<
", Alpha = " << alpha <<
", Beta = " << beta <<
", Gamma = " << gamma
222 <<
"(Gam-0 = " << gam0 <<
", Gam-1 = " << gam1 <<
", Gam-2 = " << gam2 <<
")."
228 N = alpha * beta * 0.5 / (alpha + beta);
241 if (alpha != alpha || beta != beta || sigma2 != sigma2 || gamma != gamma || H != H || H <= 0.) {
248 if (explicitoutput) {
250 msgss <<
" dh = " << dh <<
"; TOF = " << tof_h <<
", "
251 <<
"alpha = " << alpha <<
", beta = " << beta;
283 if (
name ==
"LatticeConstant") {
310 g_log.
debug(
"Function() has no new parameters to calculate. ");
314 const double INVERT_SQRT2SIGMA = 1.0 / sqrt(2.0 *
m_Sigma2);
319 const double LEFT_VALUE =
m_centre - RANGE;
320 auto iter = std::lower_bound(xValues.cbegin(), xValues.cend(), LEFT_VALUE);
322 const double RIGHT_VALUE =
m_centre + RANGE;
323 auto iter_end = std::lower_bound(iter, xValues.cend(), RIGHT_VALUE);
326 std::size_t pos(std::distance(xValues.begin(), iter));
327 for (; iter != iter_end; ++iter) {
341 g_log.
debug(
"Function() has no new parameters to calculate. ");
345 const double INVERT_SQRT2SIGMA = 1.0 / sqrt(2.0 *
m_Sigma2);
349 g_log.
debug() <<
"[F002] Peak centre = " <<
m_centre <<
"; Calcualtion Range = " << RANGE <<
".\n";
351 for (
size_t i = 0; i < nData; ++i) {
356 g_log.
debug() <<
"TOF = " << xValues[i] <<
" = " << out[i] <<
"\n";
360 g_log.
debug() <<
"TOF = " << xValues[i] <<
" out of calculation range. "
371 double H_G = sqrt(8.0 * sigma2 * M_LN2);
374 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) +
375 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);
377 H = std::pow(temp1, 0.2);
380 double gam_pv = H_L / H;
381 eta = 1.36603 * gam_pv - 0.47719 * std::pow(gam_pv, 2) + 0.11116 * std::pow(gam_pv, 3);
383 if (eta > 1 || eta < 0) {
384 g_log.
warning() <<
"Calculated eta = " << eta <<
" is out of range [0, 1].\n";
386 g_log.
debug() <<
"[DBx121] Eta = " << eta <<
"; Gamma = " << gamma <<
".\n";
395 const double beta,
const double H,
const double sigma2,
396 const double invert_sqrt2sigma,
const bool explicitoutput)
const {
398 const double u = 0.5 * alpha * (alpha * sigma2 + 2. *
x);
399 const double y = (alpha * sigma2 +
x) * invert_sqrt2sigma;
401 const double v = 0.5 * beta * (beta * sigma2 - 2. *
x);
402 const double z = (beta * sigma2 -
x) * invert_sqrt2sigma;
405 const double erfcy = gsl_sf_erfc(
y);
407 if (
fabs(erfcy) > DBL_MIN)
408 part1 = exp(u) * erfcy;
410 const double erfcz = gsl_sf_erfc(
z);
412 if (
fabs(erfcz) > DBL_MIN)
413 part2 = exp(v) * erfcz;
415 const double omega1 = (1. - eta) * N * (part1 + part2);
420 g_log.
debug() <<
"Eta = " << eta <<
"; X = " <<
x <<
" N = " << N <<
".\n";
423 const double SQRT_H_5 = sqrt(H) * .5;
424 std::complex<double> p(alpha *
x, alpha * SQRT_H_5);
425 std::complex<double> q(-beta *
x, beta * SQRT_H_5);
428 omega2 = -1.0 * N * eta * (omega2a + omega2b) * M_2_PI;
430 g_log.
debug() <<
"Exp(p) = " << exp(p) <<
", Exp(q) = " << exp(q) <<
".\n";
432 if (omega2 != omega2) {
433 g_log.
debug() <<
"Omega2 is not physical. Omega2a = " << omega2a <<
", Omega2b = " << omega2b
434 <<
", p = " << p.real() <<
", " << p.imag() <<
".\n";
436 g_log.
debug() <<
"X = " <<
x <<
" is OK. Omega 2 = " << omega2 <<
", Omega2A = " << omega2a
437 <<
", Omega2B = " << omega2b <<
"\n";
440 const double omega = omega1 + omega2;
442 if (explicitoutput || omega != omega) {
443 if (omega <= NEG_DBL_MAX || omega >= DBL_MAX || omega != omega) {
447 errss <<
"Omega = " << omega <<
" is infinity! omega1 = " << omega1 <<
", omega2 = " << omega2 <<
"\n";
448 errss <<
" u = " << u <<
", v = " << v <<
", erfc(y) = " << gsl_sf_erfc(
y) <<
", erfc(z) = " << gsl_sf_erfc(
z)
450 errss <<
" alpha = " << alpha <<
", beta = " << beta <<
" sigma2 = " << sigma2 <<
", N = " << N <<
"\n";
454 g_log.
debug() <<
"[DB] Final Value of Omega = " << omega <<
".\n";
double value
The value of the point.
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
static Kernel::Logger g_log
Logger instance.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
bool m_cellParamValueChanged
An indicator to re-calculate peak d-space position.
bool m_parameterValid
Flag to indicate whether peaks' parameters value can generate a valid peak.
double m_unitCellSize
Unit cell size.
Geometry::UnitCell m_unitCell
Unit cell.
double m_centre
Centre of the peak.
bool m_hasNewParameterValue
Flag if any parameter value changed.
double m_dcentre
Centre of the peak in d-space.
double m_fwhm
Peak's FWHM.
void setParameter(size_t, const double &value, bool explicitlySet=true) override
Set i-th parameter.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
double getParameter(size_t i) const override
Get i-th parameter.
NeutronBk2BkExpConvPVoigt : Back-to-back exponential function convoluted with pseudo-voigt for epithe...
void calculateParameters(bool explicitoutput) const override
Calculate peak parameters (alpha, beta, sigma2..)
void setParameter(size_t i, const double &value, bool explicitlySet=true) override
Override setting a new value to the i-th parameter.
void function(std::vector< double > &out, const std::vector< double > &xValues) const override
Function (local) of the vector version.
static int s_peakRadius
Default value for the peak radius.
double getPeakParameter(const std::string &) override
Get peak parameters.
void init() override
Fuction local.
std::string name() const override
Overwrite IFunction base class method: name.
void calHandEta(double sigma2, double gamma, double &H, double &eta) const
Calcualte H and Eta.
double m_Alpha
Set 2 functions to be hidden from client.
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
double calOmega(const double x, const double eta, const double N, const double alpha, const double beta, const double H, const double sigma2, const double invert_sqrt2sigma, const bool explicitoutput=false) const
Calculate peak profile I(TOF) = Omega(TOF)
void set(double _a, double _b, double _c, double _alpha, double _beta, double _gamma, const int angleunit=angDegrees)
Set lattice parameters.
double d(double h, double k, double l) const
Return d-spacing ( ) for a given h,k,l coordinate.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
Kernel::Logger g_log("ExperimentInfo")
static logger object
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.