21using namespace CurveFitting;
22using namespace Kernel;
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");
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) {
45 s = (b - a) * func(0.5 * (a + b), g, w0);
48 double x, tnm, sum, del, ddel;
50 for (it = 1, j = 1; j <
n - 1; j++)
53 del = (b - a) / (3 * tnm);
57 for (j = 1; j <= it; j++) {
58 sum += func(
x, g, w0);
60 sum += func(
x, g, w0);
63 s = (s + (b - a) * sum / tnm) / 3.0;
69void polint(
double xa[],
const double ya[],
int n,
double x,
double &
y,
double &dy) {
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++) {
78 if ((dift =
fabs(
x - xa[i])) < dif) {
86 for (
m = 1;
m <
n;
m++) {
87 for (i = 1; i <=
n -
m; i++) {
88 double den, ho, hp, w;
92 if ((den = ho - hp) == 0.0) {
93 throw std::runtime_error(
"Error in routin polint");
99 y += (dy = (2 * (ns) < (
n -
m) ? c[ns + 1] :
d[ns--]));
104double integral(
double func(
const double,
const double,
const double),
const double a,
const double b,
const double g,
108 const int JMAXP = JMAX + 1;
113 double h[JMAXP + 1], s[JMAXP];
116 for (j = 1; j <= JMAX; j++) {
117 s[j] =
midpnt(func, a, b, j, g, w0);
119 polint(&h[j - K], &s[j - K], K, 0.0, ss, dss);
123 h[j + 1] = h[j] / 9.0;
125 throw std::runtime_error(
"Too many steps in routine integrate");
132double f1(
const double x,
const double G,
const double w0) {
135 return (exp(-G * G *
x *
x / 2) * sin(w0 *
x));
140double ZFKT(
const double x,
const double G) {
142 const double q = G * G *
x *
x;
143 return (0.3333333333 + 0.6666666667 * exp(-0.5 * q) * (1 - q));
147double HKT(
const double x,
const double G,
const double F) {
149 const double q = G * G *
x *
x;
162 const double r = G * G / w / w;
165 if (
x > 0 && r > 0) {
173 const double ktb = (1 - 2 * r * (1 - exp(-q / 2) * cos(w *
x)) + 2 * r * r * w * ig);
180 const double kz =
ZFKT(
x, G);
181 return kz + F / 2 / G * (ktb - kz);
188 const auto tsmax =
static_cast<int>(std::ceil(32.768 / eps));
190 static double oldG = -1., oldV = -1., oldF = -1., oldEps = -1.;
192 const auto maxTsmax =
static_cast<int>(std::ceil(32.768 /
m_minEps));
193 static std::vector<double> gStat(maxTsmax), gDyn(maxTsmax);
195 if ((G != oldG) || (v != oldV) || (F != oldF) || (eps != oldEps)) {
200 if (G != oldG || (F != oldF)) {
207 for (
int k = 0; k < tsmax; k++) {
208 gStat[k] =
ZFKT(k * eps, G);
211 for (
int k = 0; k < tsmax; k++) {
212 gStat[k] =
HKT(k * eps, G, F);
226 double hop = v * eps;
229 for (
int k = 0; k < tsmax; k++) {
232 for (
int j = k - 1; j > 0; j--) {
233 y =
y * (1 - hop) + hop * gDyn[k - j] * gStat[j];
241 auto x = int(
fabs(t) / eps);
244 double xe = (
fabs(t) / eps) -
x;
245 return gDyn[
x] * (1 - xe) + xe * gDyn[
x + 1];
260 for (
size_t i = 0; i < nData; i++) {
261 out[i] = A *
ZFKT(xValues[i], G);
266 for (
size_t i = 0; i < nData; i++) {
267 out[i] = A *
HKT(xValues[i], G, F);
275 for (
size_t i = 0; i < nData; i++) {
314 if (attName ==
"BinWidth") {
320 throw std::invalid_argument(
"DKT: Attribute BinWidth cannot be negative.");
324 std::stringstream ss;
325 ss <<
"DKT: Attribute BinWidth too small (BinWidth < " << std::setprecision(3) <<
m_minEps <<
")";
326 throw std::invalid_argument(ss.str());
330 std::stringstream ss;
331 ss <<
"DKT: Attribute BinWidth too large (BinWidth > " << std::setprecision(3) <<
m_maxEps <<
")";
332 throw std::invalid_argument(ss.str());
342 throw std::invalid_argument(
"DynamicKuboToyabe: Unknown attribute " + attName);
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.
Mantid::API::IFunction::Attribute Attribute
Base class that represents the domain of a function.
Attribute is a non-fitting parameter.
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Represents the Jacobian in IFitFunction::functionDeriv.
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.
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.
DynamicKuboToyabe()
Constructor.
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.
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...