Keren.cpp

Go to the documentation of this file.

// 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)
 
 
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)");
}
 
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);
}
 
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);
}
 
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]);
  }
}
 
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));
}
 
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