Mantid
Loading...
Searching...
No Matches
ReflectivityMulf.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 +
9#include <cmath>
10
11using namespace Mantid::Kernel;
12
13using namespace Mantid::API;
14
15using namespace std;
16
18
19using namespace CurveFitting;
20
21DECLARE_FUNCTION(ReflectivityMulf)
22
23//----------------------------------------------------------------------------------------------
26ReflectivityMulf::ReflectivityMulf() : m_nlayer(0), m_nlayer_old(0) {}
27
28// Initialize Basic Parameters
30 m_nlayer_old = 0;
31 declareAttribute("nlayer", Attribute(0));
32 declareParameter("Theta", 2.3);
33 declareParameter("ScaleFactor", 1.0);
34 declareParameter("AirSLD", 0.0);
35 declareParameter("BulkSLD", 6.35e-6);
36 declareParameter("Roughness", 2.5);
37 declareParameter("BackGround", 1.0e-6);
38 declareParameter("Resolution", 5.0);
39}
40
41//----------------------------------------------------------------------------------------------
44void ReflectivityMulf::function1D(double *out, const double *xValues, const size_t nData) const {
45 // 1. Use a vector for all coefficient
46 vector<double> coeff(m_nlayer * 3 + 7, 0.0);
47
48 for (int i = 0; i < 7; ++i) {
49 coeff[i] = getParameter(i);
50 }
51 for (int i = 0; i < m_nlayer; ++i) {
52 // Layer SLD
53 coeff[7 + i * 3] = getParameter(7 + i * 3);
54 // Layer d
55 coeff[8 + i * 3] = getParameter(8 + i * 3);
56 // Layer Roughness
57 coeff[9 + i * 3] = getParameter(9 + i * 3);
58 }
59
60 // 2. Calculate
61 std::vector<double> dn(m_nlayer + 2);
62 std::vector<double> rnbn(m_nlayer + 2);
63 std::vector<double> zbb(m_nlayer + 2);
64 int nit = 21;
65 std::vector<double> cy(nData);
66
67 std::complex<double> c0(0.0, 0.0);
68 std::complex<double> ci(0.0, 1.0);
69 std::complex<double> cr(1.0, 0.0);
70 std::vector<std::complex<double>> rnfn(m_nlayer + 2);
71 std::vector<std::complex<double>> pfn(m_nlayer + 2);
72 std::vector<std::complex<double>> betan(m_nlayer + 2);
73 std::complex<double> rnf(0.0, 0.0);
74 std::complex<double> rnf1(0.0, 0.0);
75 std::complex<double> ac1(0.0, 0.0);
76 std::complex<double> ac2(0.0, 0.0);
77 std::complex<double> ac3(0.0, 0.0);
78 std::complex<double> ac4(0.0, 0.0);
79 std::complex<double> a12t(0.0, 0.0);
80 std::complex<double> a22t(0.0, 0.0);
81 std::complex<double> btm(0.0, 0.0);
82 std::complex<double> btm1(0.0, 0.0);
83 std::complex<double> cbtm(0.0, 0.0);
84 std::complex<double> cbtm1(0.0, 0.0);
85
86 std::complex<double> a111(0.0, 0.0);
87 std::complex<double> a112(0.0, 0.0);
88 std::complex<double> a121(0.0, 0.0);
89 std::complex<double> a122(0.0, 0.0);
90 std::complex<double> a211(0.0, 0.0);
91 std::complex<double> a212(0.0, 0.0);
92 std::complex<double> a221(0.0, 0.0);
93 std::complex<double> a222(0.0, 0.0);
94 std::complex<double> a311(0.0, 0.0);
95 std::complex<double> a312(0.0, 0.0);
96 std::complex<double> a321(0.0, 0.0);
97 std::complex<double> a322(0.0, 0.0);
98
99 double theta0 = coeff[0] * M_PI / 180.0;
100 double scalefac = coeff[1];
101 rnbn[0] = coeff[2];
102 rnbn[m_nlayer + 1] = coeff[3];
103 zbb[0] = coeff[4] * coeff[4];
104 zbb[m_nlayer + 1] = 0.0;
105 double bgd = coeff[5];
106 double pthet = coeff[6];
107
108 dn[0] = 0.0;
109 dn[m_nlayer + 1] = 0.0;
110
111 if (m_nlayer > 0) {
112 for (int i = 0; i < m_nlayer; ++i) {
113 rnbn[i + 1] = coeff[7 + i * 3];
114 dn[i + 1] = coeff[8 + i * 3];
115 zbb[i + 1] = coeff[9 + i * 3];
116 }
117 }
118
119 double tmax, tmin, x, st0, ct0, ans, gauss, f;
120 int nit1 = nit + 1;
121 std::vector<double> xnit(nit1);
122 int ii, k;
123
124 double dthet = theta0 * pthet / 100.0;
125 dthet = dthet / 2.35;
126 double dthetr = dthet * 2.51;
127 tmax = theta0 + dthet * 3;
128 tmin = theta0 - dthet * 3;
129 double dt = (tmax - tmin) / nit;
130 for (int i = 0; i < nit1; i++) {
131 xnit[i] = tmin + dt * i;
132 }
133
134 // Could be parallelized at this point
135
136 for (size_t j = 0; j < nData; ++j) {
137 double lambda = 4 * M_PI * sin(theta0) / xValues[j];
138 cy[j] = 0.0;
139 double tl = lambda * lambda;
140 double tlc = 8.0 * M_PI * M_PI / tl;
141 double con = tl / (2.0 * M_PI);
142
143 for (k = 0; k < m_nlayer + 2; k++) {
144 rnfn[k] = (1.0 - con * rnbn[k]);
145 }
146
147 for (ii = 0; ii < nit1; ii++) {
148 x = xnit[ii];
149 ct0 = cos(x);
150 st0 = sin(x);
151 pfn[0] = rnfn[0] * st0;
152
153 for (int i = 1; i < m_nlayer + 1; i++) {
154 rnf = (rnfn[i] * rnfn[i]) * cr;
155 rnf1 = (rnfn[0] * rnfn[0]) * cr;
156 pfn[i] = sqrt(rnf - (rnf1 * ct0 * ct0));
157 }
158
159 rnf = (rnfn[m_nlayer + 1] * rnfn[m_nlayer + 1]) * cr;
160 rnf1 = (rnfn[0] * rnfn[0]) * cr;
161 pfn[m_nlayer + 1] = sqrt(rnf - (rnf1 * ct0 * ct0));
162
163 for (int i = 1; i < m_nlayer + 1; i++) {
164 betan[i] = 2.0 * M_PI * dn[i] * pfn[i] / lambda;
165 }
166
167 a111 = cr;
168 if ((pfn[0] + pfn[1]) != c0) {
169 a12t = (pfn[0] - pfn[1]) / (pfn[0] + pfn[1]);
170 } else {
171 a12t = c0;
172 }
173 a112 = a12t * exp(-tlc * zbb[0] * pfn[0] * pfn[1]);
174 a121 = a112;
175 a122 = cr;
176
177 for (int i = 1; i < m_nlayer + 1; i++) {
178 btm = betan[i] * ci;
179 btm1 = -1.0 * betan[i] * ci;
180 if (imag(btm1) == 0.0)
181 btm = betan[i] * 2.0 * ci;
182 cbtm = exp(btm);
183 if (real(btm1) == 0.0)
184 cbtm1 = exp(btm1);
185 if (imag(btm1) == 0.0)
186 cbtm1 = 1.0 * cr;
187 a211 = cbtm;
188 if ((pfn[i] + pfn[i + 1]) != c0)
189 a22t = (pfn[i] - pfn[i + 1]) / (pfn[i] + pfn[i + 1]);
190 if ((pfn[i] + pfn[i + 1]) == c0)
191 a22t = c0;
192 a22t = a22t * exp(-tlc * zbb[i] * pfn[i] * pfn[i + 1]);
193 a212 = a22t * cbtm;
194 a221 = a22t * cbtm1;
195 a222 = cbtm1;
196
197 a311 = a111 * a211 + a112 * a221;
198 a312 = a111 * a212 + a112 * a222;
199 a321 = a121 * a211 + a122 * a221;
200 a322 = a121 * a212 + a122 * a222;
201
202 a111 = a311;
203 a112 = a312;
204 a121 = a321;
205 a122 = a322;
206 }
207
208 ac1 = a121;
209 ac2 = conj(ac1);
210 ac3 = a111;
211 ac4 = conj(ac3);
212 if (ac3 == c0) {
213 ans = 1.0;
214 } else {
215 ans = real((ac1 * ac2) / (ac3 * ac4));
216 }
217
218 gauss = (1.0 / dthetr) * exp(-0.5 * ((theta0 - x) / dthet) * (theta0 - x) / dthet);
219 f = ans * gauss;
220 cy[j] = cy[j] + f * dt;
221 }
222 out[j] = (cy[j] + bgd) * scalefac;
223 // g_log.information() << "cy[j]" << cy[j] << "\n";
224 }
225 // End of parallelized section
226}
227
228//----------------------------------------------------------------------------------------------
236void ReflectivityMulf::setAttribute(const std::string &attName, const API::IFunction::Attribute &att) {
237 storeAttributeValue(attName, att);
238 if (attName == "nlayer") {
239 m_nlayer = att.asInt();
240 if (m_nlayer < 0) {
241 throw std::invalid_argument("ReflectivityMulf: reflectivity number of "
242 "layers cannot be negative.");
243 }
244
245 // get the current parameter values using the old value of nlayer
246 vector<double> coeff(m_nlayer * 3 + 7, 0.0);
247 for (int i = 0; i < 7; ++i)
248 coeff[i] = getParameter(i);
249
250 if (m_nlayer > m_nlayer_old) {
251 for (int i = 0; i < m_nlayer_old; ++i) {
252 coeff[7 + i * 3] = getParameter(7 + i * 3);
253 coeff[8 + i * 3] = getParameter(8 + i * 3);
254 coeff[9 + i * 3] = getParameter(9 + i * 3);
255 }
256 // fill in the missing values if m_nlayer>m_nlayer_old
257 for (int i = m_nlayer_old; i < m_nlayer; ++i) {
258 coeff[7 + i * 3] = 0.0;
259 coeff[8 + i * 3] = 0.0;
260 coeff[9 + i * 3] = 0.0;
261 }
262 } else {
263 for (int i = 0; i < m_nlayer; ++i) {
264 coeff[7 + i * 3] = getParameter(7 + i * 3);
265 coeff[8 + i * 3] = getParameter(8 + i * 3);
266 coeff[9 + i * 3] = getParameter(9 + i * 3);
267 }
268 }
269 // now update the value of nlayer_old
271
272 // set the reflectivity layers
273 if (m_nlayer >= 0) {
275 }
276
277 declareParameter("Theta", coeff[0]);
278 declareParameter("ScaleFactor", coeff[1]);
279 declareParameter("AirSLD", coeff[2]);
280 declareParameter("BulkSLD", coeff[3]);
281 declareParameter("Roughness", coeff[4]);
282 declareParameter("BackGround", coeff[5]);
283 declareParameter("Resolution", coeff[6]);
284
285 for (int i = 0; i < m_nlayer; ++i) {
286 std::string parName = "SLD_Layer" + std::to_string(i);
287 declareParameter(parName, coeff[7 + i * 3]);
288 parName = "d_Layer" + std::to_string(i);
289 declareParameter(parName, coeff[8 + i * 3]);
290 parName = "Rough_Layer" + std::to_string(i);
291 declareParameter(parName, coeff[9 + i * 3]);
292 }
293 }
294}
295
296} // namespace Mantid::CurveFitting::Functions
const std::vector< double > * lambda
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
int asInt() const
Returns int value if attribute is a int, throws exception otherwise.
Definition: IFunction.cpp:726
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition: IFunction.cpp:1418
void storeAttributeValue(const std::string &name, const API::IFunction::Attribute &value)
Store an attribute's value.
Definition: IFunction.cpp:1464
void clearAllParameters()
Nonvirtual member which removes all declared parameters.
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.
ReflectivityMulf : Calculate the ReflectivityMulf from a simple layer model.
void function1D(double *out, const double *xValues, const size_t nData) const override
Function to calcualte reflectivity.
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
void init() override
Function initialization. Declare function parameters in this method.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)