Mantid
Loading...
Searching...
No Matches
DynamicKuboToyabe.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8
10#include "MantidAPI/Jacobian.h"
13
14#include <iomanip>
15#include <sstream>
16#include <vector>
17
19
20using namespace API;
21using namespace CurveFitting;
22using namespace Kernel;
23
24DECLARE_FUNCTION(DynamicKuboToyabe)
25
26void DynamicKuboToyabe::init() {
27 declareAttribute("BinWidth", Attribute(m_eps));
28 declareParameter("Asym", 0.2, "Amplitude at time 0");
29 declareParameter("Delta", 0.2, "Local field");
30 declareParameter("Field", 0.0, "External field");
31 declareParameter("Nu", 0.0, "Hopping rate");
32}
33
34//--------------------------------------------------------------------------------------------------------------------------------------
35// From Numerical Recipes
36
37// Midpoint method
38double midpnt(double func(const double, const double, const double), const double a, const double b, const int n,
39 const double g, const double w0) {
40 // quote & modified from numerical recipe 2nd edtion (page147)
41
42 static double s;
43
44 if (n == 1) {
45 s = (b - a) * func(0.5 * (a + b), g, w0);
46 return (s);
47 } else {
48 double x, tnm, sum, del, ddel;
49 int it, j;
50 for (it = 1, j = 1; j < n - 1; j++)
51 it *= 3;
52 tnm = it;
53 del = (b - a) / (3 * tnm);
54 ddel = del + del;
55 x = a + 0.5 * del;
56 sum = 0.0;
57 for (j = 1; j <= it; j++) {
58 sum += func(x, g, w0);
59 x += ddel;
60 sum += func(x, g, w0);
61 x += del;
62 }
63 s = (s + (b - a) * sum / tnm) / 3.0;
64 return s;
65 }
66}
67
68// Polynomial interpolation
69void polint(double xa[], const double ya[], int n, double x, double &y, double &dy) {
70 int i, m, ns = 1;
71 double dif;
72
73 dif = fabs(x - xa[1]);
74 std::vector<double> c(n + 1);
75 std::vector<double> d(n + 1);
76 for (i = 1; i <= n; i++) {
77 double dift;
78 if ((dift = fabs(x - xa[i])) < dif) {
79 ns = i;
80 dif = dift;
81 }
82 c[i] = ya[i];
83 d[i] = ya[i];
84 }
85 y = ya[ns--];
86 for (m = 1; m < n; m++) {
87 for (i = 1; i <= n - m; i++) {
88 double den, ho, hp, w;
89 ho = xa[i] - x;
90 hp = xa[i + m] - x;
91 w = c[i + 1] - d[i];
92 if ((den = ho - hp) == 0.0) { // error message!!!
93 throw std::runtime_error("Error in routin polint");
94 }
95 den = w / den;
96 d[i] = hp * den;
97 c[i] = ho * den;
98 }
99 y += (dy = (2 * (ns) < (n - m) ? c[ns + 1] : d[ns--]));
100 }
101}
102
103// Integration
104double integral(double func(const double, const double, const double), const double a, const double b, const double g,
105 const double w0) {
106
107 const int JMAX = 14;
108 const int JMAXP = JMAX + 1;
109 const int K = 5;
110
111 int j;
112 double ss, dss;
113 double h[JMAXP + 1], s[JMAXP];
114
115 h[1] = 1.0;
116 for (j = 1; j <= JMAX; j++) {
117 s[j] = midpnt(func, a, b, j, g, w0);
118 if (j >= K) {
119 polint(&h[j - K], &s[j - K], K, 0.0, ss, dss);
120 if (fabs(dss) <= fabs(ss))
121 return ss;
122 }
123 h[j + 1] = h[j] / 9.0;
124 }
125 throw std::runtime_error("Too many steps in routine integrate");
126}
127
128// End of Numerical Recipes routines
129//--------------------------------------------------------------------------------------------------------------------------------------
130
131// f1: function to integrate
132double f1(const double x, const double G, const double w0) {
133 // G = Delta in doc
134 // x = dummy time variable
135 return (exp(-G * G * x * x / 2) * sin(w0 * x));
136}
137
138// Static Zero Field Kubo Toyabe relaxation function
139// Also called Lorentzian Kubo-Toyabe
140double ZFKT(const double x, const double G) {
141 // q = Delta^2 t^2 in doc
142 const double q = G * G * x * x;
143 return (0.3333333333 + 0.6666666667 * exp(-0.5 * q) * (1 - q));
144}
145
146// Static non-zero field Kubo Toyabe relaxation function
147double HKT(const double x, const double G, const double F) {
148 // q = Delta^2 t^2 in doc
149 const double q = G * G * x * x;
150 // Muon gyromagnetic ratio * 2 * PI
151 const double gm = 2 * M_PI * PhysicalConstants::MuonGyromagneticRatio;
152 // w = omega_0 in doc
153 double w;
154 if (F > 2 * G) {
155 // Use F
156 w = gm * F;
157 } else {
158 // Use G
159 w = gm * 2 * G;
160 }
161 // r = Delta^2/omega_0^2
162 const double r = G * G / w / w;
163
164 double ig;
165 if (x > 0 && r > 0) {
166 // Compute integral
167 ig = integral(f1, 0.0, x, G, w);
168 } else {
169 // Integral is 0
170 ig = 0;
171 }
172
173 const double ktb = (1 - 2 * r * (1 - exp(-q / 2) * cos(w * x)) + 2 * r * r * w * ig);
174
175 if (F > 2 * G) {
176 // longitudinal Gaussian field
177 return ktb;
178 } else {
179 //
180 const double kz = ZFKT(x, G);
181 return kz + F / 2 / G * (ktb - kz);
182 }
183}
184
185// Dynamic Kubo-Toyabe
186double DynamicKuboToyabe::getDKT(double t, double G, double F, double v, double eps) const {
187
188 const auto tsmax = static_cast<int>(std::ceil(32.768 / eps));
189
190 static double oldG = -1., oldV = -1., oldF = -1., oldEps = -1.;
191
192 const auto maxTsmax = static_cast<int>(std::ceil(32.768 / m_minEps));
193 static std::vector<double> gStat(maxTsmax), gDyn(maxTsmax);
194
195 if ((G != oldG) || (v != oldV) || (F != oldF) || (eps != oldEps)) {
196
197 // If G or v or F or eps have changed with respect to the
198 // previous call, we need to re-do the computations
199
200 if (G != oldG || (F != oldF)) {
201
202 // But we only need to
203 // re-compute gStat if G or F have changed
204
205 // Generate static Kubo-Toyabe
206 if (F == 0) {
207 for (int k = 0; k < tsmax; k++) {
208 gStat[k] = ZFKT(k * eps, G);
209 }
210 } else {
211 for (int k = 0; k < tsmax; k++) {
212 gStat[k] = HKT(k * eps, G, F);
213 }
214 }
215 // Store new G value
216 oldG = G;
217 // Store new F value
218 oldF = F;
219 }
220
221 // Store new v value
222 oldV = v;
223 // Store new eps value
224 oldEps = eps;
225
226 double hop = v * eps;
227
228 // Generate dynamic Kubo Toyabe
229 for (int k = 0; k < tsmax; k++) {
230 double y = gStat[k];
231 // do integration
232 for (int j = k - 1; j > 0; j--) {
233 y = y * (1 - hop) + hop * gDyn[k - j] * gStat[j];
234 }
235 gDyn[k] = y;
236 }
237 }
238
239 // Interpolate table
240 // If beyond end, extrapolate
241 auto x = int(fabs(t) / eps);
242 if (x > tsmax - 2)
243 x = tsmax - 2;
244 double xe = (fabs(t) / eps) - x;
245 return gDyn[x] * (1 - xe) + xe * gDyn[x + 1];
246}
247
248// Dynamic Kubo Toyabe function
249void DynamicKuboToyabe::function1D(double *out, const double *xValues, const size_t nData) const {
250 const double &A = getParameter("Asym");
251 const double &G = fabs(getParameter("Delta"));
252 const double &F = fabs(getParameter("Field"));
253 const double &v = fabs(getParameter("Nu"));
254
255 // Zero hopping rate
256 if (v == 0.0) {
257
258 // Zero external field
259 if (F == 0.0) {
260 for (size_t i = 0; i < nData; i++) {
261 out[i] = A * ZFKT(xValues[i], G);
262 }
263 }
264 // Non-zero external field
265 else {
266 for (size_t i = 0; i < nData; i++) {
267 out[i] = A * HKT(xValues[i], G, F);
268 }
269 }
270 }
271
272 // Non-zero hopping rate
273 else {
274
275 for (size_t i = 0; i < nData; i++) {
276 out[i] = A * getDKT(xValues[i], G, F, v, m_eps);
277 }
278 }
279}
280
281//----------------------------------------------------------------------------------------------
284DynamicKuboToyabe::DynamicKuboToyabe() : m_eps(0.05), m_minEps(0.001), m_maxEps(0.1) {}
285
286//----------------------------------------------------------------------------------------------
290 calNumericalDeriv(domain, jacobian);
291}
292
293//----------------------------------------------------------------------------------------------
296void DynamicKuboToyabe::functionDeriv1D(API::Jacobian * /*jacobian*/, const double * /*xValues*/,
297 const size_t /*nData*/) {
298 throw Mantid::Kernel::Exception::NotImplementedError("functionDeriv1D is not implemented for DynamicKuboToyabe.");
299}
300
301//----------------------------------------------------------------------------------------------
307
308//----------------------------------------------------------------------------------------------
313void DynamicKuboToyabe::setAttribute(const std::string &attName, const API::IFunction::Attribute &att) {
314 if (attName == "BinWidth") {
315
316 double newVal = att.asDouble();
317
318 if (newVal < 0) {
320 throw std::invalid_argument("DKT: Attribute BinWidth cannot be negative.");
321
322 } else if (newVal < m_minEps) {
324 std::stringstream ss;
325 ss << "DKT: Attribute BinWidth too small (BinWidth < " << std::setprecision(3) << m_minEps << ")";
326 throw std::invalid_argument(ss.str());
327
328 } else if (newVal > m_maxEps) {
330 std::stringstream ss;
331 ss << "DKT: Attribute BinWidth too large (BinWidth > " << std::setprecision(3) << m_maxEps << ")";
332 throw std::invalid_argument(ss.str());
333 }
334
335 if (!nParams()) {
336 init();
337 }
338 m_eps = newVal;
340
341 } else {
342 throw std::invalid_argument("DynamicKuboToyabe: Unknown attribute " + attName);
343 }
344}
345
346} // namespace Mantid::CurveFitting::Functions
double value
The value of the point.
Definition: FitMW.cpp:51
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
Mantid::API::IFunction::Attribute Attribute
Definition: IsoRotDiff.cpp:25
#define fabs(x)
Definition: Matrix.cpp:22
Base class that represents the domain of a function.
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
Definition: IFunction.cpp:1409
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Definition: IFunction.cpp:1031
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
void clearAllParameters()
Nonvirtual member which removes all declared parameters.
void setParameter(size_t, const double &value, bool explicitlySet=true) override
Set i-th parameter.
size_t nParams() const override
Total number of parameters.
Definition: ParamFunction.h:53
double getParameter(size_t i) const override
Get i-th parameter.
Provide Dynamic Kubo Toyabe function interface to IFunction1D for muon scientists.
void init() override
Function initialization. Declare function parameters in this method.
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Function to calculate derivative numerically.
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
double getDKT(double t, double G, double F, double v, double eps) const
void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) override
Function to calculate derivative analytically.
void setActiveParameter(size_t i, double value) override
Set new value of the i-th parameter.
Marks code as not implemented yet.
Definition: Exception.h:138
double f1(const double x, const double G, const double w0)
double integral(double func(const double, const double, const double), const double a, const double b, const double g, const double w0)
double ZFKT(const double x, const double G)
double HKT(const double x, const double G, const double F)
void polint(double xa[], const double ya[], int n, double x, double &y, double &dy)
double midpnt(double func(const double, const double, const double), const double a, const double b, const int n, const double g, const double w0)
static constexpr double MuonGyromagneticRatio
Muon gyromagnetic ratio in MHz/G Taken from CalMuonDetectorPhases and DynamicKuboToyabe on 02/02/2016...