Mantid
Loading...
Searching...
No Matches
InelasticIsoRotDiff.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 +
7// Mantid Coding standars <http://www.mantidproject.org/Coding_Standards>
8// Main Module Header
10// Mantid Headers from the same project
12// Mantid headers from other projects
14#include "MantidAPI/IFunction.h"
15
16// third party library headers
17#include <boost/math/special_functions/bessel.hpp>
18// standard library headers
19#include <cmath>
20#include <limits>
21
23
24namespace {
25Mantid::Kernel::Logger g_log("IsoRotDiff");
26}
27
29
30DECLARE_FUNCTION(InelasticIsoRotDiff)
31
32
36 this->declareParameter("Height", 1.0, "scaling factor");
37 this->declareParameter("Radius", 0.98, "radius of rotation (Angstroms)");
38 this->declareParameter("Tau", 10.0, "Relaxation time, inverse of the rotational diffusion coefficient (ps)");
39 this->declareParameter("Centre", 0.0, "Shift along the X-axis");
40
41 this->declareAttribute("Q", API::IFunction::Attribute(0.3));
42 this->declareAttribute("N", API::IFunction::Attribute(25));
43}
44
49 // Ensure positive values for Height, Radius, and Diffusion constant
50 auto HeightConstraint = std::make_unique<BConstraint>(this, "Height", std::numeric_limits<double>::epsilon(), true);
51 this->addConstraint(std::move(HeightConstraint));
52 auto RadiusConstraint = std::make_unique<BConstraint>(this, "Radius", std::numeric_limits<double>::epsilon(), true);
53 this->addConstraint(std::move(RadiusConstraint));
54 auto DiffusionConstraint = std::make_unique<BConstraint>(this, "Tau", std::numeric_limits<double>::epsilon(), true);
55 this->addConstraint(std::move(DiffusionConstraint));
56}
57
64void InelasticIsoRotDiff::function1D(double *out, const double *xValues, const size_t nData) const {
65 double hbar(0.658211626); // ps*meV
66 auto H = this->getParameter("Height");
67 auto R = this->getParameter("Radius");
68 auto T = this->getParameter("Tau");
69 auto C = this->getParameter("Centre");
70 auto Q = this->getAttribute("Q").asDouble();
71 auto N = static_cast<size_t>(this->getAttribute("N").asInt()); // Number of Lorentzians
72
73 // Penalize negative parameters
74 if (R < std::numeric_limits<double>::epsilon()) {
75 for (size_t j = 0; j < nData; j++) {
76 out[j] = std::numeric_limits<double>::infinity();
77 }
78 return;
79 }
80
81 // Lorentzian intensities and HWHM
82 std::vector<double> al(N);
83 std::vector<double> HWHM(N);
84 for (size_t i = 0; i < N; i++) {
85 auto l = static_cast<unsigned int>(i + 1); // avoid annoying warnings from implicit type conversion
86 auto ld = static_cast<double>(l); // avoid annoying warnings from implicit type conversion
87 al[i] = (2 * ld + 1) * pow(boost::math::sph_bessel(l, Q * R), 2);
88 HWHM[i] = ld * (ld + 1) * hbar / T;
89 }
90
91 for (size_t j = 0; j < nData; j++) {
92 out[j] = 0.0;
93 auto E = xValues[j] - C;
94 for (size_t i = 0; i < N; i++) {
95 auto G = HWHM[i];
96 out[j] += H * al[i] * G / (G * G + E * E) / M_PI;
97 }
98 }
99}
100
101} // namespace Mantid::CurveFitting::Functions
#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
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition: IFunction.cpp:1394
virtual void addConstraint(std::unique_ptr< IConstraint > ic)
Add a constraint to function.
Definition: IFunction.cpp:376
double getParameter(size_t i) const override
Get i-th parameter.
A boundary constraint is designed to be used to set either upper or lower (or both) boundaries on a s...
Inelastic part of the IsoRotDiff function.
void init() override
overwrite IFunction base class methods
void function1D(double *out, const double *xValues, const size_t nData) const override
Calculate function values on an energy domain.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52