runtime HELMOD_DATA_PATH; add docs

CMakeLists.txt

-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 2.8.12...3.21.2)
 
 if (POLICY CMP0048)
     cmake_policy(SET CMP0048 NEW)
@@ -10,19 +10,20 @@ if(NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE Release)
 endif()
 message(">> CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}")
+message(">> CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
 
 project(libhelmod
     VERSION 2.0.0
     LANGUAGES CXX
-    DESCRIPTION "C++ interface of HelMod")
+    DESCRIPTION "C++ interface for HelMod")
 
-set(HELMOD_DATADIR "" CACHE PATH "with isotope matricies")
-if (HELMOD_DATADIR)
-    set(HELMOD_DATADIR_NOT_BUILTIN 0)
+set(HELMOD_DATA_PATH "" CACHE PATH "with isotope matricies")
+if (HELMOD_DATA_PATH)
+    set(HELMOD_DATA_PATH_NOT_BUILTIN 0)
 else()
-    set(HELMOD_DATADIR_NOT_BUILTIN 1)
+    set(HELMOD_DATA_PATH_NOT_BUILTIN 1)
 endif()
-message(">> HELMOD_DATADIR: ${HELMOD_DATADIR}")
+message(">> HELMOD_DATA_PATH: ${HELMOD_DATA_PATH}")
 
 include(GNUInstallDirs OPTIONAL)
 set(CMAKE_INSTALL_LIBDIR lib CACHE PATH "install libdir")
@@ -30,12 +31,22 @@ set(CMAKE_INSTALL_INCLUDEDIR include CACHE PATH "install includedir")
 
 add_subdirectory(src)
 
-### enable ctest
-include(CTest)
-enable_testing()
-add_subdirectory(test)
+if (BUILD_TESTING)
+    message("-- Test enabled")
+    include(CTest)
+    enable_testing()
+    add_subdirectory(test)
+else()
+    message("-- Test disabled")
+endif()
 
-message(">> CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
+option(BUILD_DOCS "Build and Install Documents (Requires Doxygen)" OFF)
+if (BUILD_DOCS)
+    message("-- Documentation enabled")
+    add_subdirectory(docs)
+else(BUILD_DOCS)
+    message("-- Documentation disabled")
+endif(BUILD_DOCS)
 
 add_custom_target(uninstall
     "${CMAKE_COMMAND}" -P "${CMAKE_SOURCE_DIR}/cmake/uninstall.cmake"

README.md

@@ -36,7 +36,7 @@ Clone this project using git, or download and decompress the source package. Aft
 ```sh
 $ mkdir build
 $ cd build
-$ cmake .. -DHELMOD_DATADIR=/usr/local/share/helmod_data
+$ cmake .. -DHELMOD_DATA_PATH=/usr/local/share/helmod_data
 
 ```
 
@@ -50,7 +50,7 @@ $ make test
 
 ```
 
-Note that the test won't pass if `HELMOD_DATADIR` was not provided during `cmake`. In that case you can still use the library by specifically pass the archive path as function argument at runtime.
+Note that the test won't pass if `HELMOD_DATA_PATH` was not provided during `cmake`. In that case you can still use the library by specifically pass the archive path as function argument at runtime.
 
 Now it's time to install the library into the system.
 
@@ -69,5 +69,5 @@ $ c++ test_mod.cc -o test_mod.exe -lhelmod
 ```
 You can also specify the fallback archive path at compile time.
 ``` sh
-$ c++ test_mod.cc -o test_mod.exe -lhelmod -DHELMOD_DATADIR=/path/to/helmod_data
+$ c++ test_mod.cc -o test_mod.exe -lhelmod -DHELMOD_DATA_PATH=/path/to/helmod_data
 ```

docs/CMakeLists.txt

+find_package(Doxygen REQUIRED)
+
+configure_file(Doxyfile.in Doxyfile)
+
+add_custom_target(docs
+    COMMAND ${DOXYGEN_EXECUTABLE} Doxyfile
+    VERBATIM)

docs/Doxyfile.in

+DOXYFILE_ENCODING      = UTF-8
+PROJECT_NAME           = "@PROJECT_NAME@"
+PROJECT_NUMBER         = "@PROJECT_VERSION@"
+PROJECT_BRIEF          = "@PROJECT_DESCRIPTION@"
+PROJECT_LOGO           =
+OUTPUT_DIRECTORY       =
+
+FULL_PATH_NAMES        = NO
+INPUT                  = "@CMAKE_BINARY_DIR@/src" "@CMAKE_SOURCE_DIR@/src"
+RECURSIVE              = YES
+MACRO_EXPANSION        = YES
+EXTRACT_ALL            = YES
+
+EXCLUDE_SYMBOLS        = __*::* *::*iterator *::*reference
+GENERATE_LATEX         = NO
+
+# USE_MATHJAX            = YES

src/CMakeLists.txt

-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 2.8.12...3.21.2)
 
 file(GLOB HELMOD_SOURCES *.cc)
 file(GLOB HELMOD_HEADERS *.h)

src/HelMod.cc

@@ -3,6 +3,7 @@
 #include "xstream"
 #include <algorithm>
 #include <filesystem>
+#include <fstream>
 #include <limits>
 #include <list>
 #include <map>
@@ -19,197 +20,204 @@ using HelMod::Base;
 using HelMod::Isotope;
 using HelMod::Modulator;
 
+std::string HelMod::data_version(const std::string &data_path) {
+  std::string fn = data_path + "/version.txt";
+  if (std::ifstream ifs;
+      std::filesystem::is_regular_file(fn) && (ifs.open(fn), ifs)) {
+    std::ifstream(fn) >> (fn = "0");
+    return fn;
+  }
+  return "0";
+}
+
 static std::map<std::string, std::map<long, std::string>> _data_ref;
 
-std::string HelMod::IsotopeFile(int Z, int A, const std::string& data_folder) {
-    long za;
-    int* pza = (int*)&za;
-    auto base_path = std::filesystem::weakly_canonical(
-        data_folder.empty() ? _HELMOD_DATADIR_STR : data_folder);
-    auto base_str = base_path.string();
-    if (!std::filesystem::is_directory(base_path)) {
-        throw std::runtime_error("invalid helmod data folder '" + base_str +
-                                 "'");
+std::string HelMod::IsotopeFile(int Z, int A, const std::string &data_path) {
+  long za;
+  int *pza = (int *)&za;
+  std::string base_str;
+  (base_str = data_path).empty() &&
+      (base_str = HelMod::fallback_data_path()).empty() &&
+      (base_str = _HELMOD_DATA_PATH_STR).empty();
+  auto base_path = std::filesystem::weakly_canonical(base_str);
+  base_str = base_path.string();
+  if (!std::filesystem::is_directory(base_path)) {
+    throw std::runtime_error("invalid helmod data path '" + base_str + "', " +
+                             "consider specifying a valid HELMOD_DATA_PATH");
+  }
+  auto entry = _data_ref.find(base_str);
+  if (entry == _data_ref.end()) {
+    // create list of contents
+    entry = _data_ref.insert({base_str, {}}).first;
+    for (const auto &p : std::filesystem::directory_iterator(base_path)) {
+      const std::string filename = p.path().filename().string();
+      // only interested in *.xdat
+      // if (!filename.end_with(".xdat")) // c++20
+      if (filename.size() < 5 ||
+          ".xdat" != filename.substr(filename.size() - 5)) {
+        continue;
+      }
+      std::stringstream ss(filename);
+      char _;
+      if (ss >> pza[0] >> _ >> pza[1]) {
+        entry->second.insert({za, filename});
+      }
     }
-    auto entry = _data_ref.find(base_str);
-    if (entry == _data_ref.end()) {
-	// create list of contents
-	entry = _data_ref.insert({base_str, {}}).first;
-	for (const auto& p: std::filesystem::directory_iterator(base_path)) {
-	    const std::string filename = p.path().filename().string();
-	    // if (!filename.end_with(".xdat")) // c++20
-	    if (filename.size() < 5 || ".xdat" != filename.substr(filename.size() - 5)) {
-		continue;
-	    }
-	    std::stringstream ss(filename);
-	    char _;
-	    if (ss >> pza[0] >> _ >> pza[1]) {
-		entry->second.insert({za, filename});
-	    }
-	}
-    }
-    pza[0] = Z;
-    pza[1] = A;
-    auto rec = entry->second.find(za);
-    return rec == entry->second.end() ? "" : (base_path / rec->second).string();
+  }
+  pza[0] = Z;
+  pza[1] = A;
+  auto rec = entry->second.find(za);
+  return rec == entry->second.end() ? "" : (base_path / rec->second).string();
 }
 
 void Base::_calculate_beta(void) {
-    _beta_inner = std::sqrt(1. - 1. / std::pow(1. + _ekin_inner / T0, 2));
-    _beta_outer = std::sqrt(1. - 1. / std::pow(1. + _ekin_outer / T0, 2));
+  _beta_inner = std::sqrt(1. - 1. / std::pow(1. + _ekin_inner / T0, 2));
+  _beta_outer = std::sqrt(1. - 1. / std::pow(1. + _ekin_outer / T0, 2));
 }
 
 std::valarray<double>
-Base::SparseMatrix::apply(const std::valarray<double>& flux_outer_rig) const {
-    const size_t n_inner = iout_starts.size();
-    std::valarray<double> flux_inner_rig(0., n_inner);
-    for (size_t i = 0; i < n_inner; ++i) {
-        const auto& row     = mat[i];
-        const double istart = iout_starts[i];
-        flux_inner_rig[i] =
-            (row * flux_outer_rig[std::slice(istart, row.size(), 1)]).sum();
-    }
-    return flux_inner_rig;
+Base::SparseMatrix::apply(const std::valarray<double> &flux_outer_rig) const {
+  const size_t n_inner = iout_starts.size();
+  std::valarray<double> flux_inner_rig(0., n_inner);
+  for (size_t i = 0; i < n_inner; ++i) {
+    const auto &row = mat[i];
+    const double istart = iout_starts[i];
+    flux_inner_rig[i] =
+        (row * flux_outer_rig[std::slice(istart, row.size(), 1)]).sum();
+  }
+  return flux_inner_rig;
 }
 
-Isotope::Isotope(const std::string& filename) {
-    yuc::ixstream ifx(filename);
-    if (!ifx) {
-        throw std::runtime_error("unable to read binary data '" + filename +
-                                 "'");
-    }
-
-    std::string isoname2;
-    ifx >> isotope_name >> isoname2;
-    isotope_name += isoname2;
+Isotope::Isotope(const std::string &filename) {
+  yuc::ixstream ifx(filename);
+  if (!ifx) {
+    throw std::runtime_error("unable to read binary data '" + filename + "'");
+  }
 
-    ifx >> *(reinterpret_cast<long*>(&isotope_number));
+  std::string isoname2;
+  ifx >> isotope_name >> isoname2;
+  isotope_name += isoname2;
 
-    ifx >> T0 >> _ekin_outer >> _ekin_inner;
-    _calculate_beta();
+  ifx >> *(reinterpret_cast<long *>(&isotope_number));
 
-    tstart = std::numeric_limits<size_t>::max();
-    tend   = std::numeric_limits<size_t>::min();
-    size_t time;
-    while (ifx >> time) {
-        tstart = std::min(tstart, time);
-        tend   = std::max(tend, time);
+  ifx >> T0 >> _ekin_outer >> _ekin_inner;
+  _calculate_beta();
 
-        _time_slices.push_back({time, {}});
-        auto& mat = _time_slices.back().second;
-        ifx >> mat.iout_starts >> mat.mat;
-    }
-}
+  tstart = std::numeric_limits<size_t>::max();
+  tend = std::numeric_limits<size_t>::min();
+  size_t time;
+  while (ifx >> time) {
+    tstart = std::min(tstart, time);
+    tend = std::max(tend, time);
 
-Modulator Isotope::ExpWindow(size_t time_start, size_t time_end) const {
-    return {*this, time_start, time_end};
+    _time_slices.push_back({time, {}});
+    auto &mat = _time_slices.back().second;
+    ifx >> mat.iout_starts >> mat.mat;
+  }
 }
 
-Modulator::Modulator(const Isotope& iso, size_t time_start, size_t time_end) {
+Modulator::Modulator(const Isotope &iso, size_t time_start, size_t time_end) {
 
-    if (time_start > iso.tend || time_end < iso.tstart) {
-        throw std::runtime_error(
-            "modulation time window out of range. avaliable: [" +          //
-            std::to_string(tstart) + ", " + std::to_string(tend) +         //
-            "],  quested: ["                                               //
-            + std::to_string(time_start) + ", " + std::to_string(time_end) //
-            + "]");
-    }
-
-    isotope_name   = iso.isotope_name;
-    isotope_number = iso.isotope_number;
-    T0             = iso.T0;
+  if (time_start > iso.tend || time_end < iso.tstart) {
+    throw std::runtime_error(
+        "modulation time window out of range. avaliable: [" +          //
+        std::to_string(tstart) + ", " + std::to_string(tend) +         //
+        "],  requested: ["                                             //
+        + std::to_string(time_start) + ", " + std::to_string(time_end) //
+        + "]");
+  }
 
-    _ekin_inner = iso._ekin_inner;
-    _ekin_outer = iso._ekin_outer;
-    _beta_inner = iso._beta_inner;
-    _beta_outer = iso._beta_outer;
+  isotope_name = iso.isotope_name;
+  isotope_number = iso.isotope_number;
+  T0 = iso.T0;
 
-    tstart = std::numeric_limits<size_t>::max();
-    tend   = std::numeric_limits<size_t>::min();
+  _ekin_inner = iso._ekin_inner;
+  _ekin_outer = iso._ekin_outer;
+  _beta_inner = iso._beta_inner;
+  _beta_outer = iso._beta_outer;
 
-    const size_t n_inner = _ekin_inner.size();
-    const size_t n_outer = _ekin_outer.size();
-    const size_t n_slice = iso._time_slices.size();
+  tstart = std::numeric_limits<size_t>::max();
+  tend = std::numeric_limits<size_t>::min();
 
-    _sparse_mat.iout_starts.resize(n_inner, n_outer); // (size, value)
-    _sparse_mat.mat.resize(n_inner);
+  const size_t n_inner = _ekin_inner.size();
+  const size_t n_outer = _ekin_outer.size();
+  const size_t n_slice = iso._time_slices.size();
 
-    auto& istarts = _sparse_mat.iout_starts;
-    std::valarray<size_t> iends((size_t)0, n_inner); // (value, count)
-    std::vector<size_t> time_ids;
+  _sparse_mat.iout_starts.resize(n_inner, n_outer); // (size, value)
+  _sparse_mat.mat.resize(n_inner);
 
-    for (size_t islice = 0; islice < n_slice; ++islice) {
-        const auto& time_slice = iso._time_slices[islice];
-        const auto& time       = time_slice.first;
-        const auto& slice      = time_slice.second;
+  auto &istarts = _sparse_mat.iout_starts;
+  std::valarray<size_t> iends((size_t)0, n_inner); // (value, count)
+  std::vector<size_t> time_ids;
 
-        if (time < time_start || time > time_end) {
-            continue;
-        }
+  for (size_t islice = 0; islice < n_slice; ++islice) {
+    const auto &time_slice = iso._time_slices[islice];
+    const auto &time = time_slice.first;
+    const auto &slice = time_slice.second;
 
-        time_ids.push_back(islice);
+    if (time < time_start || time > time_end) {
+      continue;
+    }
 
-        tstart = std::min(tstart, time);
-        tend   = std::max(tend, time);
+    time_ids.push_back(islice);
 
-        for (size_t j = 0; j < n_inner; ++j) {
-            istarts[j] = std::min(istarts[j], slice.iout_starts[j]);
-            iends[j] =
-                std::max(iends[j], slice.iout_starts[j] + slice.mat[j].size());
-        }
-    }
+    tstart = std::min(tstart, time);
+    tend = std::max(tend, time);
 
     for (size_t j = 0; j < n_inner; ++j) {
-        _sparse_mat.mat[j].resize(iends[j] - istarts[j]);
+      istarts[j] = std::min(istarts[j], slice.iout_starts[j]);
+      iends[j] = std::max(iends[j], slice.iout_starts[j] + slice.mat[j].size());
     }
+  }
 
-    const size_t n_slice_in_window = time_ids.size();
-    for (size_t i = 0; i < n_slice_in_window; ++i) {
-        const auto& slice = iso._time_slices[time_ids[i]].second;
-        std::valarray<size_t> relative_istarts = slice.iout_starts - istarts;
-        for (size_t j = 0; j < n_inner; ++j) {
-            const auto& row = slice.mat[j];
-            _sparse_mat
-                .mat[j][std::slice(relative_istarts[j], row.size(), 1)] += row;
-        }
-    }
+  for (size_t j = 0; j < n_inner; ++j) {
+    _sparse_mat.mat[j].resize(iends[j] - istarts[j]);
+  }
 
+  const size_t n_slice_in_window = time_ids.size();
+  for (size_t i = 0; i < n_slice_in_window; ++i) {
+    const auto &slice = iso._time_slices[time_ids[i]].second;
+    std::valarray<size_t> relative_istarts = slice.iout_starts - istarts;
     for (size_t j = 0; j < n_inner; ++j) {
-        _sparse_mat.mat[j] /= n_slice_in_window;
+      const auto &row = slice.mat[j];
+      _sparse_mat.mat[j][std::slice(relative_istarts[j], row.size(), 1)] += row;
     }
+  }
+
+  for (size_t j = 0; j < n_inner; ++j) {
+    _sparse_mat.mat[j] /= n_slice_in_window;
+  }
 
-    return;
+  return;
 }
 
 std::valarray<double>
-Modulator::modulate(const std::valarray<double>& ekin_lis,
-                    const std::valarray<double>& flux_lis,
-                    const std::valarray<double>& ekin_toa) const {
-    // always assume ekin_lis to envelop ekin_toa and ekin_inner
-
-    // if all energy is high enough, dont no modulation need
-    const double max_mod_e = _ekin_inner.max();
-    if (ekin_lis.min() > max_mod_e) {
-        return ekin_toa.size()
-                   ? yuc::log_interpolator(ekin_lis, flux_lis)(ekin_toa)
-                   : flux_lis;
-    }
-
-    std::valarray<double> used_ekin_toa = ekin_toa.size() ? ekin_toa : ekin_lis;
-
-    yuc::log_interpolator intf_l_e(ekin_lis, flux_lis);
-    // according to HelMod python module
-    std::valarray<double> flux_o_r = intf_l_e(_ekin_outer) / _beta_outer;
-    std::valarray<double> flux_i_e = _sparse_mat.apply(flux_o_r) * _beta_inner;
-    yuc::log_interpolator intf_i_e(_ekin_inner, flux_i_e);
-
-    if (used_ekin_toa.max() <= max_mod_e) {
-        return intf_i_e(used_ekin_toa);
-    } else {
-        for (auto& e : used_ekin_toa) {
-            e = (e > max_mod_e ? intf_l_e : intf_i_e)(e);
-        }
-        return used_ekin_toa;
+Modulator::modulate(const std::valarray<double> &ekin_lis,
+                    const std::valarray<double> &flux_lis,
+                    const std::valarray<double> &ekin_toa) const {
+  // always assume ekin_lis to envelop ekin_toa and ekin_inner
+
+  // if all energy is high enough, dont no modulation need
+  const double max_mod_e = _ekin_inner.max();
+  if (ekin_lis.min() > max_mod_e) {
+    return ekin_toa.size() ? yuc::log_interpolator(ekin_lis, flux_lis)(ekin_toa)
+                           : flux_lis;
+  }
+
+  std::valarray<double> used_ekin_toa = ekin_toa.size() ? ekin_toa : ekin_lis;
+
+  yuc::log_interpolator intf_l_e(ekin_lis, flux_lis);
+  // according to HelMod python module
+  std::valarray<double> flux_o_r = intf_l_e(_ekin_outer) / _beta_outer;
+  std::valarray<double> flux_i_e = _sparse_mat.apply(flux_o_r) * _beta_inner;
+  yuc::log_interpolator intf_i_e(_ekin_inner, flux_i_e);
+
+  if (used_ekin_toa.max() <= max_mod_e) {
+    return intf_i_e(used_ekin_toa);
+  } else {
+    for (auto &e : used_ekin_toa) {
+      e = (e > max_mod_e ? intf_l_e : intf_i_e)(e);
     }
+    return used_ekin_toa;
+  }
 }

src/HelMod.h.in

 #pragma once
 
+#include <cstdlib>
 #include <string>
 #include <valarray>
 #include <vector>
 
-#ifndef HELMOD_DATADIR
-#define _HELMOD_DATADIR_STR "@HELMOD_DATADIR@"
-#if @HELMOD_DATADIR_NOT_BUILTIN@
-#warning "HELMOD_DATADIR not specified and not built into libaray"
+#ifndef HELMOD_DATA_PATH
+#define _HELMOD_DATA_PATH_STR "@HELMOD_DATA_PATH@"
+#if @HELMOD_DATA_PATH_NOT_BUILTIN@ // HELMOD_DATA_PATH_NOT_BUILTIN?
+#warning                                                                       \
+    "HELMOD_DATA_PATH not specified and not built into libaray, runtime environment becomes mandatory"
 #endif
 #else
 #define _helmod_str(x) #x
 #define _helmod_xstr(x) _helmod_str(x)
-#define _HELMOD_DATADIR_STR _helmod_xstr(HELMOD_DATADIR)
+#define _HELMOD_DATA_PATH_STR _helmod_xstr(HELMOD_DATA_PATH)
 #endif
 
 namespace HelMod {
@@ -20,72 +22,142 @@ class Base;
 class Isotope;
 class Modulator;
 
+/** Inline function for fallback data path.
+ *
+ * Returns the first valid `HELMOD_DATA_PATH` in the following order:
+ *
+ * - the environment variable provided when the program is executed
+ * - the macro defined when the caller code was compiled
+ * - the macro defined when the package was configured
+ * - an empty string if none of above is available
+ *
+ * This function is supposed to be called internally.
+ */
+inline std::string fallback_data_path(void) {
+  auto envptr = std::getenv("HELMOD_DATA_PATH");
+  return envptr ? envptr : _HELMOD_DATA_PATH_STR;
+}
+
+/// Read from `version.txt` in `data_path`, return "0" if failed.
+std::string data_version(const std::string &data_path = fallback_data_path());
+
+/** Find the data file storing the modulation matrices for isotope (Z, A).
+ *
+ * \return The absolute path of the matched file, or an empty string if failed
+ * \param Z isotope charge number
+ * \param A isotope mass number, 0 for electron and positron
+ * \param data_path where the data files are stored
+ * - if not empty, explicitly use the given path
+ * - if given empty, use environment, or fallback to the configured path
+ * - if not given, use fallback_data_path
+ */
 std::string IsotopeFile(int Z, int A,
-                        const std::string& data_folder = _HELMOD_DATADIR_STR);
-// return "" if not found
+                        const std::string &data_path = fallback_data_path());
 
+/**
+ * \brief The base class storing infomation and provide data structures for
+ *  solar modulation.
+ */
 class Base {
-  public:
-    std::string isotope_name;
-    struct {
-        int Z, A;
-    } isotope_number;
-    // int Z, A;
-    double T0; // M0/nuc [GeV/n]
-
-    size_t tstart;
-    size_t tend;
-
-  protected:
-    // ekin is actually ekin/nuc [GeV/n]
-    std::valarray<double> _ekin_inner;
-    std::valarray<double> _ekin_outer;
-    std::valarray<double> _beta_inner;
-    std::valarray<double> _beta_outer;
-
-    void _calculate_beta(void);
-
-  protected:
-    struct SparseMatrix {
-        std::valarray<size_t> iout_starts;
-        std::vector<std::valarray<double>> mat;
-        std::valarray<double> apply(const std::valarray<double>&) const;
-    };
-
-    friend class Isotope;
-    friend class Modulator;
-};
-
-class Isotope : public Base {
-  public:
-    Isotope(const std::string& filename);
-    Isotope(int Z, int A, const std::string& data_folder = _HELMOD_DATADIR_STR)
-        : Isotope(IsotopeFile(Z, A, data_folder)) {}
-    Modulator ExpWindow(size_t time_start, size_t time_end) const;
-
-  protected:
-    std::vector<std::pair<size_t, SparseMatrix>> _time_slices;
-
-    friend class Modulator;
+public:
+  std::string isotope_name;
+  struct {
+    /// charge number
+    int Z;
+    /// mass number, 0 for electron and positron
+    int A;
+  } isotope_number;
+  /// rest mass per nucleon [GeV/n]
+  double T0;
+
+  /// start date in `YYYYMMDD`
+  size_t tstart;
+  // end date in `YYYYMMDD`
+  size_t tend;
+
+protected:
+  std::valarray<double> _ekin_inner; // [GeV/n]
+  std::valarray<double> _ekin_outer; // [GeV/n]
+  std::valarray<double> _beta_inner;
+  std::valarray<double> _beta_outer;
+
+  void _calculate_beta(void);
+
+protected:
+  struct SparseMatrix {
+    std::valarray<size_t> iout_starts;
+    std::vector<std::valarray<double>> mat;
+    std::valarray<double> apply(const std::valarray<double> &) const;
+  };
+
+  friend class Isotope;
+  friend class Modulator;
 };
 
+/// Apply modulation matrices to cosmic-ray spectra
 class Modulator : public Base {
-  public:
-    Modulator(const Modulator&) = default;
-    Modulator(const Isotope&, //
-                 size_t time_start, size_t time_end);
-
-    std::valarray<double>
-    modulate(const std::valarray<double>& ekin_lis,
-             const std::valarray<double>& flux_lis,
-             const std::valarray<double>& ekin_toa = {}) const;
-
-  protected:
-    SparseMatrix _sparse_mat;
-
-  protected:
-    Modulator() = default;
+public:
+  /** Construct for an isotope between a time window.
+   * It is convenient to create Isotope objects using *initializer_lists*,
+   * e.g., modulator for proton flux between AMS-02 time window:
+   * ``` C++
+   * HelMod::Modulator mod_ams02_proton({1, 1}, 20110519, 20131126);
+   * ```
+   */
+  Modulator(const Isotope &, size_t time_start, size_t time_end);
+  Modulator(const Modulator &) = default;
+
+  /** Modulate a local interstellar cosmic-ray spectum.
+   *
+   * \param ekin_lis [GeV/n]
+   *        The kinetic energy per nucleon of input LIS.
+   *
+   * \param flux_lis [~1/m^2/sr/s/(GeV/n)]
+   *        The kinetic energy per nucleon of input LIS.  Other units are valid
+   *        as long as they are linearly proportional to *per kinetic energy*.
+   *
+   * \param ekin_toa [GeV/n]
+   *        The kinetic energy per nucleon of output modulated flux,
+   *        if empty, align with input `ekin_lis`
+   *
+   * \return The modulated flux aligned with `ekin_toa`,
+   *         in the same unit as `flux_lis`.
+   */
+  std::valarray<double>
+  modulate(const std::valarray<double> &ekin_lis,
+           const std::valarray<double> &flux_lis,
+           const std::valarray<double> &ekin_toa = {}) const;
+
+protected:
+  SparseMatrix _sparse_mat;
+
+protected:
+  Modulator() = default;
+
+  friend class Isotope;
+};
 
-    friend class Isotope;
+/// Access modulation matrices for a single isotope.
+class Isotope : public Base {
+public:
+  /// Construct with data file
+  Isotope(const std::string &filename);
+
+  /// Construct with isotope numbers (Z, A)
+  Isotope(int Z, int A, const std::string &data_folder = fallback_data_path())
+      : Isotope(IsotopeFile(Z, A, data_folder)) {}
+
+  /** Construct a Modulator object for this isotope in a time window.
+   * \param time_start start date in `YYYYMMDD`
+   * \param time_end     end date in `YYYYMMDD`
+   */
+  Modulator exp_window(size_t time_start, size_t time_end) const {
+    return {*this, time_start, time_end};
+  }
+
+protected:
+  std::vector<std::pair<size_t, SparseMatrix>> _time_slices;
+
+  friend class Modulator;
 };
 }; // namespace HelMod

test/CMakeLists.txt

-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 2.8.12...3.21.2)
 
 file(GLOB Test_SOURCES test_mod.cc)
 file(GLOB Test_DATA proton_lis.tsv)