diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 9543836..98a3668 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -255,12 +255,12 @@ namespace serif::composition { return true; } - bool Composition::isValidSymbol(const std::string& symbol) const { - return chemSpecies::species.count(symbol) > 0; + bool Composition::isValidSymbol(const std::string& symbol) { + return chemSpecies::species.contains(symbol); } double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { - if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { + if (!m_registeredSymbols.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); throw std::runtime_error("Symbol is not registered."); } @@ -276,7 +276,7 @@ namespace serif::composition { } m_finalized = false; - double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); + const double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); m_compositions.at(symbol).setMassFraction(mass_fraction); return old_mass_fraction; @@ -432,9 +432,9 @@ namespace serif::composition { Composition mixedComposition(mixedSymbols); for (const auto& symbol : mixedSymbols) { - double thisMassFrac, otherMassFrac = 0.0; + double otherMassFrac = 0.0; - thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; + const double thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); @@ -449,7 +449,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -462,7 +462,7 @@ namespace serif::composition { std::unordered_map Composition::getMassFraction() const { std::unordered_map mass_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { mass_fractions[symbol] = getMassFraction(symbol); } return mass_fractions; @@ -474,7 +474,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -487,7 +487,7 @@ namespace serif::composition { std::unordered_map Composition::getNumberFraction() const { std::unordered_map number_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { number_fractions[symbol] = getNumberFraction(symbol); } return number_fractions; @@ -498,7 +498,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -527,40 +527,29 @@ namespace serif::composition { throw std::runtime_error("Composition not finalized. Cannot retrieve mean atomic mass number."); } - double mean_A = 0.0; + double zSum = 0.0; // Loop through all registered species in the composition. for (const auto &val: m_compositions | std::views::values) { - const CompositionEntry& entry = val; - const chemSpecies::Species& species = entry.isotope(); - - const double mass_fraction = entry.mass_fraction(); - const double particle_mass_g = species.mass(); - const int mass_number = species.a(); - - // Avoid division by zero, though a valid species should have a positive mass. - if (particle_mass_g > 0) { - // Calculate the number fraction for this species. - const double number_fraction = (mass_fraction / particle_mass_g) * m_meanParticleMass; - mean_A += number_fraction * mass_number; - } + zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); } + const double mean_A = m_meanParticleMass * zSum; return mean_A; } Composition Composition::subset(const std::vector& symbols, std::string method) const { - std::array methods = {"norm", "none"}; + const std::array methods = {"norm", "none"}; - if (std::find(methods.begin(), methods.end(), method) == methods.end()) { - std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; + if (std::ranges::find(methods, method) == methods.end()) { + const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method); throw std::runtime_error(errorMessage); } Composition subsetComposition; for (const auto& symbol : symbols) { - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } else { @@ -569,7 +558,7 @@ namespace serif::composition { subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); } if (method == "norm") { - bool isNorm = subsetComposition.finalize(true); + const bool isNorm = subsetComposition.finalize(true); if (!isNorm) { LOG_ERROR(m_logger, "Subset composition is invalid."); throw std::runtime_error("Subset composition is invalid."); @@ -578,14 +567,14 @@ namespace serif::composition { return subsetComposition; } - void Composition::setCompositionMode(bool massFracMode) { + void Composition::setCompositionMode(const bool massFracMode) { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized."); } bool okay = true; - for (auto& [_, entry] : m_compositions) { + for (auto &entry: m_compositions | std::views::values) { if (massFracMode) { okay = entry.setMassFracMode(m_meanParticleMass); } else { @@ -599,6 +588,53 @@ namespace serif::composition { m_massFracMode = massFracMode; } + CanonicalComposition Composition::getCanonicalComposition(bool harsh) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + CanonicalComposition canonicalComposition; + constexpr std::array canonicalH = { + "H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7" + }; + constexpr std::array canonicalHe = { + "He-3", "He-4", "He-5", "He-6", "He-7", "He-8", "He-9", "He-10" + }; + for (const auto& symbol : canonicalH) { + if (hasSymbol(symbol)) { + canonicalComposition.X += getMassFraction(symbol); + } + } + for (const auto& symbol : canonicalHe) { + if (hasSymbol(symbol)) { + canonicalComposition.Y += getMassFraction(symbol); + } + } + + for (const auto& symbol : getRegisteredSymbols()) { + const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH); + const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe); + + if (isHSymbol || isHeSymbol) { + continue; // Skip canonical H and He symbols + } + + canonicalComposition.Z += getMassFraction(symbol); + } + + const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y); + if (std::abs(Z - canonicalComposition.Z) > 1e-6) { + if (!harsh) { + LOG_WARNING(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + } + else { + LOG_ERROR(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + throw std::runtime_error("Validation composition Z (X-Y = " + std::to_string(Z) + ") is different than canonical composition Z (" + std::to_string(canonicalComposition.Z) + ") (∑a_i where a_i != H/He)."); + } + } + return canonicalComposition; + } + bool Composition::hasSymbol(const std::string& symbol) const { return m_compositions.count(symbol) > 0; } diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 457f96d..861cc47 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -34,6 +34,20 @@ #include "atomicSpecies.h" namespace serif::composition { + struct CanonicalComposition { + double X = 0.0; ///< Mass fraction of Hydrogen. + double Y = 0.0; ///< Mass fraction of Helium. + double Z = 0.0; ///< Mass fraction of Metals. + + friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) { + os << ""; + return os; + } + }; + /** * @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use. */ @@ -220,7 +234,7 @@ namespace serif::composition { * @param symbol The symbol to check. * @return True if the symbol is valid, false otherwise. */ - bool isValidSymbol(const std::string& symbol) const; + static bool isValidSymbol(const std::string& symbol); /** * @brief Checks if the given mass fractions are valid. @@ -466,6 +480,16 @@ namespace serif::composition { */ void setCompositionMode(bool massFracMode); + /** + * @brief Gets the current canonical composition (X, Y, Z). + * @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements. + * @return True if mass fraction mode, false if number fraction mode. + * + * @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed. + * @throws std::runtime_error if harsh is true and the canonical composition is not valid. + */ + [[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const; + /** * @brief Overloaded output stream operator for Composition. * @param os The output stream.