CrystalFieldFunction.cpp

    propFun.setHamiltonian(hamiltonian, nre);
  } else {
    throw std::runtime_error("Physical property type not understood: " +
                             propName);
  }
}

/// Update m_spectrum function.
void CrystalFieldFunction::updateTargetFunction() const {
  if (!m_target) {
    buildTargetFunction();
    return;
  }
  m_dirtyTarget = false;
  if (isMultiSite()) {
    updateMultiSite();
  } else {
    updateSingleSite();
  }
  m_target->checkFunction();
}

/// Update the target function in a single site case.
void CrystalFieldFunction::updateSingleSite() const {
  if (isMultiSpectrum()) {
    updateSingleSiteMultiSpectrum();
  } else {
    updateSingleSiteSingleSpectrum();
  }
}

/// Update the target function in a multi site case.
void CrystalFieldFunction::updateMultiSite() const {
  if (isMultiSpectrum()) {
    updateMultiSiteMultiSpectrum();
  } else {
    updateMultiSiteSingleSpectrum();
  }
}

/// Update the target function in a single site - single spectrum case.
void CrystalFieldFunction::updateSingleSiteSingleSpectrum() const {
  auto fwhmVariation = m_control.getAttribute("FWHMVariation").asDouble();
  auto peakShape = m_control.getAttribute("PeakShape").asString();
  bool fixAllPeaks = m_control.getAttribute("FixAllPeaks").asBool();
  auto xVec = m_control.getAttribute("FWHMX").asVector();
  auto yVec = m_control.getAttribute("FWHMX").asVector();
  auto &FWHMs = m_control.FWHMs();
  auto defaultFWHM = FWHMs.empty() ? 0.0 : FWHMs[0];
  size_t indexShift = hasBackground() ? 1 : 0;

  FunctionDomainGeneral domain;
  FunctionValues values;
  m_source->function(domain, values);
  m_target->setAttributeValue("NumDeriv", true);
  auto &spectrum = dynamic_cast<CompositeFunction &>(*m_target);
  CrystalFieldUtils::updateSpectrumFunction(
      spectrum, peakShape, values, indexShift, xVec, yVec, fwhmVariation,
      defaultFWHM, fixAllPeaks);
}

/// Update the target function in a single site - multi spectrum case.
void CrystalFieldFunction::updateSingleSiteMultiSpectrum() const {
  DoubleFortranVector energies;
  ComplexFortranMatrix waveFunctions;
  ComplexFortranMatrix hamiltonian;
  ComplexFortranMatrix hamiltonianZeeman;
  int nre = 0;
  auto &peakCalculator = dynamic_cast<CrystalFieldPeaksBase &>(*m_source);
  peakCalculator.calculateEigenSystem(energies, waveFunctions, hamiltonian,
                                      hamiltonianZeeman, nre);
  hamiltonian += hamiltonianZeeman;
  size_t iFirst = hasBackground() ? 1 : 0;

  auto &fun = dynamic_cast<MultiDomainFunction &>(*m_target);
  auto &temperatures = m_control.temperatures();
  auto &FWHMs = m_control.FWHMs();
  for (size_t iSpec = 0; iSpec < temperatures.size(); ++iSpec) {
    updateSpectrum(*fun.getFunction(iSpec), nre, energies, waveFunctions,
                   temperatures[iSpec], FWHMs[iSpec], iSpec, iFirst);
  }

  for (auto &prop : m_mapPrefixes2PhysProps) {
    updatePhysprop(nre, energies, waveFunctions, hamiltonian, *prop.second);
  }
}

/// Update the target function in a multi site - single spectrum case.
void CrystalFieldFunction::updateMultiSiteSingleSpectrum() const {
  auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
  auto peakShape = getAttribute("PeakShape").asString();
  bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();
  auto xVec = m_control.getAttribute("FWHMX").asVector();
  auto yVec = m_control.getAttribute("FWHMX").asVector();
  auto &FWHMs = m_control.FWHMs();
  auto defaultFWHM = FWHMs.empty() ? 0.0 : FWHMs[0];

  size_t spectrumIndexShift = hasBackground() ? 1 : 0;
  auto &compSource = compositeSource();
  for (size_t ionIndex = 0; ionIndex < compSource.nFunctions(); ++ionIndex) {
    FunctionDomainGeneral domain;
    FunctionValues values;
    compSource.getFunction(ionIndex)->function(domain, values);

    auto &ionSpectrum = dynamic_cast<CompositeFunction &>(
        *m_target->getFunction(ionIndex + spectrumIndexShift));
    CrystalFieldUtils::updateSpectrumFunction(ionSpectrum, peakShape, values, 0,
                                              xVec, yVec, fwhmVariation,
                                              defaultFWHM, fixAllPeaks);
  }
}

/// Update the target function in a multi site - multi spectrum case.
void CrystalFieldFunction::updateMultiSiteMultiSpectrum() const {
  auto &compSource = compositeSource();
  for (size_t ionIndex = 0; ionIndex < compSource.nFunctions(); ++ionIndex) {
    DoubleFortranVector energies;
    ComplexFortranMatrix waveFunctions;
    ComplexFortranMatrix hamiltonian;
    ComplexFortranMatrix hamiltonianZeeman;
    int nre = 0;
    auto &peakCalculator = dynamic_cast<CrystalFieldPeaksBase &>(
        *compSource.getFunction(ionIndex));
    peakCalculator.calculateEigenSystem(energies, waveFunctions, hamiltonian,
                                        hamiltonianZeeman, nre);
    hamiltonian += hamiltonianZeeman;
    size_t iFirst = ionIndex == 0 && hasBackground() ? 1 : 0;

    auto &temperatures = m_control.temperatures();
    auto &FWHMs = m_control.FWHMs();
    for (size_t iSpec = 0; iSpec < temperatures.size(); ++iSpec) {
      auto &spectrum =
          dynamic_cast<CompositeFunction &>(*m_target->getFunction(iSpec));
      auto &ionSpectrum =
          dynamic_cast<CompositeFunction &>(*spectrum.getFunction(ionIndex));
      updateSpectrum(ionSpectrum, nre, energies, waveFunctions,
                     temperatures[iSpec], FWHMs[iSpec], iSpec, iFirst);
    }

    std::string prefix("ion");
    prefix.append(std::to_string(ionIndex)).append(".");
    auto prefixSize = prefix.size();
    for (auto prop : m_mapPrefixes2PhysProps) {
      if (prop.first.substr(0, prefixSize) == prefix) {
        updatePhysprop(nre, energies, waveFunctions, hamiltonian, *prop.second);
      }
    }
  }
}

/// Update a function for a single spectrum.
/// @param spectrum :: A Spectrum function to update.
/// @param nre :: An id of the ion.
/// @param energies :: A vector with energies.
/// @param waveFunctions :: A matrix with wave functions.
/// @param temperature :: A temperature of the spectrum.
/// @param fwhm :: A full width at half maximum to set to each peak.
/// @param iSpec :: An index of the created spectrum in m_target composite
/// function.
/// @param iFirst :: An index of the first peak in spectrum composite function.
void CrystalFieldFunction::updateSpectrum(
    API::IFunction &spectrum, int nre, const DoubleFortranVector &energies,
    const ComplexFortranMatrix &waveFunctions, double temperature, double fwhm,
    size_t iSpec, size_t iFirst) const {
  const auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
  const auto peakShape = getAttribute("PeakShape").asString();
  const bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();
  auto xVec = m_control.getFunction(iSpec)->getAttribute("FWHMX").asVector();
  auto yVec = m_control.getFunction(iSpec)->getAttribute("FWHMX").asVector();
  auto intensityScaling =
      m_control.getFunction(iSpec)->getParameter("IntensityScaling");
  FunctionValues values;
  calcExcitations(nre, energies, waveFunctions, temperature, values,
                  intensityScaling);
  auto &composite = dynamic_cast<API::CompositeFunction &>(spectrum);
  CrystalFieldUtils::updateSpectrumFunction(composite, peakShape, values,
                                            iFirst, xVec, yVec, fwhmVariation,
                                            fwhm, fixAllPeaks);
}

/// Make maps between parameter names and indices
void CrystalFieldFunction::makeMaps() const {
  m_mapNames2Indices.clear();
  m_mapIndices2Names.resize(nParams());
  if (isMultiSite()) {
    if (isMultiSpectrum()) {
      makeMapsMultiSiteMultiSpectrum();
    } else {
      makeMapsMultiSiteSingleSpectrum();
    }
  } else {
    if (isMultiSpectrum()) {
      makeMapsSingleSiteMultiSpectrum();
    } else {
      makeMapsSingleSiteSingleSpectrum();
    }
  }
}

/// Make parameter names from names of a function and map them to indices
/// @param fun :: A function to get parameter names from.
/// @param iFirst :: An index that maps to the first parameter of fun.
/// @param prefix :: A prefix to add to all parameters
size_t
CrystalFieldFunction::makeMapsForFunction(const IFunction &fun, size_t iFirst,
                                          const std::string &prefix) const {
  auto n = fun.nParams();
  for (size_t i = 0; i < n; ++i) {
    size_t j = i + iFirst;
    auto name(prefix);
    name.append(fun.parameterName(i));
    m_mapNames2Indices[name] = j;
    m_mapIndices2Names[j] = name;
  }
  return n;
}

/// Parameter-index map for single-site single-spectrum
void CrystalFieldFunction::makeMapsSingleSiteSingleSpectrum() const {
  size_t i = makeMapsForFunction(*m_source, 0, "");

  size_t peakIndex = 0;
  // If there is a background it's the first function in m_target
  if (hasBackground()) {
    auto &background = *m_target->getFunction(0);
    i += makeMapsForFunction(background, i, BACKGROUND_PREFIX + ".");
    peakIndex = 1;
  }
  // All other functions are peaks.
  for (size_t ip = peakIndex; ip < m_target->nFunctions(); ++ip) {
    std::string prefix(PEAK_PREFIX);
    prefix.append(std::to_string(ip - peakIndex)).append(".");
    i += makeMapsForFunction(*m_target->getFunction(ip), i, prefix);
  }
}

/// Parameter-index map for single-site multi-spectrum
void CrystalFieldFunction::makeMapsSingleSiteMultiSpectrum() const {
  size_t i = 0;
  // Intensity scalings for each spectrum
  for (size_t j = 0; j < m_control.nFunctions(); ++j) {
    std::string prefix(SPECTRUM_PREFIX);
    prefix.append(std::to_string(j)).append(".");
    i += makeMapsForFunction(*m_control.getFunction(j), i, prefix);
  }
  // Crystal field parameters
  i += makeMapsForFunction(*m_source, i, "");

  size_t peakIndex = 0;
  for (size_t iSpec = 0; iSpec < m_target->nFunctions(); ++iSpec) {
    if (auto spectrum = dynamic_cast<const CompositeFunction *>(
            m_target->getFunction(iSpec).get())) {
      // This is a normal spectrum
      std::string spectrumPrefix(SPECTRUM_PREFIX);
      spectrumPrefix.append(std::to_string(iSpec)).append(".");
      // If there is a background it's the first function in spectrum
      if (hasBackground()) {
        auto &background = *spectrum->getFunction(0);
        i += makeMapsForFunction(background, i,
                                 spectrumPrefix + BACKGROUND_PREFIX + ".");
        peakIndex = 1;
      }
      // All other functions are peaks.
      for (size_t ip = peakIndex; ip < spectrum->nFunctions(); ++ip) {
        std::string prefix(spectrumPrefix);
        prefix.append(PEAK_PREFIX)
            .append(std::to_string(ip - peakIndex))
            .append(".");
        i += makeMapsForFunction(*spectrum->getFunction(ip), i, prefix);
      }
    } else {
      // This is a physical property function
      std::string prefix(m_control.physProps()[iSpec - nSpectra()]);
      prefix.append(".");
      i += makeMapsForFunction(*m_target->getFunction(iSpec), i, prefix);
    }
  }
}

/// Parameter-index map for multi-site single-spectrum
void CrystalFieldFunction::makeMapsMultiSiteSingleSpectrum() const {
  size_t i = 0;
  // Intensity scalings for each ion
  auto &crystalField = compositeSource();
  for (size_t ion = 0; ion < crystalField.nFunctions(); ++ion) {
    std::string prefix(ION_PREFIX);
    prefix.append(std::to_string(ion)).append(".");
    i += makeMapsForFunction(*crystalField.getFunction(ion), i, prefix);
  }
  // Spectrum split into an optional background and groups of peaks for
  // each ion
  size_t ionIndex = 0;
  // If there is a background it's the first function in spectrum
  if (hasBackground()) {
    auto &background = *m_target->getFunction(0);
    i += makeMapsForFunction(background, i, BACKGROUND_PREFIX + ".");
    ionIndex = 1;
  }
  // All other functions are ion spectra.
  for (size_t ion = ionIndex; ion < m_target->nFunctions(); ++ion) {
    std::string ionPrefix(ION_PREFIX);
    ionPrefix.append(std::to_string(ion - ionIndex)).append(".");
    // All other functions are peaks.
    auto &spectrum =
        dynamic_cast<const CompositeFunction &>(*m_target->getFunction(ion));
    for (size_t ip = 0; ip < spectrum.nFunctions(); ++ip) {
      std::string prefix(ionPrefix);
      prefix.append(PEAK_PREFIX).append(std::to_string(ip)).append(".");
      i += makeMapsForFunction(*spectrum.getFunction(ip), i, prefix);
    }
  }
}

/// Parameter-index map for multi-site multi-spectrum
void CrystalFieldFunction::makeMapsMultiSiteMultiSpectrum() const {
  size_t i = 0;
  // Intensity scalings for each spectrum
  for (size_t j = 0; j < m_control.nFunctions(); ++j) {
    std::string prefix(SPECTRUM_PREFIX);
    prefix.append(std::to_string(j)).append(".");
    i += makeMapsForFunction(*m_control.getFunction(j), i, prefix);
  }
  // Intensity scalings for each ion
  auto &crystalField = compositeSource();
  for (size_t ion = 0; ion < crystalField.nFunctions(); ++ion) {
    std::string prefix(ION_PREFIX);
    prefix.append(std::to_string(ion)).append(".");
    i += makeMapsForFunction(*crystalField.getFunction(ion), i, prefix);
  }

  // The spectra (background and peak) parameters
  for (size_t iSpec = 0; iSpec < nSpectra(); ++iSpec) {
    auto &spectrum =
        dynamic_cast<const CompositeFunction &>(*m_target->getFunction(iSpec));
    std::string spectrumPrefix(SPECTRUM_PREFIX);
    spectrumPrefix.append(std::to_string(iSpec)).append(".");

    // All other functions are ion spectra.
    for (size_t ion = 0; ion < crystalField.nFunctions(); ++ion) {
      auto &ionSpectrum =
          dynamic_cast<const CompositeFunction &>(*spectrum.getFunction(ion));
      size_t peakIndex = 0;
      if (ion == 0 && hasBackground()) {
        peakIndex = 1;
        std::string prefix(spectrumPrefix);
        prefix.append(BACKGROUND_PREFIX).append(".");
        i += makeMapsForFunction(*ionSpectrum.getFunction(0), i, prefix);
      }
      std::string ionPrefix(ION_PREFIX);
      ionPrefix.append(std::to_string(ion)).append(".").append(spectrumPrefix);
      // All other functions are peaks.
      for (size_t ip = peakIndex; ip < ionSpectrum.nFunctions(); ++ip) {
        std::string prefix(ionPrefix);
        prefix.append(PEAK_PREFIX)
            .append(std::to_string(ip - peakIndex))
            .append(".");
        i += makeMapsForFunction(*ionSpectrum.getFunction(ip), i, prefix);
      }
    }
  }
  // The phys prop parameters
  for (size_t iSpec = nSpectra(); iSpec < m_target->nFunctions(); ++iSpec) {
    auto &spectrum =
        dynamic_cast<const CompositeFunction &>(*m_target->getFunction(iSpec));
    std::string physPropPrefix(spectrum.getFunction(0)->name());
    physPropPrefix.append(".");
    for (size_t ion = 0; ion < crystalField.nFunctions(); ++ion) {
      std::string prefix(ION_PREFIX);
      prefix.append(std::to_string(ion)).append(".").append(physPropPrefix);
      i += makeMapsForFunction(*spectrum.getFunction(ion), i, prefix);
    }
  }
}

} // namespace Functions
} // namespace CurveFitting
} // namespace Mantid