From 08d8c4c14041f69bb46c02d4d736ac0442b610c4 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Fri, 16 Aug 2024 12:00:27 +0100 Subject: [PATCH] move_pkt_withtime: Use template parameter for CLIGHT_ASSUMED --- gammapkt.cc | 12 ++++++------ grid.cc | 12 +++++------- grid.h | 4 ++-- vectors.h | 29 +++++++++++++++-------------- 4 files changed, 28 insertions(+), 29 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 1134e30c6..d66b989a5 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -409,8 +409,8 @@ static auto thomson_angle() -> double { return mu; } -[[nodiscard]] static auto scatter_dir(const std::array dir_in, const double cos_theta) - -> std::array +[[nodiscard]] static auto scatter_dir(const std::array dir_in, + const double cos_theta) -> std::array // 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); diff --git a/grid.cc b/grid.cc index b0107847c..132be6ed6 100644 --- a/grid.cc +++ b/grid.cc @@ -2462,13 +2462,11 @@ template return -1.; } -static auto get_coordboundary_distances_cylindrical2d(const std::array pkt_pos, - const std::array pkt_dir, - const std::array pktposgridcoord, - const std::array pktvelgridcoord, int cellindex, - const double tstart, const std::array cellcoordmax, - bool neglect_exp) - -> std::tuple, std::array> { +static auto get_coordboundary_distances_cylindrical2d( + const std::array pkt_pos, const std::array pkt_dir, + const std::array pktposgridcoord, const std::array pktvelgridcoord, int cellindex, + const double tstart, const std::array cellcoordmax, + bool neglect_exp) -> std::tuple, std::array> { // 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; diff --git a/grid.h b/grid.h index 95db53000..f09a6f019 100644 --- a/grid.h +++ b/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 dir, std::array pos, double tstart, int cellindex, - enum cell_boundary *pkt_last_cross, bool neglect_exp = false) - -> std::tuple; + enum cell_boundary *pkt_last_cross, + bool neglect_exp = false) -> std::tuple; [[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float // mass fraction of an element (all isotopes combined) diff --git a/vectors.h b/vectors.h index e70be76f1..b21332c49 100644 --- a/vectors.h +++ b/vectors.h @@ -34,15 +34,15 @@ template } template -[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, const std::array y) - -> double +[[nodiscard]] [[gnu::const]] constexpr auto dot(const std::array x, + const std::array 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 x, const double t) - -> std::array +[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span x, + const double t) -> std::array // Routine for getting velocity vector of the flow at a position with homologous expansion. { return std::array{x[0] / t, x[1] / t, x[2] / t}; @@ -59,8 +59,8 @@ template return std::array{vec[0] * scalefactor, vec[1] * scalefactor, vec[2] * scalefactor}; } -[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, const std::array vel) - -> std::array +[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array dir1, + const std::array vel) -> std::array // aberation of angles in special relativity // dir1: direction unit vector in frame1 // vel: velocity of frame2 relative to frame1 @@ -141,17 +141,17 @@ template return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time)); } +template constexpr auto move_pkt_withtime(std::span pos_rf, const std::array 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 pos_rf, const std::array 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 +constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> 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); } [[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 n_rf, double *Q, double *U, const std::array v) - -> std::array { +constexpr auto frame_transform(const std::array n_rf, double *Q, double *U, + const std::array v) -> std::array { // Meridian frame in the RF const auto [ref1_rf, ref2_rf] = meridian(n_rf);