From 2bc0d624e70a707582c11803f960b4f129dca3da Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Sun, 1 Sep 2024 23:08:17 +0100 Subject: [PATCH] Replace doppler_packet_nucmf_on_nurf() with calculate_doppler_nucmf_on_nurf() --- gammapkt.cc | 12 ++++----- packet.cc | 2 +- rpkt.cc | 14 +++++----- vectors.h | 76 ++++++++++++++++++++++++++--------------------------- vpkt.cc | 4 +-- 5 files changed, 52 insertions(+), 56 deletions(-) diff --git a/gammapkt.cc b/gammapkt.cc index 7f3e2f053..112cb05ac 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -277,7 +277,7 @@ auto get_chi_compton_rf(const Packet &pkt) -> double { const double chi_cmf = sigma_cmf * grid::get_nnetot(grid::get_cell_modelgridindex(pkt.where)); // convert between frames - const double chi_rf = chi_cmf * doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double chi_rf = chi_cmf * calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); assert_testmodeonly(std::isfinite(chi_rf)); @@ -460,7 +460,7 @@ void compton_scatter(Packet &pkt) { // It now has a rest frame direction and a co-moving frequency. // Just need to set the rest frame energy. - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); pkt.nu_rf = pkt.nu_cmf / dopplerfactor; pkt.e_rf = pkt.e_cmf / dopplerfactor; @@ -566,7 +566,7 @@ auto get_chi_photo_electric_rf(const Packet &pkt) -> double { // Now convert between frames. - const double chi_rf = chi_cmf * doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double chi_rf = chi_cmf * calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); return chi_rf; } @@ -622,7 +622,7 @@ auto sigma_pair_prod_rf(const Packet &pkt) -> double { // Now need to convert between frames. - double chi_rf = chi_cmf * doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + double chi_rf = chi_cmf * calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); if (chi_rf < 0) { printout("Negative pair production sigma. Setting to zero. Abort? %g\n", chi_rf); @@ -717,7 +717,7 @@ void pair_prod(Packet &pkt) { pkt.dir = angle_ab(dir_cmf, vel_vec); - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); pkt.nu_rf = pkt.nu_cmf / dopplerfactor; pkt.e_rf = pkt.e_cmf / dopplerfactor; @@ -1066,7 +1066,7 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) { // that it's now a gamma ray. pkt.prop_time = pkt.tdecay; - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); pkt.nu_rf = pkt.nu_cmf / dopplerfactor; pkt.e_rf = pkt.e_cmf / dopplerfactor; diff --git a/packet.cc b/packet.cc index c06d5e8a5..12d8a4768 100644 --- a/packet.cc +++ b/packet.cc @@ -77,7 +77,7 @@ void place_pellet(const double e0, const int cellindex, const int pktnumber, Pac // pellet packet is moving with the homologous flow, so dir is proportional to pos pkt.dir = vec_norm(pkt.pos); // assign dir = pos / vec_len(pos) - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); pkt.e_rf = pkt.e_cmf / dopplerfactor; pkt.trueemissiontype = EMTYPE_NOTSET; diff --git a/rpkt.cc b/rpkt.cc index 068bef2c6..8127663e8 100644 --- a/rpkt.cc +++ b/rpkt.cc @@ -59,7 +59,7 @@ auto get_nu_cmf_abort(const std::array pos, const std::array bool { } else if (thickcell) [[unlikely]] { // In the case of optically thick cells, we treat the packets in grey approximation to speed up the calculation - const double chi_grey = - grid::get_kappagrey(mgi) * grid::get_rho(mgi) * doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double chi_grey = grid::get_kappagrey(mgi) * grid::get_rho(mgi) * + calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); edist = tau_next / chi_grey; pkt.next_trans = -1; @@ -740,7 +740,7 @@ auto do_rpkt_step(Packet &pkt, const double t2) -> bool { const auto nu_cmf_abort = get_nu_cmf_abort(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, abort_dist); const auto d_nu_on_d_l = (nu_cmf_abort - pkt.nu_cmf) / abort_dist; - const auto doppler = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const auto doppler = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); if constexpr (EXPANSIONOPACITIES_ON) { std::tie(edist, pkt.next_trans, event_is_boundbound) = get_event_expansion_opacity( @@ -1111,7 +1111,7 @@ __host__ __device__ void emit_rpkt(Packet &pkt) { // Finally we want to put in the rest frame energy and frequency. And record // that it's now a r-pkt. - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time); pkt.nu_rf = pkt.nu_cmf / dopplerfactor; pkt.e_rf = pkt.e_cmf / dopplerfactor; diff --git a/vectors.h b/vectors.h index 39014eb9a..d0a5958c6 100644 --- a/vectors.h +++ b/vectors.h @@ -71,44 +71,15 @@ template return vec_norm(dir2); } -// Doppler factor -// arguments: -// dir_rf: the rest frame direction (unit vector) of light propagation -// vel_rf: velocity of the comoving frame relative to the rest frame -// returns: the ratio f = nu_cmf / nu_rf -[[nodiscard]] constexpr auto doppler_nucmf_on_nurf(const std::array &dir_rf, - const std::array &vel_rf) -> double { - assert_testmodeonly(dot(vel_rf, vel_rf) / CLIGHTSQUARED >= 0.); - assert_testmodeonly(dot(vel_rf, vel_rf) / CLIGHTSQUARED < 1.); - - const double ndotv = dot(dir_rf, vel_rf); - double dopplerfactor = 1. - (ndotv / CLIGHT); - - if (USE_RELATIVISTIC_DOPPLER_SHIFT) { - const double betasq = dot(vel_rf, vel_rf) / CLIGHTSQUARED; - assert_testmodeonly(betasq >= 0.); // v < c - assert_testmodeonly(betasq < 1.); // v < c - dopplerfactor = dopplerfactor / std::sqrt(1 - betasq); - } - - assert_testmodeonly(std::isfinite(dopplerfactor)); - assert_testmodeonly(dopplerfactor > 0); - - return dopplerfactor; -} - -[[nodiscard]] constexpr auto doppler_squared_nucmf_on_nurf(const std::array &pos_rf, - const std::array &dir_rf, - const double prop_time) -> double -// Doppler factor squared, either to first order v/c or fully relativisitic -// depending on USE_RELATIVISTIC_DOPPLER_SHIFT -// -// arguments: +// Doppler factor squared, either to first order v/c or fully relativisitic depending on USE_RELATIVISTIC_DOPPLER_SHIFT +// Arguments: // pos_rf: the rest frame position of the packet // dir_rf: the rest frame direction (unit vector) of light propagation // prop_time: the propagation time of the packet // returns: the ratio f = (nu_cmf / nu_rf) ^ 2 -{ +[[nodiscard]] constexpr auto doppler_squared_nucmf_on_nurf(const std::array &pos_rf, + const std::array &dir_rf, + const double prop_time) -> double { // velocity of the comoving frame relative to the rest frame const auto vel_rf = get_velocity(pos_rf, prop_time); @@ -126,10 +97,37 @@ template return dopplerfactorsq; } -[[nodiscard]] constexpr auto doppler_packet_nucmf_on_nurf(const std::array &pos_rf, - const std::array &dir_rf, - const double prop_time) -> double { - return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time)); +// Doppler factor either to first order v/c or fully relativisitic depending on USE_RELATIVISTIC_DOPPLER_SHIFT +// Arguments: +// pos_rf: the rest frame position of the packet +// dir_rf: the rest frame direction (unit vector) of light propagation +// prop_time: the propagation time of the packet +// returns: the ratio f = nu_cmf / nu_rf +[[nodiscard]] constexpr auto calculate_doppler_nucmf_on_nurf(const std::array &pos_rf, + const std::array &dir_rf, + const double prop_time) -> double { + // dir_rf: the rest frame direction (unit vector) of light propagation + // vel_rf: velocity of the comoving frame relative to the rest frame + // returns: the ratio f = nu_cmf / nu_rf + const auto vel_rf = get_velocity(pos_rf, prop_time); + + assert_testmodeonly(dot(vel_rf, vel_rf) / CLIGHTSQUARED >= 0.); + assert_testmodeonly(dot(vel_rf, vel_rf) / CLIGHTSQUARED < 1.); + + const double ndotv = dot(dir_rf, vel_rf); + double dopplerfactor = 1. - (ndotv / CLIGHT); + + if (USE_RELATIVISTIC_DOPPLER_SHIFT) { + const double betasq = dot(vel_rf, vel_rf) / CLIGHTSQUARED; + assert_testmodeonly(betasq >= 0.); // v < c + assert_testmodeonly(betasq < 1.); // v < c + dopplerfactor = dopplerfactor / std::sqrt(1 - betasq); + } + + assert_testmodeonly(std::isfinite(dopplerfactor)); + assert_testmodeonly(dopplerfactor > 0); + + return dopplerfactor; } // Move a packet along a straight line (specified by current dir vector). The distance moved is in the rest frame. @@ -147,7 +145,7 @@ constexpr auto move_pkt_withtime(std::array &pos_rf, const std::array // During motion, rest frame energy and frequency are conserved. // But need to update the co-moving ones. - const double dopplerfactor = doppler_packet_nucmf_on_nurf(pos_rf, dir_rf, prop_time); + const double dopplerfactor = calculate_doppler_nucmf_on_nurf(pos_rf, dir_rf, prop_time); nu_cmf = nu_rf * dopplerfactor; e_cmf = e_rf * dopplerfactor; diff --git a/vpkt.cc b/vpkt.cc index b2acb55ea..047c1147d 100644 --- a/vpkt.cc +++ b/vpkt.cc @@ -908,8 +908,6 @@ auto vpkt_call_estimators(Packet &pkt, const enum packet_type type_before_rpkt) const double t_current = pkt.prop_time; - const auto vel_vec = get_velocity(pkt.pos, pkt.prop_time); - std::stringstream vpkt_contrib_row; bool any_dir_escaped = false; @@ -927,7 +925,7 @@ auto vpkt_call_estimators(Packet &pkt, const enum packet_type type_before_rpkt) if (t_arrive >= VSPEC_TIMEMIN_input && t_arrive <= VSPEC_TIMEMAX_input) { // time selection - const double doppler = doppler_nucmf_on_nurf(obsdir, vel_vec); + const double doppler = calculate_doppler_nucmf_on_nurf(pkt.pos, obsdir, pkt.prop_time); const double nu_rf = pkt.nu_cmf / doppler; const double e_rf = pkt.e_cmf / doppler;