-
Notifications
You must be signed in to change notification settings - Fork 122
/
Keren.cpp
129 lines (115 loc) · 4.63 KB
/
Keren.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
// 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 +
#include "MantidCurveFitting/Functions/Keren.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidKernel/PhysicalConstants.h"
#include <cmath>
namespace {
constexpr double TWOPI = 2.0 * M_PI;
}
namespace Mantid::CurveFitting::Functions {
DECLARE_FUNCTION(Keren)
/**
* Initialise parameters for function
*/
void Keren::init() {
declareParameter("A", 1.0, "Polarization at time zero");
declareParameter("Delta", 0.2, "Distribution width of local fields (MHz)");
declareParameter("Field", 50.0, "Longitudinal field (Gauss)");
declareParameter("Fluct", 0.2, "Hopping rate (inverse correlation time, MHz)");
}
/**
* Set the ith active parameter's value
* We use the Larmor frequency rather than the field
* omega_L = 2 * pi * gamma_mu * B
* @param i :: [input] Index of parameter
* @param value :: [input, output] Value of ith parameter
* @throws std::runtime_error if param is inactive (fixed)
*/
void Keren::setActiveParameter(size_t i, double value) {
if (!isActive(i)) {
throw std::runtime_error("Attempt to use an inactive parameter");
}
if (parameterName(i) == "Field") {
// Value passed in is omega_L, return B
setParameter(i, value / (PhysicalConstants::MuonGyromagneticRatio * TWOPI), false);
} else
setParameter(i, value, false);
}
/**
* Get the ith active parameter's value
* We use the Larmor frequency rather than the field
* omega_L = 2 * pi * gamma_mu * B
* @param i :: [input] Index of parameter
* @returns :: Value of ith parameter
* @throws std::runtime_error if param is inactive (fixed)
*/
double Keren::activeParameter(size_t i) const {
if (!isActive(i)) {
throw std::runtime_error("Attempt to use an inactive parameter");
}
if (parameterName(i) == "Field") {
// Fit to omega_L = 2 * pi * gamma_mu * B
return getParameter(i) * PhysicalConstants::MuonGyromagneticRatio * TWOPI;
} else
return getParameter(i);
}
/**
* Calculate the value of the function for each given X value
* @param out :: [output] Array to be filled with function values
* @param xValues :: [input] Array of X values to calculate function at
* @param nData :: [input] Number of X values
*/
void Keren::function1D(double *out, const double *xValues, const size_t nData) const {
// Get parameters
const double initial = getParameter("A");
const double delta = getParameter("Delta");
const double fluct = getParameter("Fluct");
const double larmor = getParameter("Field") * PhysicalConstants::MuonGyromagneticRatio * TWOPI;
for (size_t i = 0; i < nData; i++) {
out[i] = initial * polarization(delta, larmor, fluct, xValues[i]);
}
}
/**
* Calculate the muon polarization P_z(t) at time t as a fraction of the initial
* polarization (P_0(t) = 1)
* Equation (3) in the paper
* @param delta :: [input] Delta, distribution width of local fields in MHz
* @param larmor :: [input] omega_L, Larmor frequency in MHz
* @param fluct :: [input] nu, hopping rate in MHz
* @param time :: [input] t, time in microseconds
* @returns :: Polarization P_z(t) (dimensionless)
*/
double Keren::polarization(const double delta, const double larmor, const double fluct, const double time) const {
return exp(-1.0 * relaxation(delta, larmor, fluct, time));
}
/**
* Calculate the relaxation Gamma(t)t at time t
* Equation (4) in the paper
* @param delta :: [input] Delta, distribution width of local fields in MHz
* @param larmor :: [input] omega_L, Larmor frequency in MHz
* @param fluct :: [input] nu, hopping rate in MHz
* @param time :: [input] t, time in microseconds
* @returns :: Relaxation Gamma(t)*t (dimensionless)
*/
double Keren::relaxation(const double delta, const double larmor, const double fluct, const double time) const {
// Useful shortcuts
const double deltaSq = delta * delta;
const double omegaSq = larmor * larmor;
const double nuSq = fluct * fluct;
const double omegaT = larmor * time;
const double nuT = fluct * time;
const double expon = exp(-1.0 * nuT);
// 2*Delta^2 / (omega_L^2 + nu^2)^2
const double prefactor = (2.0 * deltaSq) / ((omegaSq + nuSq) * (omegaSq + nuSq));
const double term1 = (omegaSq + nuSq) * nuT;
const double term2 = omegaSq - nuSq;
const double term3 = 1.0 - expon * cos(omegaT);
const double term4 = 2.0 * fluct * larmor * expon * sin(omegaT);
return prefactor * (term1 + term2 * term3 - term4);
}
} // namespace Mantid::CurveFitting::Functions