code cleanup

Author: YihanWangAstro

CHANGELOG.md

@@ -158,7 +158,7 @@ data.add_flux(nu_min=1e14, nu_max=1e15, num_points=5, t=times, flux=flux, err=er
   - **ConfigParams**: `fwd_SSC` → `fwd_ssc`, `rvs_SSC` → `rvs_ssc`, `IC_cooling` → `ssc_cooling`, `KN` → `kn`
   - **PyRadiation**: `IC_cooling` → `ssc_cooling`, `SSC` → `ssc`, `KN` → `kn`
   - **Model Constructor**: `forward_rad` → `fwd_rad`, `reverse_rad` → `rvs_rad`
-  - All parameter names now follow consistent snake_case convention
+  - All parameter names now follow the consistent snake_case convention
 
 #### **Enhanced Code Documentation**
 - **Comprehensive Comments**: Added detailed documentation to all major classes and structures

include/forward-shock.hpp

@@ -19,7 +19,7 @@
  *          based on if the Ejecta class supports mass and energy injection methods.
  * <!-- ************************************************************************************** -->
  */
-template <typename Ejecta, typename Medium>
+template <typename Ejecta>
 struct ForwardState {
     static constexpr bool mass_inject = HasDmdt<Ejecta>;   ///< whether Ejecta class has dmdt method
     static constexpr bool energy_inject = HasDedt<Ejecta>; ///< whether Ejecta class has dedt method
@@ -62,7 +62,7 @@ struct ForwardState {
 template <typename Ejecta, typename Medium>
 class ForwardShockEqn {
   public:
-    using State = ForwardState<Ejecta, Medium>;
+    using State = ForwardState<Ejecta>;
 
     /**
      * <!-- ************************************************************************************** -->

include/inverse-compton.h

@@ -148,7 +148,7 @@ struct ICPhoton {
 
     /**
      * <!-- ************************************************************************************** -->
-     * @brief Computes the base-2 logarithm of the photon specific intensity at a given frequency.
+     * @brief Computes the base-2 logarithm of the photon-specific intensity at a given frequency.
      * @param log2_nu The base-2 logarithm of the frequency
      * @return The base-2 logarithm of the photon specific intensity at the given frequency
      * <!-- ************************************************************************************** -->
@@ -162,20 +162,16 @@ struct ICPhoton {
   private:
     void generate_grid();
 
-    void fill_integral_table();
-
     Array nu0; // input frequency nu0 grid value
 
     Array I_nu_dnu; // input I_nu
 
     Array gamma; // gamma grid boundary values
 
-    Array N_e_dgamma; // electron column density
+    Array dN_e; // electron column density
 
     MeshGrid IC_tab;
 
-    Real min_nu_quested{con::inf};
-
     static constexpr size_t gamma_grid_per_order{7}; // Number of frequency bins
 
     static constexpr size_t nu_grid_per_order{5}; // Number of gamma bins
@@ -203,13 +199,13 @@ inline constexpr Real IC_x0 = 0.47140452079103166;
  * @tparam Photons Type of the photon distribution
  * @param electrons The electron grid
  * @param photons The photon grid
- * @param kn flag for Klein-Nishina
+ * @param KN flag for Klein-Nishina
  * @return A 3D grid of IC photons
  * <!-- ************************************************************************************** -->
  */
 template <typename Electrons, typename Photons>
 ICPhotonGrid<Electrons, Photons> generate_IC_photons(ElectronGrid<Electrons> const& electrons,
-                                                     PhotonGrid<Photons> const& photons, bool kn = true) noexcept;
+                                                     PhotonGrid<Photons> const& photons, bool KN = true) noexcept;
 
 /**
  * <!-- ************************************************************************************** -->
@@ -245,9 +241,6 @@ template <typename Electrons, typename Photons>
 ICPhoton<Electrons, Photons>::ICPhoton(Electrons const& electrons, Photons const& photons, bool KN) noexcept
     : photons(photons), electrons(electrons), KN(KN) {}
 
-template <typename Electrons, typename Photons>
-void ICPhoton<Electrons, Photons>::fill_integral_table() {}
-
 template <typename Electrons, typename Photons>
 void ICPhoton<Electrons, Photons>::generate_grid() {
     const Real gamma_min = std::min(electrons.gamma_m, electrons.gamma_c);
@@ -258,24 +251,20 @@ void ICPhoton<Electrons, Photons>::generate_grid() {
     const Real nu_max = photons.nu_M * 10;
     const size_t nu_size = static_cast<size_t>(std::log10(nu_max / nu_min) * nu_grid_per_order);
 
-    Array dnu0;
-
-    Array dgamma;
+    Array dnu0, dgamma;
 
     logspace_boundary_center(std::log2(nu_min), std::log2(nu_max), nu_size, nu0, dnu0);
 
     logspace_boundary_center(std::log2(gamma_min), std::log2(gamma_max), gamma_size, gamma, dgamma);
 
-    // I_nu = Array({nu_size}, -1);
-    // column_den = Array({gamma_size}, -1);
     IC_tab = MeshGrid({gamma_size, nu_size}, -1);
     generated = true;
 
     I_nu_dnu = Array::from_shape({nu_size});
-    N_e_dgamma = Array::from_shape({gamma_size});
+    dN_e = Array::from_shape({gamma_size});
 
     for (size_t i = 0; i < gamma_size; ++i) {
-        N_e_dgamma(i) = electrons.compute_column_den(gamma(i)) * dgamma(i);
+        dN_e(i) = electrons.compute_column_den(gamma(i)) * dgamma(i);
     }
 
     for (size_t j = 0; j < nu_size; ++j) {
@@ -299,7 +288,7 @@ Real ICPhoton<Electrons, Photons>::compute_I_nu(Real nu) {
     for (size_t i = gamma_size; i-- > 0;) {
         const Real gamma_i = gamma(i);
         const Real upscatter = 4 * IC_x0 * gamma_i * gamma_i;
-        const Real Ndgamma = N_e_dgamma(i);
+        const Real Ndgamma = dN_e(i);
 
         if (nu > upscatter * nu0.back())
             break;
@@ -324,6 +313,7 @@ Real ICPhoton<Electrons, Photons>::compute_I_nu(Real nu) {
                 break;
             }
         }
+
         if (extrapolate) {
             IC_I_nu += Ndgamma *
                        (IC_tab(i, 0) + (IC_tab(i, 0) - IC_tab(i, 1)) / (nu0(1) - nu0(0)) * (nu0(0) - nu / upscatter));

include/mesh.h

@@ -164,6 +164,8 @@ Coord auto_grid(Ejecta const& jet, Array const& t_obs, Real theta_cut, Real thet
  * @tparam Ejecta Type of the jet/ejecta class
  * @param jet The jet/ejecta object
  * @param gamma_cut Lorentz factor cutoff value
+ * @param phi_resol
+ * @param theta_resol
  * @param is_axisymmetric Flag for axisymmetric jets
  * @return Angle (in radians) at which the jet's Lorentz factor drops to gamma_cut
  * <!-- ************************************************************************************** -->

include/observer.h

@@ -333,10 +333,7 @@ MeshGrid Observer::specific_flux(Array const& t_obs, Array const& nu_obs, Photon
 
             InterpState state;
             for (size_t k = 0; k < t_grid - 1 && t_idx < t_obs_len; k++) {
-                const Real t_hi = time(i, j, k + 1);
-                if (t_hi < t_obs(t_idx)) {
-                    continue;
-                } else {
+                if (const Real t_hi = time(i, j, k + 1); t_hi >= t_obs(t_idx)) {
                     const size_t idx_start = t_idx;
                     iterate_to(t_hi, t_obs, t_idx);
                     const size_t idx_end = t_idx;
@@ -395,7 +392,6 @@ Array Observer::specific_flux_series(Array const& t_obs, Array const& nu_obs, Ph
             while (t_idx < t_obs_len && k < t_grid - 1) {
                 if (time(i, j, k + 1) < t_obs(t_idx)) {
                     k++;
-                    continue;
                 } else {
                     if (set_boundaries(state, eff_i, i, j, k, lg2_nu_src(t_idx), photons)) [[likely]] {
                         //F_nu(t_idx) += loglog_interpolate(state, lg2_t_obs(t_idx), lg2_t(i, j, k));

include/reverse-shock.hpp

@@ -18,7 +18,7 @@
  *          injection capabilities.
  * <!-- ************************************************************************************** -->
  */
-template <typename Ejecta, typename Medium>
+template <typename Ejecta>
 struct ReverseState {
     static constexpr bool mass_inject = HasDmdt<Ejecta>;   ///< Whether ejecta has mass injection
     static constexpr bool energy_inject = HasDedt<Ejecta>; ///< Whether ejecta has energy injection
@@ -56,21 +56,21 @@ struct ReverseState {
 template <typename Ejecta, typename Medium>
 class FRShockEqn {
   public:
-    using State = ReverseState<Ejecta, Medium>;
+    using State = ReverseState<Ejecta>;
 
     /**
      * <!-- ************************************************************************************** -->
      * @brief Constructor for the FRShockEqn class.
      * @details Initializes the forward-reverse shock equation with the given medium, ejecta, and parameters.
      * @param medium The medium through which the shock propagates
-     * @param jet The ejecta driving the shock
+     * @param ejecta The ejecta driving the shock
      * @param phi Azimuthal angle
      * @param theta Polar angle
      * @param rad_fwd Radiation params for forward shock
      * @param rad_rvs Radiation params for reverse shock
      * <!-- ************************************************************************************** -->
      */
-    FRShockEqn(Medium const& medium, Ejecta const& jet, Real phi, Real theta, RadParams const& rad_fwd,
+    FRShockEqn(Medium const& medium, Ejecta const& ejecta, Real phi, Real theta, RadParams const& rad_fwd,
                RadParams const& rad_rvs);
 
     Medium const& medium;    ///< Reference to the medium properties

include/utilities.h

@@ -346,12 +346,12 @@ struct Tables {
     std::array<double, BUCKETS> log2_c{};
     std::array<uint64_t, BUCKETS> c_bits{}; // optional (not used in hot path)
     Tables() {
-        const int shift = 52 - BUCKET_BITS;
-        const uint64_t expbits = uint64_t(EXP_BIAS) << 52;
+        constexpr int shift = 52 - BUCKET_BITS;
+        constexpr uint64_t expbits = static_cast<uint64_t>(EXP_BIAS) << 52;
         for (size_t i = 0; i < BUCKETS; ++i) {
-            uint64_t frac_mid = (uint64_t(i) << shift) | (uint64_t(1) << (shift - 1));
+            const uint64_t frac_mid = (static_cast<uint64_t>(i) << shift) | (static_cast<uint64_t>(1) << (shift - 1));
             uint64_t bits = expbits | frac_mid; // center c_i bits
-            double c = std::bit_cast<double>(bits);
+            const double c = std::bit_cast<double>(bits);
             c_bits[i] = bits;
             inv_c[i] = 1.0 / c;
             log2_c[i] = std::log2(c);
@@ -366,17 +366,17 @@ inline const Tables& tables() {
 
 inline double decompose_normals(double x, int& e) {
     uint64_t bits = std::bit_cast<uint64_t>(x);
-    uint64_t expo = (bits >> 52) & EXP_MASK;
+    const uint64_t expo = (bits >> 52) & EXP_MASK;
     if (expo == 0 || expo == EXP_MASK) {
         int ee;
         double m = std::frexp(x, &ee); // [0.5,1)
         m *= 2.0;
         e = ee - 1;
         return m; // -> [1,2)
     }
-    e = int(expo) - EXP_BIAS;
-    bits = (bits & FRAC_MASK) | (uint64_t(EXP_BIAS) << 52); // force mantissa exponent
-    return std::bit_cast<double>(bits);                     // m in [1,2)
+    e = static_cast<int>(expo) - EXP_BIAS;
+    bits = (bits & FRAC_MASK) | (static_cast<uint64_t>(EXP_BIAS) << 52); // force mantissa exponent
+    return std::bit_cast<double>(bits);                                  // m in [1,2)
 }
 
 /**

pybind/error_handling.h

@@ -12,13 +12,13 @@
 
 namespace afterglow {
 
-    class ValidationError : public std::invalid_argument {
+    class ValidationError final : public std::invalid_argument {
       public:
         explicit ValidationError(const std::string& message) : std::invalid_argument(message) {}
         explicit ValidationError(const char* message) : std::invalid_argument(message) {}
     };
 
-    class LogicError : public std::logic_error {
+    class LogicError final : public std::logic_error {
       public:
         explicit LogicError(const std::string& message) : std::logic_error(message) {}
         explicit LogicError(const char* message) : std::logic_error(message) {}

pybind/mcmc.h

@@ -268,7 +268,7 @@ struct MultiBandModel {
      * @param data Complete multi-wavelength observational dataset for fitting
      * <!-- ************************************************************************************** -->
      */
-    MultiBandModel(MultiBandData data);
+    explicit MultiBandModel(MultiBandData data);
 
     /**
      * <!-- ************************************************************************************** -->
@@ -334,15 +334,15 @@ struct MultiBandModel {
      * @param t_min Minimum time for photon grid generation [seconds]
      * @param t_max Maximum time for photon grid generation [seconds]
      * @param obs Observer object for flux calculations
-     * @param f_photons Forward shock synchrotron photon grid
-     * @param r_photons Reverse shock synchrotron photon grid
-     * @param f_IC_photons Forward shock inverse Compton photon grid
-     * @param r_IC_photons Reverse shock inverse Compton photon grid
+     * @param fwd_photons Forward shock synchrotron photon grid
+     * @param rvs_photons Reverse shock synchrotron photon grid
+     * @param fwd_IC_photons Forward shock inverse Compton photon grid
+     * @param rvs_IC_photons Reverse shock inverse Compton photon grid
      * <!-- ************************************************************************************** -->
      */
     template <typename ICPhotonGrid>
-    void generate_photons(Params const& param, double t_min, double t_max, Observer& obs, SynPhotonGrid& f_photons,
-                          SynPhotonGrid& r_photons, ICPhotonGrid& f_IC_photons, ICPhotonGrid& r_IC_photons);
+    void generate_photons(Params const& param, double t_min, double t_max, Observer& obs, SynPhotonGrid& fwd_photons,
+                          SynPhotonGrid& rvs_photons, ICPhotonGrid& fwd_IC_photons, ICPhotonGrid& rvs_IC_photons);
 
     /**
      * <!-- ************************************************************************************** -->

pybind/pymodel.cpp

@@ -199,8 +199,8 @@ void save_photon_details(PhotonGrid const& photons, PyShock& details) {
     }
 }
 
-void PyModel::single_evo_details(Shock const& shock, Coord const& coord, Array const& t_obs, Observer& obs,
-                                 PyRadiation const& rad, PyShock& details) const {
+void PyModel::single_evo_details(Shock const& shock, Coord const& coord, Observer& obs, PyRadiation const& rad,
+                                 PyShock& details) const {
     obs.observe(coord, shock, obs_setup.lumi_dist, obs_setup.z);
 
     details.t_obs = obs.time / unit::sec;
@@ -239,7 +239,7 @@ auto PyModel::details(Real t_min, Real t_max) const -> PyDetails {
 
         save_shock_details(fwd_shock, details.fwd);
 
-        single_evo_details(fwd_shock, coord, t_obs, observer, fwd_rad, details.fwd);
+        single_evo_details(fwd_shock, coord, observer, fwd_rad, details.fwd);
 
         return details;
     } else {
@@ -250,9 +250,9 @@ auto PyModel::details(Real t_min, Real t_max) const -> PyDetails {
 
         save_shock_details(rvs_shock, details.rvs);
 
-        single_evo_details(fwd_shock, coord, t_obs, observer, fwd_rad, details.fwd);
+        single_evo_details(fwd_shock, coord, observer, fwd_rad, details.fwd);
 
-        single_evo_details(rvs_shock, coord, t_obs, observer, rvs_rad, details.rvs);
+        single_evo_details(rvs_shock, coord, observer, rvs_rad, details.rvs);
 
         return details;
     }
@@ -284,8 +284,8 @@ auto PyModel::flux_density(PyArray const& t, PyArray const& nu) -> PyFlux {
     const Array t_obs = t * unit::sec;
     const Array nu_obs = nu * unit::Hz;
 
-    auto flux_func = [](Observer& obs, Array const& t, Array const& nu, auto& photons) -> XTArray {
-        return obs.specific_flux_series(t, nu, photons) / unit::flux_den_cgs;
+    auto flux_func = [](Observer& obs, Array const& time, Array const& freq, auto& photons) -> XTArray {
+        return obs.specific_flux_series(time, freq, photons) / unit::flux_den_cgs;
     };
 
     auto result = compute_emission(t_obs, nu_obs, flux_func);
@@ -300,8 +300,8 @@ auto PyModel::flux(PyArray const& t, double nu_min, double nu_max, size_t num_nu
     const Array nu_obs = xt::logspace(std::log10(nu_min * unit::Hz), std::log10(nu_max * unit::Hz), num_nu);
     const Array t_obs = t * unit::sec;
 
-    auto flux_func = [](Observer& obs, Array const& t, Array const& nu, auto& photons) -> XTArray {
-        return obs.flux(t, nu, photons) / unit::flux_cgs;
+    auto flux_func = [](Observer& obs, Array const& time, Array const& freq, auto& photons) -> XTArray {
+        return obs.flux(time, freq, photons) / unit::flux_cgs;
     };
 
     auto result = compute_emission(t_obs, nu_obs, flux_func);
@@ -377,8 +377,8 @@ auto PyModel::flux_density_exposures(PyArray const& t, PyArray const& nu, PyArra
 
     const auto sampling = generate_exposure_sampling(t, nu, expo_time, num_points);
 
-    auto flux_func = [](Observer& obs, Array const& t, Array const& nu, auto& photons) -> XTArray {
-        return obs.specific_flux_series(t, nu, photons) / unit::flux_den_cgs;
+    auto flux_func = [](Observer& obs, Array const& time, Array const& freq, auto& photons) -> XTArray {
+        return obs.specific_flux_series(time, freq, photons) / unit::flux_den_cgs;
     };
 
     auto result = compute_emission(sampling.t_obs_sorted, sampling.nu_obs_sorted, flux_func);
@@ -395,8 +395,8 @@ auto PyModel::flux_density_grid(PyArray const& t, PyArray const& nu) -> PyFlux {
     const Array t_obs = t * unit::sec;
     const Array nu_obs = nu * unit::Hz;
 
-    auto flux_func = [](Observer& obs, Array const& t, Array const& nu, auto& photons) -> XTArray {
-        return obs.specific_flux(t, nu, photons) / unit::flux_den_cgs;
+    auto flux_func = [](Observer& obs, Array const& time, Array const& freq, auto& photons) -> XTArray {
+        return obs.specific_flux(time, freq, photons) / unit::flux_den_cgs;
     };
 
     auto result = compute_emission(t_obs, nu_obs, flux_func);