initial commit

.gitignore

+build
+deprecated
+*.disabled
+
+# config
+.ycm_extra_conf.py
+.cache
+
+# Editor #
+.*.sw?
+*~
+
+# Compiled source #
+###################
+*.com
+*.class
+*.dll
+*.exe
+*.o
+*.so
+
+# Packages #
+############
+# it's better to unpack these files and commit the raw source
+# git has its own built in compression methods
+*.7z
+*.dmg
+*.gz
+*.iso
+*.jar
+*.rar
+*.tar
+*.zip
+
+# Logs and databases #
+######################
+*.log
+*.sql
+*.sqlite
+
+# OS generated files #
+######################
+.DS_Store
+.DS_Store?
+._*
+.Spotlight-V100
+.Trashes
+ehthumbs.db
+Thumbs.db

CMakeLists.txt

+cmake_minimum_required(VERSION 2.8)
+
+if (POLICY CMP0048)
+    cmake_policy(SET CMP0048 NEW)
+endif (POLICY CMP0048)
+
+set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall -Wextra")
+set(CMAKE_CXX_FLAGS_RELEASE "-O2")
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+endif()
+message(">> CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}")
+
+project(libhelmod
+    VERSION 2.0.0
+    LANGUAGES CXX
+    DESCRIPTION "C++ interface of HelMod")
+
+set(HELMOD_DATADIR "" CACHE PATH "with isotope matricies")
+if (HELMOD_DATADIR)
+    set(HELMOD_DATADIR_NOT_BUILTIN 0)
+else()
+    set(HELMOD_DATADIR_NOT_BUILTIN 1)
+endif()
+message(">> HELMOD_DATADIR: ${HELMOD_DATADIR}")
+
+include(GNUInstallDirs OPTIONAL)
+set(CMAKE_INSTALL_LIBDIR lib CACHE PATH "install libdir")
+set(CMAKE_INSTALL_INCLUDEDIR include CACHE PATH "install includedir")
+
+add_subdirectory(src)
+
+### enable ctest
+include(CTest)
+enable_testing()
+add_subdirectory(test)
+
+message(">> CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
+
+add_custom_target(uninstall
+    "${CMAKE_COMMAND}" -P "${CMAKE_SOURCE_DIR}/cmake/uninstall.cmake"
+)

README.md

+C++ interface for the [HelMod](http://helmod.org) offline calculator.
+
+See upper stream pages [History and Citation](http://www.helmod.org/index.php/history-and-citation) and [Bibliography](http://www.helmod.org/index.php/publications).
+
+## Instruction
+
+### 1. Prepare modulator archive
+
+In order to use this C++ library to calculate solar modulation of cosmic rays (CR), modulator archive must be provided. The archive stores time dependent matrices that are used to transform a local interstellar CR spectrum (LIS) into modulated spectrum at the top of atmosphere (TOA).
+
+This library reads a reorganized version of such archive which is more disk-friendly than the official one. You can [download the latest archive here](https://cloud.itp.ac.cn/d/2927e56f3364406c87b5/).
+
+Let's say the tarball `helmod_data-4.1.2.tar.gz` has been download and is stored in the current working directory. As an example, to extract the archive file into `/usr/local/share/`, execute
+
+```sh
+$ sudo mkdir -pv /usr/local/share
+$ sudo tar -xzvf ./helmod_data-4.1.2.tar.gz -C /usr/local/share
+
+```
+
+A directory containing all the modulation matrices will be created as `/usr/local/share/helmod_data-4.1.2/`.  It's recommended to create a symbolic link as an entry point that later provided to the library.
+
+```sh
+$ sudo ln -sv helmod_data-4.1.2 /usr/local/share/helmod_data
+$ # check the link
+$ ls -ld /usr/local/share/helmod_data*
+
+```
+
+### 2. Compile the library
+
+This library is built using `cmake`. Make sure it's presented on your computer.
+
+Clone this project using git, or download and decompress the source package. After changing your working directory into the project root directory (containing the README you currently reading), create a new directory to build inside.
+
+```sh
+$ mkdir build
+$ cd build
+$ cmake .. -DHELMOD_DATADIR=/usr/local/share/helmod_data
+
+```
+
+Note that without specifying `CMAKE_INSTALL_PREFIX`, the library will later be installed to `/usr/local`.
+
+If the cmake command terminated successfully, you can now build and test the library.
+
+```sh
+$ make
+$ 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.
+
+Now it's time to install the library into the system.
+
+```sh
+$ sudo make install
+
+```
+
+### 3. Use in C++ program
+
+See `test/test_mod.cc` for a demonstration of modulate proton flux during AMS02 mission window.
+
+To compile this code, run
+``` sh
+$ 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
+```

cmake/uninstall.cmake

+set(MANIFEST "${CMAKE_CURRENT_BINARY_DIR}/install_manifest.txt")
+
+if(NOT EXISTS ${MANIFEST})
+    message(FATAL_ERROR "Cannot find install manifest: '${MANIFEST}'")
+endif()
+
+file(STRINGS ${MANIFEST} files)
+foreach(file ${files})
+    file(GLOB parent_content ${file}/* ${file}/.*)
+    list(LENGTH parent_content parent_len)
+    while(parent_len EQUAL 0)
+	message(STATUS "Removing: '${file}'")
+	file(REMOVE_RECURSE ${file})
+	if(EXISTS ${file})
+	    MESSAGE(STATUS "'${file}' not removed")
+	endif()
+	get_filename_component(file ${file} DIRECTORY)
+	file(GLOB parent_content ${file}/* ${file}/.*)
+	list(LENGTH parent_content parent_len)
+    endwhile()
+
+endforeach(file)

src/CMakeLists.txt

+cmake_minimum_required(VERSION 2.8)
+
+file(GLOB HELMOD_SOURCES *.cc)
+file(GLOB HELMOD_HEADERS *.h)
+
+configure_file(HelMod.h.in HelMod.h @ONLY)
+
+add_library(helmod SHARED ${HELMOD_SOURCES})
+target_include_directories(helmod PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+				  PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
+
+set_target_properties(helmod PROPERTIES
+    LINKER_LANGUAGE CXX
+    CXX_STANDARD 17
+    PUBLIC_HEADER "${HELMOD_HEADERS}"
+    )
+
+install(TARGETS helmod)
+
+install(FILES ${CMAKE_CURRENT_BINARY_DIR}/HelMod.h
+    DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

src/HelMod.cc

+#include "HelMod.h"
+#include "loginterp"
+#include "xstream"
+#include <algorithm>
+#include <filesystem>
+#include <limits>
+#include <list>
+#include <map>
+#include <sstream>
+#include <stdexcept>
+#include <valarray>
+
+// labeling convention:
+// i(nner), o(uter): helmode matrix grid
+// t(oa)  , l(is)  : user interface grid
+// e(kn)  , r(ig)  : energy quantity
+
+using HelMod::Base;
+using HelMod::Isotope;
+using HelMod::Modulator;
+
+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 +
+                                 "'");
+    }
+    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();
+}
+
+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));
+}
+
+std::valarray<double>
+Base::SparseMatrix::apply(const std::valarray<double>& flux_outer_rig) const {
+    const size_t n_inner = iout_starts.size();
+    // if (flux_outer_rig.size() != n_inner) {
+    //     throw std::logic_error("wrong size of `flux_outer', "
+    //                            "should interpolate onto `_ekin_outer' first");
+    // }
+    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;
+
+    ifx >> *(reinterpret_cast<long*>(&isotope_number));
+
+    ifx >> T0 >> _ekin_outer >> _ekin_inner;
+    _calculate_beta();
+
+    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);
+
+        _time_slices.push_back({time, {}});
+        auto& mat = _time_slices.back().second;
+        ifx >> mat.iout_starts >> mat.mat;
+    }
+}
+
+Modulator Isotope::ExpWindow(size_t time_start, size_t time_end) const {
+    return {*this, time_start, 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;
+
+    _ekin_inner = iso._ekin_inner;
+    _ekin_outer = iso._ekin_outer;
+    _beta_inner = iso._beta_inner;
+    _beta_outer = iso._beta_outer;
+
+    tstart = std::numeric_limits<size_t>::max();
+    tend   = std::numeric_limits<size_t>::min();
+
+    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();
+
+    _sparse_mat.iout_starts.resize(n_inner, n_outer); // (size, value)
+    _sparse_mat.mat.resize(n_inner);
+
+    auto& istarts = _sparse_mat.iout_starts;
+    std::valarray<size_t> iends((size_t)0, n_inner); // (value, count)
+    std::vector<size_t> time_ids;
+
+    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;
+
+        if (time < time_start || time > time_end) {
+            continue;
+        }
+
+        time_ids.push_back(islice);
+
+        tstart = std::min(tstart, time);
+        tend   = std::max(tend, time);
+
+        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());
+        }
+    }
+
+    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) {
+            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;
+}
+
+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;
+    }
+}

src/HelMod.h.in

+#pragma once
+
+#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"
+#endif
+#else
+#define _helmod_str(x) #x
+#define _helmod_xstr(x) _helmod_str(x)
+#define _HELMOD_DATADIR_STR _helmod_xstr(HELMOD_DATADIR)
+#endif
+
+namespace HelMod {
+class Base;
+class Isotope;
+class Modulator;
+
+std::string IsotopeFile(int Z, int A,
+                        const std::string& data_folder = _HELMOD_DATADIR_STR);
+// return "" if not found
+
+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;
+};
+
+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;
+
+    friend class Isotope;
+};
+}; // namespace HelMod

src/loginterp

+// version 2.0
+#pragma once
+
+#include <cmath>
+#include <set>
+
+namespace yuc {
+class log_interpolator {
+  protected:
+    struct data_point {
+        double x, y;
+        bool operator<(const data_point& d) const { return x < d.x; }
+    };
+
+  protected:
+    std::set<data_point> data;
+
+  public:
+    log_interpolator(void) {}
+    template <typename V> log_interpolator(const V& xa, const V& ya) {
+        insert(xa, ya);
+    }
+
+  public:
+    void clear(void) { data.clear(); }
+    void insert(double x, double y) {
+        if (x >= 0 && y >= 0) {
+            data.insert({std::log(x), std::log(y)});
+        }
+    }
+    template <typename V> void insert(const V& xa, const V& ya) {
+        for (size_t i = 0; i < xa.size(); ++i) {
+            this->insert(xa[i], ya[i]);
+        }
+    }
+
+  public:
+    double operator()(double x) const {
+        if (x < 0 || data.size() < 2) {
+            return NAN;
+        }
+        x = std::log(x);
+
+        auto it1 = data.upper_bound({x}), it2 = std::prev(it1);
+        if (it1 == data.begin()) {
+            it2 = std::next(it1);
+        } else if (it1 == data.end()) {
+            it1 = std::prev(it2);
+        }
+
+        const double x1 = it1->x, y1 = it1->y;
+        const double x2 = it2->x, y2 = it2->y;
+        if (x1 == x2) {
+            return std::exp(0.5 * (y1 + y2));
+        }
+        if (y1 == -HUGE_VAL || y2 == -HUGE_VAL) {
+            if (y1 == y2) {
+                return 0.;
+            } else {
+                // use log-linear interpolation if one y is zero
+                return std::exp(y1) +
+                       (x - x1) * (std::exp(y2) - std::exp(y1)) / (x2 - x1);
+            }
+        }
+        return std::exp(y1 + (x - x1) * (y2 - y1) / (x2 - x1));
+        // (power-law extrapolating behavior)
+    }
+    template <typename V> V operator()(const V& xs) const {
+        V ys(xs.size());
+        for (size_t i = 0; i < xs.size(); ++i) {
+            ys[i] = this->operator()(xs[i]);
+        }
+        return ys;
+    }
+};
+};
+
+// vi:ft=cpp

src/xstream

+#pragma once
+#include <fstream>
+#include <string>
+#include <vector>
+#include <valarray>
+
+// Data is always stored little_endian
+#if BYTE_ORDER == BIG_ENDIAN
+#define OXSTREAM_NUMERIC_WRITE write_reverse
+#define IXSTREAM_NUMERIC_READ read_reverse
+#else
+#define OXSTREAM_NUMERIC_WRITE write
+#define IXSTREAM_NUMERIC_READ read
+#endif
+
+namespace yuc {
+class oxstream : protected std::ofstream {
+  protected:
+    static constexpr size_t cell_witdh = 8;
+
+  public:
+    using std::ofstream::ofstream;
+