move_pkt_withtime: Use template parameter for CLIGHT_ASSUMED

gammapkt.cc

@@ -409,8 +409,8 @@ static auto thomson_angle() -> double {
   return mu;
 }
 
-[[nodiscard]] static auto scatter_dir(const std::array<double, 3> dir_in, const double cos_theta)
-    -> std::array<double, 3>
+[[nodiscard]] static auto scatter_dir(const std::array<double, 3> dir_in,
+                                      const double cos_theta) -> std::array<double, 3>
 // Routine for scattering a direction through angle theta.
 {
   // begin with setting the direction in coordinates where original direction
@@ -1000,15 +1000,15 @@ void wollaeger_thermalisation(Packet &pkt) {
   bool end_packet = false;
   while (!end_packet) {
     // distance to the next cell
-    const auto [sdist, snext] =
-        grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross, true);
+    const auto [sdist, snext] = grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where,
+                                                        &pkt_copy.last_cross, true);
     const double s_cont = sdist * t_current * t_current * t_current / std::pow(pkt_copy.prop_time, 3);
     const int mgi = grid::get_cell_modelgridindex(pkt_copy.where);
     if (mgi != grid::get_npts_model()) {
       tau += grid::get_rho(mgi) * s_cont * mean_gamma_opac;  // contribution to the integral
     }
     // move packet copy now
-    move_pkt_withtime(pkt_copy, sdist, true);
+    move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist);
 
     grid::change_cell(pkt_copy, snext);
     end_packet = (pkt_copy.type == TYPE_ESCAPE);
@@ -1070,7 +1070,7 @@ void guttman_thermalisation(Packet &pkt) {
         column_densities[i] += grid::get_rho_tmin(mgi) * s_cont;  // contribution to the integral
       }
       // move packet copy now
-      move_pkt_withtime(pkt_copy, sdist, true);
+      move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist);
 
       grid::change_cell(pkt_copy, snext);
       end_packet = (pkt_copy.type == TYPE_ESCAPE);

grid.cc

@@ -2462,13 +2462,11 @@ template <size_t S1>
   return -1.;
 }
 
-static auto get_coordboundary_distances_cylindrical2d(const std::array<double, 3> pkt_pos,
-                                                      const std::array<double, 3> pkt_dir,
-                                                      const std::array<double, 3> pktposgridcoord,
-                                                      const std::array<double, 3> pktvelgridcoord, int cellindex,
-                                                      const double tstart, const std::array<double, 3> cellcoordmax,
-                                                      bool neglect_exp)
-    -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
+static auto get_coordboundary_distances_cylindrical2d(
+    const std::array<double, 3> pkt_pos, const std::array<double, 3> pkt_dir,
+    const std::array<double, 3> pktposgridcoord, const std::array<double, 3> pktvelgridcoord, int cellindex,
+    const double tstart, const std::array<double, 3> cellcoordmax,
+    bool neglect_exp) -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
   // to get the cylindrical intersection, get the spherical intersection with Z components set to zero, and the
   // propagation speed set to the xy component of the 3-velocity
   double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP;

grid.h

@@ -130,8 +130,8 @@ void write_grid_restart_data(int timestep);
 [[nodiscard]] auto get_ndo_nonempty(int rank) -> int;
 [[nodiscard]] auto get_totmassradionuclide(int z, int a) -> double;
 [[nodiscard]] auto boundary_distance(std::array<double, 3> dir, std::array<double, 3> pos, double tstart, int cellindex,
-                                     enum cell_boundary *pkt_last_cross, bool neglect_exp = false)
-    -> std::tuple<double, int>;
+                                     enum cell_boundary *pkt_last_cross,
+                                     bool neglect_exp = false) -> std::tuple<double, int>;
 
 [[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float
 // mass fraction of an element (all isotopes combined)

vectors.h

@@ -34,15 +34,15 @@ template <size_t VECDIM>
 }
 
 template <size_t S1, size_t S2>
-[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array<double, S1> x, const std::array<double, S2> y)
-    -> double
+[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array<double, S1> x,
+                                                const std::array<double, S2> y) -> double
 // vector dot product
 {
   return std::inner_product(x.begin(), x.end(), y.begin(), 0.);
 }
 
-[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span<const double, 3> x, const double t)
-    -> std::array<double, 3>
+[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span<const double, 3> x,
+                                                        const double t) -> std::array<double, 3>
 // Routine for getting velocity vector of the flow at a position with homologous expansion.
 {
   return std::array<double, 3>{x[0] / t, x[1] / t, x[2] / t};
@@ -59,8 +59,8 @@ template <size_t S1, size_t S2>
   return std::array<double, 3>{vec[0] * scalefactor, vec[1] * scalefactor, vec[2] * scalefactor};
 }
 
-[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array<double, 3> dir1, const std::array<double, 3> vel)
-    -> std::array<double, 3>
+[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array<double, 3> dir1,
+                                                     const std::array<double, 3> vel) -> std::array<double, 3>
 // aberation of angles in special relativity
 //   dir1: direction unit vector in frame1
 //   vel: velocity of frame2 relative to frame1
@@ -141,17 +141,17 @@ template <size_t S1, size_t S2>
   return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time));
 }
 
+template <double CLIGHT_ASSUMED = CLIGHT_PROP>
 constexpr auto move_pkt_withtime(std::span<double, 3> pos_rf, const std::array<double, 3> dir_rf, double &prop_time,
                                  const double nu_rf, double &nu_cmf, const double e_rf, double &e_cmf,
-                                 const double distance, bool neglect_exp = false) -> double
+                                 const double distance) -> double
 /// Subroutine to move a packet along a straight line (specified by current
 /// dir vector). The distance moved is in the rest frame.
 {
-  double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP;
   assert_always(distance >= 0);
 
   const double nu_cmf_old = nu_cmf;
-  prop_time += distance / CLIGHT_PROP_DIST_CALC;
+  prop_time += distance / CLIGHT_ASSUMED;
 
   pos_rf[0] += (dir_rf[0] * distance);
   pos_rf[1] += (dir_rf[1] * distance);
@@ -171,9 +171,10 @@ constexpr auto move_pkt_withtime(std::span<double, 3> pos_rf, const std::array<d
   return dopplerfactor;
 }
 
-constexpr auto move_pkt_withtime(Packet &pkt, const double distance, bool neglect_exp = false) -> double {
-  return move_pkt_withtime(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf, distance,
-                           neglect_exp);
+template <double CLIGHT_ASSUMED = CLIGHT_PROP>
+constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
+  return move_pkt_withtime<CLIGHT_ASSUMED>(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf,
+                                           distance);
 }
 
 [[nodiscard]] [[gnu::const]] constexpr auto get_arrive_time(const Packet &pkt) -> double
@@ -310,8 +311,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance, bool neglec
 }
 
 // Routine to transform the Stokes Parameters from RF to CMF
-constexpr auto frame_transform(const std::array<double, 3> n_rf, double *Q, double *U, const std::array<double, 3> v)
-    -> std::array<double, 3> {
+constexpr auto frame_transform(const std::array<double, 3> n_rf, double *Q, double *U,
+                               const std::array<double, 3> v) -> std::array<double, 3> {
   // Meridian frame in the RF
   const auto [ref1_rf, ref2_rf] = meridian(n_rf);