IPowderDiffPeakFunction.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 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IPowderDiffPeakFunction.h"
#include "MantidAPI/Jacobian.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/Exception.h"
 
#include <boost/lexical_cast.hpp>
#include <cmath>
 
const double IGNOREDCHANGE = 1.0E-9;
 
namespace Mantid::API {
 
int IPowderDiffPeakFunction::s_peakRadius = 5;
 
//----------------------------------------------------------------------------------------------
IPowderDiffPeakFunction::IPowderDiffPeakFunction()
    : m_centre(0.), m_dcentre(0.), m_fwhm(0.), m_hasNewParameterValue(false), m_cellParamValueChanged(false),
      m_sortedProfileParameterNames(), m_unitCell(), m_unitCellSize(0.), m_parameterValid(false), mH(0), mK(0), mL(0),
      mHKLSet(false), LATTICEINDEX(9999), HEIGHTINDEX(9999) {
  // Set peak's radius from configuration
  auto peakRadius = Kernel::ConfigService::Instance().getValue<int>("curvefitting.peakRadius");
  if (peakRadius.is_initialized() && peakRadius.get() != s_peakRadius) {
    setPeakRadius(peakRadius.get());
  }
}
 
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::setParameter(size_t i, const double &value, bool explicitlySet) {
  double origparamvalue = getParameter(i);
  if (fabs(origparamvalue - value) > IGNOREDCHANGE) {
    m_hasNewParameterValue = true;
  }
  ParamFunction::setParameter(i, value, explicitlySet);
}
 
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::setParameter(const std::string &name, const double &value, bool explicitlySet) {
  double origparamvalue = getParameter(name);
  if (fabs(origparamvalue - value) > IGNOREDCHANGE) {
    m_hasNewParameterValue = true;
  }
  ParamFunction::setParameter(name, value, explicitlySet);
}
 
//----------------------------------------------------------------------------------------------
double IPowderDiffPeakFunction::centre() const {
  // Re-calcualte peak parameters if required
  if (m_hasNewParameterValue)
    calculateParameters(false);
 
  return m_centre;
}
 
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::setHeight(const double h) { setParameter(HEIGHTINDEX, h); }
 
//----------------------------------------------------------------------------------------------
double IPowderDiffPeakFunction::height() const {
  double height = this->getParameter(HEIGHTINDEX);
  return height;
}
 
//----------------------------------------------------------------------------------------------
double IPowderDiffPeakFunction::fwhm() const {
  if (m_hasNewParameterValue)
    calculateParameters(false);
 
  return m_fwhm;
}
 
//----------------------------------------------------------------------------------------------
double IPowderDiffPeakFunction::getMaximumValue(const std::vector<double> &xValues, size_t &indexmax) const {
  // Calculate function with given data points
  std::vector<double> out(xValues.size(), 0.);
  function(out, xValues);
 
  // Find out the maximum value
  double max = -DBL_MAX;
  indexmax = 0;
  for (size_t i = 0; i < out.size(); ++i) {
    if (out[i] > max) {
      max = out[i];
      indexmax = i;
    }
  }
 
  return max;
}
 
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::setMillerIndex(int h, int k, int l) {
  // Check validity and set flag
  if (mHKLSet) {
    // Throw exception if tried to reset the miller index
    std::stringstream errss;
    errss << "Profile function " << name() << "cannot have (HKL) reset.";
    throw std::runtime_error(errss.str());
  } else {
    // Set flag
    mHKLSet = true;
  }
 
  // Set value
  mH = static_cast<int>(h);
  mK = static_cast<int>(k);
  mL = static_cast<int>(l);
 
  // Check value valid or not
  if (mH * mH + mK * mK + mL * mL < 1.0E-8) {
    std::stringstream errmsg;
    errmsg << "H = K = L = 0 is not allowed";
    throw std::invalid_argument(errmsg.str());
  }
}
 
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::getMillerIndex(int &h, int &k, int &l) {
  h = static_cast<int>(mH);
  k = static_cast<int>(mK);
  l = static_cast<int>(mL);
}
 
//----------------------------------------------------------------------------------------------
void IPowderDiffPeakFunction::setPeakRadius(const int &r) {
  if (r > 0) {
    s_peakRadius = r;
 
    // std::string setting = boost::lexical_cast<std::string>(r);
    // Kernel::ConfigService::Instance().setString("curvefitting.peakRadius",setting);
  }
}
 
//----------------------------------------------------------------------------------------------
bool IPowderDiffPeakFunction::hasProfileParameter(std::string paramname) {
  auto candname = lower_bound(m_sortedProfileParameterNames.begin(), m_sortedProfileParameterNames.end(), paramname);
  if (candname == m_sortedProfileParameterNames.end())
    return false;
 
  return *candname != paramname;
}
 
//-------------------------  External Functions
//---------------------------------------------------
// FIXME : This is the same function used for ThermalNeutron... ...
std::complex<double> E1(std::complex<double> z) {
  std::complex<double> exp_e1;
 
  double rz = real(z);
  double az = abs(z);
 
  if (fabs(az) < 1.0E-8) {
    // If z = 0, then the result is infinity... diverge!
    std::complex<double> r(1.0E300, 0.0);
    exp_e1 = r;
  } else if (az <= 10.0 || (rz < 0.0 && az < 20.0)) {
    // Some interesting region, equal to integrate to infinity, converged
    // cout << "[DB] Type 1\n";
 
    std::complex<double> r(1.0, 0.0);
    exp_e1 = r;
    std::complex<double> cr = r;
 
    for (size_t k = 1; k <= 150; ++k) {
      auto dk = double(k);
      cr = -cr * dk * z / ((dk + 1.0) * (dk + 1.0));
      exp_e1 += cr;
      if (abs(cr) < abs(exp_e1) * 1.0E-15) {
        // cr is converged to zero
        break;
      }
    } // ENDFOR k
 
    const double el = 0.5772156649015328;
    exp_e1 = -el - log(z) + (z * exp_e1);
  } else {
    // Rest of the region
    std::complex<double> ct0(0.0, 0.0);
    for (int k = 120; k > 0; --k) {
      std::complex<double> dk(double(k), 0.0);
      ct0 = dk / (10.0 + dk / (z + ct0));
    } // ENDFOR k
 
    exp_e1 = 1.0 / (z + ct0);
    exp_e1 = exp_e1 * exp(-z);
    if (rz < 0.0 && fabs(imag(z)) < 1.0E-10) {
      std::complex<double> u(0.0, 1.0);
      exp_e1 = exp_e1 - (M_PI * u);
    }
  }
 
  return exp_e1;
}
 
} // namespace Mantid::API