code cleanup

YihanWangAstro · YihanWangAstro · commit 7695610d984e · 2025-10-12T00:10:27.000-05:00
diff --git a/include/mesh.h b/include/mesh.h
@@ -412,7 +412,7 @@ Coord auto_grid(Ejecta const& jet, Array const& t_obs, Real theta_cut, Real thet
     Real theta_max = std::min(jet_edge, theta_cut);
 
     size_t theta_num = std::max<size_t>(static_cast<size_t>((theta_max - theta_min) * 180 / con::pi * theta_resol), 56);
-    const size_t uniform_theta_num = static_cast<size_t>(theta_num * 0.3);
+    const size_t uniform_theta_num = static_cast<size_t>(static_cast<Real>(theta_num) * 0.3);
     size_t adaptive_theta_num = theta_num - uniform_theta_num;
 
     const Array uniform_theta = xt::linspace(theta_min, theta_max, uniform_theta_num);
@@ -422,7 +422,7 @@ Coord auto_grid(Ejecta const& jet, Array const& t_obs, Real theta_cut, Real thet
     // coord.theta = uniform_theta;
 
     const size_t phi_num = std::max<size_t>(static_cast<size_t>(360 * phi_resol), 1);
-    const size_t uniform_phi_num = static_cast<size_t>(phi_num * 0.3);
+    const size_t uniform_phi_num = static_cast<size_t>(static_cast<Real>(phi_num) * 0.3);
     const size_t adaptive_phi_num = phi_num - uniform_phi_num;
 
     const Array uniform_phi = xt::linspace(0., 2 * con::pi, uniform_phi_num);
diff --git a/include/observer.h b/include/observer.h
@@ -387,9 +387,8 @@ Array Observer::specific_flux_series(Array const& t_obs, Array const& nu_obs, Ph
             size_t t_idx = 0;
             iterate_to(time(i, j, 0), t_obs, t_idx);
 
-            size_t k = 0;
             InterpState state;
-            while (t_idx < t_obs_len && k < t_grid - 1) {
+            for (size_t k = 0; t_idx < t_obs_len && k < t_grid - 1;) {
                 if (time(i, j, k + 1) < t_obs(t_idx)) {
                     k++;
                 } else {
diff --git a/include/simple-shock.hpp b/include/simple-shock.hpp
@@ -16,7 +16,7 @@
  * @brief Represents the state vector for the simple shock equation.
  * <!-- ************************************************************************************** -->
  */
-template <typename Ejecta, typename Medium>
+template <typename Ejecta>
 struct SimpleState {
     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
@@ -55,7 +55,7 @@ struct SimpleState {
 template <typename Ejecta, typename Medium>
 class SimpleShockEqn {
   public:
-    using State = SimpleState<Ejecta, Medium>;
+    using State = SimpleState<Ejecta>;
 
     /**
      * <!-- ************************************************************************************** -->
diff --git a/include/synchrotron.h b/include/synchrotron.h
@@ -231,17 +231,17 @@ void update_electrons_4Y(SynElectronGrid& electrons, Shock const& shock);
 
 /**
  * <!-- ************************************************************************************** -->
- * @brief Calculates cooling Lorentz factor based on comoving time, magnetic field, and IC parameters
+ * @brief Calculates a cooling Lorentz factor based on comoving time, magnetic field, and IC parameters
  * @details Accounts for synchrotron and inverse Compton cooling using an iterative approach
- *          to handle the Lorentz factor dependent IC cooling.
- * @param t_com Comoving time
+ *          to handle the Lorentz factor-dependent IC cooling.
+ * @param t_comv Comoving time
  * @param B Magnetic field
  * @param Ys Inverse Compton Y parameters
  * @param p Power-law index for electron energy distribution
  * @return The cooling Lorentz factor
  * <!-- ************************************************************************************** -->
  */
-Real compute_gamma_c(Real t_com, Real B, InverseComptonY const& Ys, Real p);
+Real compute_gamma_c(Real t_comv, Real B, InverseComptonY const& Ys, Real p);
 
 /**
  * <!-- ************************************************************************************** -->
diff --git a/pybind/mcmc.cpp b/pybind/mcmc.cpp
@@ -29,15 +29,15 @@ std::vector<size_t> MultiBandData::logscale_screen(PyArray const& data, size_t n
     const double log_start = std::log10(static_cast<double>(data(0)));
     const double log_end = std::log10(static_cast<double>(data(total_size - 1)));
     const double log_range = log_end - log_start;
-    const size_t total_points = static_cast<size_t>(std::ceil(log_range * num_order)) + 1;
+    const size_t total_points = static_cast<size_t>(std::ceil(log_range * static_cast<double>(num_order))) + 1;
 
     std::vector<size_t> indices;
     indices.reserve(total_points);
 
     // Always include first point
     indices.push_back(0);
 
-    const double step = log_range / (total_points - 1);
+    const double step = log_range / static_cast<double>(total_points - 1);
 
     for (size_t i = 1; i < total_points - 1; ++i) {
         const double log_target = log_start + i * step;
@@ -59,7 +59,7 @@ std::vector<size_t> MultiBandData::logscale_screen(PyArray const& data, size_t n
         }
     }
 
-    // Always include last point
+    // Always include the last point
     if (total_size > 1) {
         indices.push_back(total_size - 1);
     }
@@ -152,7 +152,7 @@ Medium MultiBandModel::select_medium(Params const& param) const {
 }
 
 void MultiBandData::add_flux_density(double nu, PyArray const& t, PyArray const& Fv_obs, PyArray const& Fv_err,
-                                     std::optional<PyArray> weights) {
+                                     std::optional<PyArray> const& weights) {
     AFTERGLOW_REQUIRE(t.size() == Fv_obs.size() && t.size() == Fv_err.size(), "light curve array inconsistent length!");
 
     Array w = xt::ones<Real>({t.size()});
@@ -169,7 +169,7 @@ void MultiBandData::add_flux_density(double nu, PyArray const& t, PyArray const&
 }
 
 void MultiBandData::add_flux(double nu_min, double nu_max, size_t num_points, PyArray const& t, PyArray const& Fv_obs,
-                             PyArray const& Fv_err, std::optional<PyArray> weights) {
+                             PyArray const& Fv_err, const std::optional<PyArray>& weights) {
     AFTERGLOW_REQUIRE(t.size() == Fv_obs.size() && t.size() == Fv_err.size(), "light curve array inconsistent length!");
     AFTERGLOW_REQUIRE(is_ascending(t), "Time array must be in ascending order!");
     AFTERGLOW_REQUIRE(nu_min < nu_max, "nu_min must be less than nu_max!");
@@ -195,7 +195,7 @@ void MultiBandData::add_flux(double nu_min, double nu_max, size_t num_points, Py
 }
 
 void MultiBandData::add_spectrum(double t, PyArray const& nu, PyArray const& Fv_obs, PyArray const& Fv_err,
-                                 std::optional<PyArray> weights) {
+                                 const std::optional<PyArray>& weights) {
     AFTERGLOW_REQUIRE(nu.size() == Fv_obs.size() && nu.size() == Fv_err.size(), "spectrum array inconsistent length!");
 
     Array w = xt::ones<Real>({nu.size()});
diff --git a/pybind/mcmc.h b/pybind/mcmc.h
@@ -82,7 +82,7 @@ struct MultiBandData {
      * <!-- ************************************************************************************** -->
      */
     void add_flux_density(double nu, PyArray const& t, PyArray const& Fv_obs, PyArray const& Fv_err,
-                          std::optional<PyArray> weights = std::nullopt);
+                          std::optional<PyArray> const& weights = std::nullopt);
 
     /**
      * <!-- ************************************************************************************** -->
@@ -100,7 +100,7 @@ struct MultiBandData {
      * <!-- ************************************************************************************** -->
      */
     void add_flux(double nu_min, double nu_max, size_t num_points, PyArray const& t, PyArray const& Fv_obs,
-                  PyArray const& Fv_err, std::optional<PyArray> weights = std::nullopt);
+                  PyArray const& Fv_err, const std::optional<PyArray>& weights = std::nullopt);
 
     /**
      * <!-- ************************************************************************************** -->
@@ -116,7 +116,7 @@ struct MultiBandData {
      * <!-- ************************************************************************************** -->
      */
     void add_spectrum(double t, PyArray const& nu, PyArray const& Fv_obs, PyArray const& Fv_err,
-                      std::optional<PyArray> weights = std::nullopt);
+                      const std::optional<PyArray>& weights = std::nullopt);
 
     /**
      * <!-- ************************************************************************************** -->
@@ -150,7 +150,7 @@ struct MultiBandData {
      * @return size_t Total number of observational data points
      * <!-- ************************************************************************************** -->
      */
-    size_t data_points_num() const;
+    [[nodiscard]] size_t data_points_num() const;
 
     Array times;            ///< Consolidated array of all observation times [seconds]
     Array frequencies;      ///< Consolidated array of all observation frequencies [Hz]
diff --git a/pybind/pymodel.cpp b/pybind/pymodel.cpp
@@ -15,7 +15,7 @@
 #include "xtensor/misc/xsort.hpp"
 
 Ejecta PyTophatJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading, Real duration,
-                   std::optional<PyMagnetar> magnetar) {
+                   const std::optional<PyMagnetar>& magnetar) {
     Ejecta jet;
     jet.eps_k = math::tophat(theta_c, E_iso);
     jet.Gamma0 = math::tophat_plus_one(theta_c, Gamma0 - 1);
@@ -30,7 +30,7 @@ Ejecta PyTophatJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading, Real d
 }
 
 Ejecta PyGaussianJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading, Real duration,
-                     std::optional<PyMagnetar> magnetar) {
+                     const std::optional<PyMagnetar>& magnetar) {
     Ejecta jet;
     jet.eps_k = math::gaussian(theta_c, E_iso);
     jet.Gamma0 = math::gaussian_plus_one(theta_c, Gamma0 - 1);
@@ -45,7 +45,7 @@ Ejecta PyGaussianJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading, Real
 }
 
 Ejecta PyPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real k_e, Real k_g, bool spreading, Real duration,
-                     std::optional<PyMagnetar> magnetar) {
+                     const std::optional<PyMagnetar>& magnetar) {
     Ejecta jet;
     jet.eps_k = math::powerlaw(theta_c, E_iso, k_e);
     jet.Gamma0 = math::powerlaw_plus_one(theta_c, Gamma0 - 1, k_g);
@@ -70,7 +70,7 @@ Ejecta PyPowerLawWing(Real theta_c, Real E_iso_w, Real Gamma0_w, Real k_e, Real
 }
 
 Ejecta PyStepPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real E_iso_w, Real Gamma0_w, Real k_e, Real k_g,
-                         bool spreading, Real duration, std::optional<PyMagnetar> magnetar) {
+                         bool spreading, Real duration, const std::optional<PyMagnetar>& magnetar) {
     Ejecta jet;
     jet.eps_k = math::step_powerlaw(theta_c, E_iso, E_iso_w, k_e);
     jet.Gamma0 = math::step_powerlaw_plus_one(theta_c, Gamma0 - 1, Gamma0_w - 1, k_g);
@@ -86,7 +86,7 @@ Ejecta PyStepPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real E_iso_w, Re
 }
 
 Ejecta PyTwoComponentJet(Real theta_c, Real E_iso, Real Gamma0, Real theta_w, Real E_iso_w, Real Gamma0_w,
-                         bool spreading, Real duration, std::optional<PyMagnetar> magnetar) {
+                         bool spreading, Real duration, const std::optional<PyMagnetar>& magnetar) {
     Ejecta jet;
     jet.eps_k = math::two_component(theta_c, theta_w, E_iso, E_iso_w);
 
@@ -319,7 +319,7 @@ auto PyModel::generate_exposure_sampling(PyArray const& t, PyArray const& nu, Py
     // Generate time-frequency samples within each exposure window
     for (size_t i = 0, j = 0; i < t.size() && j < total_points; ++i) {
         const Real t_start = t(i);
-        const Real dt = expo_time(i) / (num_points - 1);
+        const Real dt = expo_time(i) / static_cast<Real>(num_points - 1);
 
         for (size_t k = 0; k < num_points && j < total_points; ++k, ++j) {
             t_obs(j) = t_start + k * dt;
@@ -375,15 +375,15 @@ auto PyModel::flux_density_exposures(PyArray const& t, PyArray const& nu, PyArra
                       "time, frequency, and exposure time arrays must have the same size");
     AFTERGLOW_REQUIRE(num_points >= 2, "num_points must be at least 2 to sample within each exposure time");
 
-    const auto sampling = generate_exposure_sampling(t, nu, expo_time, num_points);
+    const auto [t_obs_sorted, nu_obs_sorted, idx_sorted] = generate_exposure_sampling(t, nu, expo_time, num_points);
 
     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);
+    auto result = compute_emission(t_obs_sorted, nu_obs_sorted, flux_func);
 
-    average_exposure_flux(result, sampling.idx_sorted, t.size(), num_points);
+    average_exposure_flux(result, idx_sorted, t.size(), num_points);
 
     result.calc_total();
     return result;
diff --git a/pybind/pymodel.h b/pybind/pymodel.h
@@ -62,7 +62,7 @@ struct PyMagnetar {
  * <!-- ************************************************************************************** -->
  */
 Ejecta PyTophatJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading = false, Real duration = 1,
-                   std::optional<PyMagnetar> magnetar = std::nullopt);
+                   const std::optional<PyMagnetar>& magnetar = std::nullopt);
 
 /**
  * <!-- ************************************************************************************** -->
@@ -82,7 +82,7 @@ Ejecta PyTophatJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading = false
  * <!-- ************************************************************************************** -->
  */
 Ejecta PyGaussianJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading = false, Real duration = 1,
-                     std::optional<PyMagnetar> magnetar = std::nullopt);
+                     const std::optional<PyMagnetar>& magnetar = std::nullopt);
 
 /**
  * <!-- ************************************************************************************** -->
@@ -104,7 +104,7 @@ Ejecta PyGaussianJet(Real theta_c, Real E_iso, Real Gamma0, bool spreading = fal
  * <!-- ************************************************************************************** -->
  */
 Ejecta PyPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real k_e, Real k_g, bool spreading = false,
-                     Real duration = 1, std::optional<PyMagnetar> magnetar = std::nullopt);
+                     Real duration = 1, const std::optional<PyMagnetar>& magnetar = std::nullopt);
 
 /**
  * <!-- ************************************************************************************** -->
@@ -147,7 +147,7 @@ Ejecta PyPowerLawWing(Real theta_c, Real E_iso_w, Real Gamma0_w, Real k_e, Real
  * <!-- ************************************************************************************** -->
  */
 Ejecta PyStepPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real E_iso_w, Real Gamma0_w, Real k_e, Real k_g,
-                         bool spreading, Real duration, std::optional<PyMagnetar> magnetar);
+                         bool spreading, Real duration, const std::optional<PyMagnetar>& magnetar);
 
 /**
  * <!-- ************************************************************************************** -->
@@ -169,7 +169,8 @@ Ejecta PyStepPowerLawJet(Real theta_c, Real E_iso, Real Gamma0, Real E_iso_w, Re
  * <!-- ************************************************************************************** -->
  */
 Ejecta PyTwoComponentJet(Real theta_c, Real E_iso, Real Gamma0, Real theta_w, Real E_iso_w, Real Gamma0_w,
-                         bool spreading = false, Real duration = 1, std::optional<PyMagnetar> magnetar = std::nullopt);
+                         bool spreading = false, Real duration = 1,
+                         const std::optional<PyMagnetar>& magnetar = std::nullopt);
 
 /**
  * <!-- ************************************************************************************** -->
@@ -399,13 +400,13 @@ class PyModel {
      * @param axisymmetric Whether to assume axisymmetric jet structure (default: true)
      * <!-- ************************************************************************************** -->
      */
-    PyModel(Ejecta jet, Medium medium, PyObserver observer, PyRadiation fwd_rad,
-            std::optional<PyRadiation> rvs_rad = std::nullopt,
-            std::tuple<Real, Real, Real> resolutions = std::make_tuple(0.3, 1, 10), Real rtol = 1e-5,
+    PyModel(Ejecta jet, Medium medium, const PyObserver& observer, const PyRadiation& fwd_rad,
+            const std::optional<PyRadiation>& rvs_rad = std::nullopt,
+            const std::tuple<Real, Real, Real>& resolutions = std::make_tuple(0.3, 1, 10), Real rtol = 1e-5,
             bool axisymmetric = true)
         : jet_(std::move(jet)),
           medium_(std::move(medium)),
-          obs_setup(std::move(observer)),
+          obs_setup(observer),
           fwd_rad(fwd_rad),
           rvs_rad_opt(rvs_rad),
           phi_resol(std::get<0>(resolutions)),
diff --git a/script/mcmc.ipynb b/script/mcmc.ipynb
diff --git a/script/quick.ipynb b/script/quick.ipynb
@@ -14,7 +14,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 2,
    "id": "9fb5c06d",
    "metadata": {},
    "outputs": [],
@@ -37,7 +37,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "id": "85de1437",
    "metadata": {},
    "outputs": [
diff --git a/src/dynamics/forward-shock.tpp b/src/dynamics/forward-shock.tpp
@@ -17,8 +17,7 @@ ForwardShockEqn<Ejecta, Medium>::ForwardShockEqn(Medium const& medium, Ejecta co
       theta0(theta),
       rad(rad_params),
       dOmega0(1 - std::cos(theta)),
-      theta_s(theta_s),
-      m_jet0(0) {
+      theta_s(theta_s){
     m_jet0 = ejecta.eps_k(phi, theta0) / ejecta.Gamma0(phi, theta0) / con::c2;
     if constexpr (HasSigma<Ejecta>) {
         m_jet0 /= 1 + ejecta.sigma0(phi, theta0);
@@ -175,23 +174,23 @@ void grid_solve_fwd_shock(size_t i, size_t j, View const& t, Shock& shock, FwdEq
     Real t0 = std::min(t.front(), 1 * unit::sec);
     eqn.set_init_state(state, t0);
 
-    // Early exit if initial Lorentz factor is below cutoff
+    // Early exit if the initial Lorentz factor is below cutoff
     if (state.Gamma <= con::Gamma_cut) {
         set_stopping_shock(i, j, shock, state);
         return;
     }
 
-    // Set up ODE solver with adaptive step size control
+    // Set up the ODE solver with adaptive step size control
     auto stepper = make_dense_output(rtol, rtol, runge_kutta_dopri5<typename FwdEqn::State>());
 
     stepper.initialize(state, t0, 0.01 * t0);
 
-    // Solve ODE and update shock state at each requested time point
+    // Solve ODE and update the shock state at each requested time point
     for (size_t k = 0; stepper.current_time() <= t.back();) {
         // Advance solution by one adaptive step
         stepper.do_step(eqn);
 
-        // Update shock state for all time points that have been passed in this step
+        // Update shock state for all-time points that have been passed in this step
         while (k < t.size() && stepper.current_time() > t(k)) {
             stepper.calc_state(t(k), state);
             save_fwd_shock_state(i, j, k, eqn, state, shock);
diff --git a/src/dynamics/simple-shock.tpp b/src/dynamics/simple-shock.tpp
@@ -17,8 +17,7 @@ SimpleShockEqn<Ejecta, Medium>::SimpleShockEqn(Medium const& medium, Ejecta cons
       theta0(theta),
       rad(rad_params),
       dOmega0(1 - std::cos(theta0)),
-      theta_s(theta_s),
-      m_jet0(0) {
+      theta_s(theta_s) {
     m_jet0 = ejecta.eps_k(phi, theta0) / ejecta.Gamma0(phi, theta0) / con::c2;
     if constexpr (HasSigma<Ejecta>) {
         m_jet0 /= 1 + ejecta.sigma0(phi, theta0);
diff --git a/src/radiation/inverse-compton.cpp b/src/radiation/inverse-compton.cpp
diff --git a/src/radiation/prompt.cpp b/src/radiation/prompt.cpp
diff --git a/tests/benchmark/benchmark.cpp b/tests/benchmark/benchmark.cpp
