FickDiffusionSQE.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/FickDiffusionSQE.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/Jacobian.h"
#include <cmath>
#include <limits>
 
namespace {
Mantid::Kernel::Logger g_log("FickDiffusionSQE");
}
 
namespace Mantid::CurveFitting::Functions {
 
DECLARE_FUNCTION(FickDiffusionSQE)
 
 
void FickDiffusionSQE::declareParameters() {
  this->declareParameter("Height", 1.0, "scaling factor");
  this->declareParameter("Centre", 0.0, "Shift along the X-axis");
  this->declareParameter("DiffCoeff", 2.3, "Diffusion coefficient (10^(-5)cm^2/s)");
}
 
void FickDiffusionSQE::function1D(double *out, const double *xValues, const size_t nData) const {
  auto H = this->getParameter("Height");
  auto C = this->getParameter("Centre");
  auto Q = this->getAttribute("Q").asDouble();
  auto D = this->getParameter("DiffCoeff");
 
  // Penalize negative parameters, just in case they show up
  // when calculating the numeric derivative
  if (H < std::numeric_limits<double>::epsilon() || D < std::numeric_limits<double>::epsilon()) {
    for (size_t j = 0; j < nData; j++) {
      out[j] = std::numeric_limits<double>::infinity();
    }
    return;
  }
 
  // Lorentzian intensities and HWHM
  D *= 0.10; // conversion from 10^{-5}cm^2/s to Angstrom^2/ps, the internal
             // units used
  auto G = D * Q * Q;
  for (size_t j = 0; j < nData; j++) {
    auto E = xValues[j] - C;
    out[j] += H * G / (G * G + E * E) / M_PI;
  }
}
 
void FickDiffusionSQE::functionDeriv1D(Mantid::API::Jacobian *jacobian, const double *xValues, const size_t nData) {
  const double deltaF{0.1}; // increase parameter by this fraction
  const size_t nParam = this->nParams();
  // cutoff defines the smallest change in the parameter when calculating the
  // numerical derivative
  std::map<std::string, double> cutoff;
  cutoff["DiffCoeff"] = 0.2; // 0.2x10^(-5)cm^2/s
  cutoff["Centre"] = 0.0001; // 0.1micro-eV
  std::vector<double> out(nData);
  this->applyTies();
  this->function1D(out.data(), xValues, nData);
 
  for (size_t iP = 0; iP < nParam; iP++) {
    if (this->isActive(iP)) {
      std::vector<double> derivative(nData);
      const double pVal = this->activeParameter(iP);
      const std::string pName = this->parameterName(iP);
      if (pName == "Height") {
        // exact derivative
        this->setActiveParameter(iP, 1.0);
        this->applyTies();
        this->function1D(derivative.data(), xValues, nData);
      } else {
        // numerical derivative
        double delta = cutoff[pName] > fabs(pVal * deltaF) ? cutoff[pName] : pVal * deltaF;
        this->setActiveParameter(iP, pVal + delta);
        this->applyTies();
        this->function1D(derivative.data(), xValues, nData);
        for (size_t i = 0; i < nData; i++) {
          derivative[i] = (derivative[i] - out[i]) / delta;
        }
      }
      this->setActiveParameter(iP, pVal); // restore the value of the parameter
      // fill the jacobian for this parameter
      for (size_t i = 0; i < nData; i++) {
        jacobian->set(i, iP, derivative[i]);
      }
    }
  }
}
 
} // namespace Mantid::CurveFitting::Functions