/
Voigt.cpp
221 lines (196 loc) · 7.56 KB
/
Voigt.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------------------------------
#include "MantidCurveFitting/Functions/Voigt.h"
#include "MantidAPI/FunctionFactory.h"
#include <cmath>
#include <limits>
namespace Mantid::CurveFitting::Functions {
using namespace CurveFitting;
DECLARE_FUNCTION(Voigt)
namespace {
/// @cond
// Coefficients for approximation using 4 lorentzians
const size_t NLORENTZIANS = 4;
const double COEFFA[NLORENTZIANS] = {-1.2150, -1.3509, -1.2150, -1.3509};
const double COEFFB[NLORENTZIANS] = {1.2359, 0.3786, -1.2359, -0.3786};
const double COEFFC[NLORENTZIANS] = {-0.3085, 0.5906, -0.3085, 0.5906};
const double COEFFD[NLORENTZIANS] = {0.0210, -1.1858, -0.0210, 1.1858};
const char *LORENTZ_AMP = "LorentzAmp";
const char *LORENTZ_POS = "LorentzPos";
const char *LORENTZ_FWHM = "LorentzFWHM";
const char *GAUSSIAN_FWHM = "GaussianFWHM";
const double SQRTLN2 = std::sqrt(M_LN2);
const double SQRTPI = std::sqrt(M_PI);
///@endcond
} // namespace
/**
* Declare the active parameters for the function
*/
void Voigt::declareParameters() {
declareParameter(LORENTZ_AMP, 0.0, "Value of the Lorentzian amplitude");
declareParameter(LORENTZ_POS, 0.0, "Position of the Lorentzian peak");
declareParameter(LORENTZ_FWHM, 0.0, "Value of the full-width half-maximum for the Lorentzian");
declareParameter(GAUSSIAN_FWHM, 0.0, "Value of the full-width half-maximum for the Gaussian");
}
/**
* Calculate Voigt function for each x value
* @param out :: The values of the function at each x point
* @param xValues :: The X values
* @param nData :: The number of X values to evaluate
*/
void Voigt::functionLocal(double *out, const double *xValues, const size_t nData) const {
calculateFunctionAndDerivative(xValues, nData, out, nullptr);
}
/**
* Derivatives of function with respect to active parameters
* @param out :: The Jacobian matrix containing the partial derivatives for each
* x value
* @param xValues :: The X values
* @param nData :: The number of X values to evaluate
*/
void Voigt::functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) {
calculateFunctionAndDerivative(xValues, nData, nullptr, out);
}
/**
* Calculates both function & derivative together
* @param xValues :: The X values
* @param nData :: The number of X values to evaluate
* @param functionValues :: Calculated y values
* @param derivatives :: The Jacobian matrix containing the partial derivatives
* for each x value (allowed null)
*/
void Voigt::calculateFunctionAndDerivative(const double *xValues, const size_t nData, double *functionValues,
API::Jacobian *derivatives) const {
const double a_L = getParameter(LORENTZ_AMP);
const double lorentzPos = getParameter(LORENTZ_POS);
const double gamma_L = getParameter(LORENTZ_FWHM);
const double gamma_G = getParameter(GAUSSIAN_FWHM);
const double rtln2oGammaG = SQRTLN2 / gamma_G;
const double prefactor = (a_L * SQRTPI * gamma_L * SQRTLN2 / gamma_G);
for (size_t i = 0; i < nData; ++i) {
const double xoffset = xValues[i] - lorentzPos;
const double X = xoffset * 2.0 * rtln2oGammaG;
const double Y = gamma_L * rtln2oGammaG;
double fx(0.0), dFdx(0.0), dFdy(0.0);
for (size_t j = 0; j < NLORENTZIANS; ++j) {
const double ymA(Y - COEFFA[j]);
const double xmB(X - COEFFB[j]);
const double alpha = COEFFC[j] * ymA + COEFFD[j] * xmB;
const double beta = ymA * ymA + xmB * xmB;
const double ratioab = alpha / beta;
fx += ratioab;
dFdx += (COEFFD[j] / beta) - 2.0 * (X - COEFFB[j]) * ratioab / beta;
dFdy += (COEFFC[j] / beta) - 2.0 * (Y - COEFFA[j]) * ratioab / beta;
}
if (functionValues) {
functionValues[i] = prefactor * fx;
}
if (derivatives) {
derivatives->set(i, 0, prefactor * fx / a_L);
derivatives->set(i, 1, -prefactor * dFdx * 2.0 * rtln2oGammaG);
derivatives->set(i, 2, prefactor * (fx / gamma_L + dFdy * rtln2oGammaG));
derivatives->set(i, 3, -prefactor * (fx + (rtln2oGammaG) * (2.0 * xoffset * dFdx + gamma_L * dFdy)) / gamma_G);
}
}
}
/**
* Returns the value of the "LorentzPos" parameter.
* @return value of centre of peak
*/
double Voigt::centre() const { return getParameter(LORENTZ_POS); }
/**
* Return the value of the "LorentzAmp" parameter
* @return value of height of peak
*/
double Voigt::height() const {
if (getParameter(LORENTZ_AMP) == 0.0 || getParameter(LORENTZ_FWHM) == 0.0 || getParameter(GAUSSIAN_FWHM) == 0.0) {
return 0.0;
}
double pos = getParameter(LORENTZ_POS);
double h;
functionLocal(&h, &pos, 1);
return h;
}
/**
* Gives the FWHM of the peak. This is estimated as
* 0.5*(LorentzFWHM + GaussianFWHM)
* @return value of FWHM of peak
*/
double Voigt::fwhm() const { return (getParameter(LORENTZ_FWHM) + getParameter(GAUSSIAN_FWHM)); }
/**
* Set the centre of the peak, the LorentzPos parameter
* @param value :: The new value for the centre of the peak
*/
void Voigt::setCentre(const double value) { this->setParameter(LORENTZ_POS, value); }
/**
* Set the height of the peak. Sets LorentzAmp parameter to 1.5*value
* @param value :: The new value for the centre of the peak
*/
void Voigt::setHeight(const double value) {
auto lorentzFwhm = getParameter(LORENTZ_FWHM);
if (lorentzFwhm == 0.0) {
lorentzFwhm = std::numeric_limits<double>::epsilon();
setParameter(LORENTZ_FWHM, lorentzFwhm);
}
auto lorentzAmp = getParameter(LORENTZ_AMP);
if (lorentzAmp == 0.0) {
lorentzAmp = std::numeric_limits<double>::epsilon();
setParameter(LORENTZ_AMP, lorentzAmp);
}
auto gaussFwhm = getParameter(GAUSSIAN_FWHM);
if (gaussFwhm == 0.0) {
setParameter(GAUSSIAN_FWHM, std::numeric_limits<double>::epsilon());
}
auto h = height();
this->setParameter(LORENTZ_AMP, lorentzAmp * value / h);
}
/**
* Set the FWHM of the peak
* @param value :: The new value for the FWHM of the peak
*/
void Voigt::setFwhm(const double value) {
auto lorentzFwhm = getParameter(LORENTZ_FWHM);
if (lorentzFwhm == 0.0) {
lorentzFwhm = std::numeric_limits<double>::epsilon();
}
auto gaussFwhm = getParameter(GAUSSIAN_FWHM);
if (gaussFwhm == 0.0) {
gaussFwhm = std::numeric_limits<double>::epsilon();
}
auto ratio = lorentzFwhm / (lorentzFwhm + gaussFwhm);
this->setParameter(LORENTZ_FWHM, ratio * value);
this->setParameter(GAUSSIAN_FWHM, (1.0 - ratio) * value);
}
/**
* Returns the integral intensity of the peak
*/
double Voigt::intensity() const {
if (getParameter(GAUSSIAN_FWHM) == 0.0) {
return 0.0;
}
return M_PI * getParameter(LORENTZ_AMP) * getParameter(LORENTZ_FWHM) / 2.0;
}
/**
* Sets the integral intensity of the peak
* @param value :: The new value for the intensity.
*/
void Voigt::setIntensity(const double value) {
auto lorentzFWHM = getParameter(LORENTZ_FWHM);
if (lorentzFWHM == 0.0) {
lorentzFWHM = std::numeric_limits<double>::epsilon();
setParameter(LORENTZ_FWHM, lorentzFWHM);
}
auto gaussFwhm = getParameter(GAUSSIAN_FWHM);
if (gaussFwhm == 0.0) {
setParameter(GAUSSIAN_FWHM, std::numeric_limits<double>::epsilon());
}
setParameter(LORENTZ_AMP, 2.0 * value / (M_PI * lorentzFWHM));
}
} // namespace Mantid::CurveFitting::Functions