.. _program_listing_file_src_extinct_Extinct.hpp:

Program Listing for File Extinct.hpp
====================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_extinct_Extinct.hpp>` (``src/extinct/Extinct.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef EXTINCT_HPP
   #define EXTINCT_HPP
   
   #include "../interpolation/Interpolator1D.hpp"
   #include "../phot/FilterTabulated.hpp"
   #include "../utils/Constants.hpp"
   #include "../utils/HDF5Utils.hpp"
   #include "../utils/TOMLUtils.hpp"
   #include "hdf5.h" // NOLINT(misc-include-cleaner) -- see HDF5Utils.hpp's own comment on including hdf5.h wholesale
   #include <algorithm>
   #include <cmath>
   #include <cstddef>
   #include <filesystem>
   #include <stdexcept>
   #include <string>
   #include <toml.hpp>
   #include <utility>
   #include <vector>
   
   namespace io
   {
       class SimControls;
   } // namespace io
   
   namespace extinct
   {
       inline static const std::string defaultRegistry = // NOLINT(bugprone-throwing-static-initialization,cert-err58-cpp) -- built from fixed string literals, so the (theoretically throwing) path conversion can never actually throw here
           (std::filesystem::path("data") / std::filesystem::path("extinct")
           / std::filesystem::path("extinct.toml")); 
       inline static const std::string defaultVRegistry = // NOLINT(bugprone-throwing-static-initialization,cert-err58-cpp) -- see defaultRegistry above
           (std::filesystem::path("data") / std::filesystem::path("filters")
           / std::filesystem::path("V_filter.toml")); 
       class Extinct
       {
       public:
   
           Extinct(const std::string& extinctName,
               const io::SimControls& controls,
               const std::string& registryName = defaultRegistry) :
               controls_(controls)
           {
               loadCurveImpl(extinctName, registryName);
           }
   
           // Copyable (rebinding controls_ to the same referent) but not
           // assignable (controls_ can't be reseated), matching Specsyn's
           // own identical copy/move declarations exactly -- see its
           // comment. Never actually copied/assigned in practice: every
           // Extinct is held via unique_ptr (see SimControls's own
           // extinct_).
           Extinct(const Extinct&) = default;
           Extinct(Extinct&&) = default;
           auto operator=(const Extinct&) -> Extinct& = delete;
           auto operator=(Extinct&&) -> Extinct& = delete;
           ~Extinct() = default;
   
           void loadCurve(const std::string& extinctName,
               const std::string& registryName = defaultRegistry)
           {
               Extinct tmp(*this);
               tmp.loadCurveImpl(extinctName, registryName);
               commitFrom(tmp);
           }
   
           void rebuildCache();
   
           // Observers
   
           [[nodiscard]] auto controls() const -> const io::SimControls& { return controls_; }
   
           [[nodiscard]] auto wlDat() const -> const std::vector<double>& { return wlDat_; }
   
           [[nodiscard]] auto extinctDat() const -> const std::vector<double>& { return extinctDat_; }
   
           [[nodiscard]] auto wl() const -> const std::vector<double>& { return wl_; }
   
           [[nodiscard]] auto extinct() const -> const std::vector<double>& { return extinct_; }
   
           [[nodiscard]] auto wlObs() const -> std::vector<double>;
   
           [[nodiscard]] auto wlOffset() const { return wlOffset_; }
   
           [[nodiscard]] auto applyExtinction(const double A_V, // NOLINT(readability-identifier-naming) -- see above
               const std::vector<double>& spec) const -> std::vector<double>
           {
               if (wl_.empty())
               {
                   throw std::runtime_error(
                       "Extinct::applyExtinction: this Extinct has no wavelength "
                       "grid -- controls_.specsyn() was null the last time "
                       "rebuildCache() ran, so there is nothing to extinguish "
                       "spec against");
               }
               std::vector<double> result(wl_.size());
               for (std::size_t i = 0; i < wl_.size(); i++)
               {
                   result[i] = spec[wlOffset_ + i] * std::exp(-A_V * extinct_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinct_ both have size wl_.size(), and i is bounded by wl_.size(); spec has the same size as controls_.specsyn()->wl() (this class's own contract), so wlOffset_ + i stays in bounds
               }
               return result;
           }
   
           [[nodiscard]] auto applyExtinctionCts(const std::vector<double>& spec,
               const bool nebular = false) const -> std::vector<double>
           {
               if (wl_.empty())
               {
                   throw std::runtime_error(
                       "Extinct::applyExtinctionCts: this Extinct has no "
                       "wavelength grid -- controls_.specsyn() was null the "
                       "last time rebuildCache() ran, so there is nothing to "
                       "extinguish spec against");
               }
               const auto& fac = nebular ? extinctionFacCtsNeb_ : extinctionFacCts_;
               std::vector<double> result(wl_.size());
               for (std::size_t i = 0; i < wl_.size(); i++)
               {
                   result[i] = spec[wlOffset_ + i] * fac[i]; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and fac (extinctionFacCts_ or extinctionFacCtsNeb_) all have size wl_.size() (by computeExtinctionFacCts()'s own contract), and i is bounded by wl_.size(); spec has the same size as the wl originally passed to the constructor, so wlOffset_ + i stays in bounds -- see applyExtinction()'s own identical indexing
               }
               return result;
           }
   
           [[nodiscard]] auto applyExtinctionLines(const double A_V, // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
               const std::vector<double>& lineLum) const -> std::vector<double>
           {
               if (wl_.empty())
               {
                   throw std::runtime_error(
                       "Extinct::applyExtinctionLines: this Extinct has no "
                       "wavelength grid -- controls_.specsyn() was null the "
                       "last time rebuildCache() ran, so extinctLines_ could "
                       "not be meaningfully normalized");
               }
               std::vector<double> result(extinctLines_.size());
               for (std::size_t i = 0; i < extinctLines_.size(); i++)
               {
                   result[i] = lineLum[i] * std::exp(-A_V * extinctLines_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctLines_ both have size extinctLines_.size(), and i is bounded by extinctLines_.size(); lineLum has the same size, by this method's own contract
               }
               return result;
           }
   
           [[nodiscard]] auto applyExtinctionCtsLines(const std::vector<double>& lineLum) const -> std::vector<double>
           {
               if (wl_.empty())
               {
                   throw std::runtime_error(
                       "Extinct::applyExtinctionCtsLines: this Extinct has no "
                       "wavelength grid -- controls_.specsyn() was null the "
                       "last time rebuildCache() ran, so extinctLines_ could "
                       "not be meaningfully normalized");
               }
               std::vector<double> result(extinctLines_.size());
               for (std::size_t i = 0; i < extinctLines_.size(); i++)
               {
                   result[i] = lineLum[i] * extinctionFacCtsLines_[i]; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctionFacCtsLines_ both have size extinctLines_.size() (extinctionFacCtsLines_ by computeExtinctionFacCtsLines()'s own contract), and i is bounded by extinctLines_.size(); lineLum has the same size, by this method's own contract
               }
               return result;
           }
   
       private:
   
           void commitFrom(Extinct& tmp)
           {
               wlDat_ = std::move(tmp.wlDat_);
               extinctDat_ = std::move(tmp.extinctDat_);
               wl_ = std::move(tmp.wl_);
               extinct_ = std::move(tmp.extinct_);
               wlOffset_ = tmp.wlOffset_;
               extinctLines_ = std::move(tmp.extinctLines_);
               extinctionFacCts_ = std::move(tmp.extinctionFacCts_);
               extinctionFacCtsNeb_ = std::move(tmp.extinctionFacCtsNeb_);
               extinctionFacCtsLines_ = std::move(tmp.extinctionFacCtsLines_);
           }
   
           void loadCurveImpl(const std::string& extinctName,
               const std::string& registryName)
           {
               // Locate and parse the registry file
               const auto [registry, registryPath] =
                   utils::parseTOMLFile(registryName, "Extinct");
   
               // Validate that the registry actually lists this curve
               const auto curves = utils::getStringArrayField(registry, "curves");
               if (std::ranges::find(curves, extinctName) == curves.end())
               {
                   throw std::runtime_error(
                       "Extinct: registry " + registryPath.string() +
                       " has no extinction curve '" + extinctName + "'");
               }
   
               // The registry's top-level "file" entry names the HDF5 file
               // holding the actual curve data, relative to the directory
               // containing the registry itself
               const auto h5Name = registry["file"].value<std::string>(); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- toml::table::operator[] is a keyed lookup, not a bounds-checkable container index; a missing key just yields a null node_view, handled by value<std::string>() returning nullopt
               if (!h5Name.has_value())
               {
                   throw std::runtime_error(
                       "Extinct: registry " + registryPath.string() +
                       " is missing required 'file' field");
               }
               const auto h5Path = registryPath.parent_path() / h5Name.value();
   
               // NOLINTBEGIN(misc-include-cleaner) -- see HDF5Utils.hpp's own comment
               const hid_t file = H5Fopen(h5Path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
               if (file < 0)
               {
                   throw std::runtime_error(
                       "Extinct: unable to open HDF5 file " + h5Path.string());
               }
               const hid_t grp = H5Gopen2(file, extinctName.c_str(), H5P_DEFAULT);
               if (grp < 0)
               {
                   H5Fclose(file);
                   throw std::runtime_error(
                       "Extinct: unable to open group " + extinctName +
                       " in HDF5 file " + h5Path.string());
               }
   
               wlDat_ = utils::readDataset1D(grp, "wavelength", "Extinct");
               extinctDat_ = utils::readDataset1D(grp, "kappa", "Extinct");
   
               H5Gclose(grp);
               H5Fclose(file);
               // NOLINTEND(misc-include-cleaner)
   
               rebuildCacheImpl();
           }
   
           static void normalize(const std::vector<double>& wl,
               std::vector<double>& extinct,
               std::vector<double>& extinctLines,
               const std::string& vRegistry = defaultVRegistry)
           {
               const phot::FilterTabulated vFilt("Generic", "Johnson", "V", vRegistry);
   
               // V filter response, converted from wavelength to
               // frequency and reversed back into increasing order (wl
               // increases, so nu = c/wl decreases), then rebuilt as a
               // frequency-space interpolator
               std::vector<double> filtNu(vFilt.wl().size());
               std::ranges::transform(vFilt.wl(), filtNu.begin(),
                   [](const double w) -> double { return utils::c / (w * utils::Angstrom); });
               std::ranges::reverse(filtNu);
               const std::vector<double> filtResp(vFilt.responseData().rbegin(), vFilt.responseData().rend());
               const interp::Interpolator1D<1> respInterp(filtNu, filtResp);
   
               // Same wavelength-to-frequency transformation for the
               // extinction curve
               std::vector<double> extNu(wl.size());
               std::ranges::transform(wl, extNu.begin(),
                   [](const double w) -> double { return utils::c / (w * utils::Angstrom); });
               std::ranges::reverse(extNu);
               const std::vector<double> extRev(extinct.rbegin(), extinct.rend());
               const interp::Interpolator1D<1> extInterp(extNu, extRev);
   
               const double norm = (std::log(100.0) / 5.0) *
                   respInterp.integ(respInterp.xMin(), respInterp.xMax()) /
                   respInterp.integ(extInterp, respInterp.xMin(), respInterp.xMax());
   
               for (double& e : extinct) { e *= norm; }
               for (double& e : extinctLines) { e *= norm; }
           }
   
           [[nodiscard]] auto extinctFac(const double A_V) const -> std::vector<double> // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
           {
               std::vector<double> result(wl_.size());
               for (std::size_t i = 0; i < wl_.size(); i++)
               {
                   result[i] = std::exp(-A_V * extinct_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinct_ both have size wl_.size(), and i is bounded by wl_.size()
               }
               return result;
           }
   
           [[nodiscard]] auto extinctFacLines(const double A_V) const -> std::vector<double> // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
           {
               std::vector<double> result(extinctLines_.size());
               for (std::size_t i = 0; i < extinctLines_.size(); i++)
               {
                   result[i] = std::exp(-A_V * extinctLines_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctLines_ both have size extinctLines_.size(), and i is bounded by extinctLines_.size()
               }
               return result;
           }
   
           void initExtinctLines(const interp::Interpolator1D<1>& interp);
   
           using ExtinctFacFn = std::vector<double> (Extinct::*)(double) const;
   
           [[nodiscard]] auto expectedExtinctFac(ExtinctFacFn fac, std::size_t nInt,
               bool nebular) const -> std::vector<double>;
   
           void computeExtinctionFacCts();
   
           void computeExtinctionFacCtsLines();
   
           void rebuildCacheImpl();
   
           const io::SimControls& controls_; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members) -- deliberately a live reference, not a copy: see this member's own comment for why. Only ever used through the same non-copyable, non-movable ownership pattern (unique_ptr in SimControls's own extinct_) as every other class with a reference member in this codebase (e.g. Specsyn's own controls_), so the usual objection (disabling implicit copy/move assignment) doesn't apply in practice.
   
           std::vector<double> wlDat_;      
           std::vector<double> extinctDat_; 
           std::vector<double> wl_;         
           std::vector<double> extinct_;    
           std::size_t wlOffset_ = 0;       
           std::vector<double> extinctLines_;
   
           std::vector<double> extinctionFacCts_;
   
           std::vector<double> extinctionFacCtsNeb_;
   
           std::vector<double> extinctionFacCtsLines_;
       };
   
   } // namespace extinct
   
   #endif // EXTINCT_HPP