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);

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);