Skip to content

Commit

Permalink
Replace doppler_packet_nucmf_on_nurf() with calculate_doppler_nucmf_o…
Browse files Browse the repository at this point in the history 
…n_nurf()
  • Loading branch information
@lukeshingles
lukeshingles committed Sep 1, 2024
1 parent 9de7eca commit 2bc0d62
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 56 deletions.
12 changes: 6 additions & 6 deletions gammapkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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));

Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion packet.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
14 changes: 7 additions & 7 deletions rpkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ auto get_nu_cmf_abort(const std::array<double, 3> pos, const std::array<double,
pos[1] + (dir[1] * half_abort_dist) + (dir[1] * half_abort_dist),
pos[2] + (dir[2] * half_abort_dist) + (dir[2] * half_abort_dist)};

const double nu_cmf_abort = nu_rf * doppler_packet_nucmf_on_nurf(abort_pos, dir, abort_time);
const double nu_cmf_abort = nu_rf * calculate_doppler_nucmf_on_nurf(abort_pos, dir, abort_time);

return nu_cmf_abort;
}
Expand Down Expand Up @@ -417,15 +417,15 @@ void electron_scatter_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;
}

void rpkt_event_continuum(Packet &pkt, const Rpkt_continuum_absorptioncoeffs &chi_rpkt_cont, const int modelgridindex) {
const double nu = pkt.nu_cmf;

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);
const double chi_cont = chi_rpkt_cont.total * dopplerfactor;
const double chi_escatter = chi_rpkt_cont.ffescat * dopplerfactor;
const double chi_ff = chi_rpkt_cont.ffheat * dopplerfactor;
Expand Down Expand Up @@ -727,8 +727,8 @@ auto do_rpkt_step(Packet &pkt, const double t2) -> 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;
Expand All @@ -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(
Expand Down Expand Up @@ -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;

Expand Down
76 changes: 37 additions & 39 deletions vectors.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,44 +71,15 @@ template <size_t S1, size_t S2>
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<double, 3> &dir_rf,
const std::array<double, 3> &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<double, 3> &pos_rf,
const std::array<double, 3> &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<double, 3> &pos_rf,
const std::array<double, 3> &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);

Expand All @@ -126,10 +97,37 @@ template <size_t S1, size_t S2>
return dopplerfactorsq;
}

[[nodiscard]] constexpr auto doppler_packet_nucmf_on_nurf(const std::array<double, 3> &pos_rf,
const std::array<double, 3> &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<double, 3> &pos_rf,
const std::array<double, 3> &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.
Expand All @@ -147,7 +145,7 @@ constexpr auto move_pkt_withtime(std::array<double, 3> &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;
Expand Down
4 changes: 1 addition & 3 deletions vpkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;

Expand Down

0 comments on commit 2bc0d62

Please sign in to comment.