Skip to content

Commit

Permalink
Remove const and pure attributes
Browse files Browse the repository at this point in the history
  • Loading branch information
lukeshingles committed Sep 1, 2024
1 parent 50974c3 commit a815dca
Showing 1 changed file with 25 additions and 29 deletions.
54 changes: 25 additions & 29 deletions vectors.h
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,14 @@

// return the the magnitude of a vector
template <size_t VECDIM>
[[nodiscard]] [[gnu::const]] constexpr auto vec_len(const std::array<double, VECDIM> vec) -> double {
[[nodiscard]] constexpr auto vec_len(const std::array<double, VECDIM> vec) -> double {
const double squaredlen = std::accumulate(vec.begin(), vec.end(), 0., [](auto a, auto b) { return a + b * b; });

return std::sqrt(squaredlen);
}

// get a normalized copy of vec_in
[[nodiscard]] [[gnu::const]] constexpr auto vec_norm(const std::array<double, 3> vec_in) {
[[nodiscard]] constexpr auto vec_norm(const std::array<double, 3> vec_in) {
const double magnitude = vec_len(vec_in);
const std::array<double, 3> vec_out{vec_in[0] / magnitude, vec_in[1] / magnitude, vec_in[2] / magnitude};

Expand All @@ -33,35 +33,32 @@ template <size_t VECDIM>

// vector dot product
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]] constexpr auto dot(const std::array<double, S1> x, const std::array<double, S2> y) -> double {
// if len(x) < len(y), the extra elements of y are ignored
return std::inner_product(x.begin(), x.end(), y.begin(), 0.);
}

// Get velocity vector of the flow at a position with homologous expansion.
[[nodiscard]] [[gnu::pure]] constexpr auto get_velocity(std::span<const double, 3> x,
const double t) -> std::array<double, 3> {
[[nodiscard]] constexpr auto get_velocity(std::span<const double, 3> x, const double t) -> std::array<double, 3> {
return std::array<double, 3>{x[0] / t, x[1] / t, x[2] / t};
}

[[nodiscard]] [[gnu::const]] constexpr auto cross_prod(const std::array<double, 3> vec_a,
const std::array<double, 3> vec_b) {
[[nodiscard]] constexpr auto cross_prod(const std::array<double, 3> vec_a, const std::array<double, 3> vec_b) {
return std::array<double, 3>{(vec_a[1] * vec_b[2]) - (vec_b[1] * vec_a[2]),
(vec_a[2] * vec_b[0]) - (vec_b[2] * vec_a[0]),
(vec_a[0] * vec_b[1]) - (vec_b[0] * vec_a[1])};
}

[[nodiscard]] [[gnu::const]] constexpr auto vec_scale(const std::array<double, 3> vec, const double scalefactor) {
[[nodiscard]] constexpr auto vec_scale(const std::array<double, 3> vec, const double scalefactor) {
return std::array<double, 3>{vec[0] * scalefactor, vec[1] * scalefactor, vec[2] * scalefactor};
}

// aberation of angles in special relativity
// dir1: direction unit vector in frame1
// vel: velocity of frame2 relative to frame1
// dir2: direction vector in frame2
[[nodiscard]] [[gnu::const]] constexpr auto angle_ab(const std::array<double, 3> dir1,
const std::array<double, 3> vel) -> std::array<double, 3> {
[[nodiscard]] constexpr auto angle_ab(const std::array<double, 3> dir1,
const std::array<double, 3> vel) -> std::array<double, 3> {
const double vsqr = dot(vel, vel) / CLIGHTSQUARED;
const double gamma_rel = 1. / std::sqrt(1 - vsqr);

Expand All @@ -80,8 +77,8 @@ template <size_t S1, size_t S2>
// 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
[[gnu::const]] [[nodiscard]] constexpr auto doppler_nucmf_on_nurf(const std::array<double, 3> dir_rf,
const std::array<double, 3> vel_rf) -> double {
[[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.);

Expand All @@ -101,9 +98,9 @@ template <size_t S1, size_t S2>
return dopplerfactor;
}

[[gnu::const]] [[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
[[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
//
Expand All @@ -130,9 +127,9 @@ template <size_t S1, size_t S2>
return dopplerfactorsq;
}

[[gnu::pure]] [[nodiscard]] constexpr auto doppler_packet_nucmf_on_nurf(const std::span<const double, 3> pos_rf,
const std::array<double, 3> dir_rf,
const double prop_time) -> double {
[[nodiscard]] constexpr auto doppler_packet_nucmf_on_nurf(const std::span<const 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));
}

Expand Down Expand Up @@ -167,16 +164,16 @@ 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
[[nodiscard]] constexpr auto get_arrive_time(const Packet &pkt) -> double
// We know that a packet escaped at "escape_time". However, we have
// to allow for travel time. Use the formula in Leon's paper. The extra
// distance to be travelled beyond the reference surface is ds = r_ref (1 - mu).
{
return pkt.escape_time - (dot(pkt.pos, pkt.dir) / CLIGHT_PROP);
}

[[nodiscard]] [[gnu::const]] constexpr auto get_escapedirectionbin(const std::array<double, 3> dir_in,
const std::array<double, 3> syn_dir) -> int {
[[nodiscard]] constexpr auto get_escapedirectionbin(const std::array<double, 3> dir_in,
const std::array<double, 3> syn_dir) -> int {
constexpr auto xhat = std::array<double, 3>{1.0, 0.0, 0.0};

// sometimes dir vectors aren't accurately normalised
Expand Down Expand Up @@ -220,10 +217,9 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
}

// Rotation angle from the scattering plane
[[nodiscard]] [[gnu::const]] constexpr auto get_rot_angle(const std::array<double, 3> n1,
const std::array<double, 3> n2,
const std::array<double, 3> ref1,
const std::array<double, 3> ref2) -> double {
[[nodiscard]] constexpr auto get_rot_angle(const std::array<double, 3> n1, const std::array<double, 3> n2,
const std::array<double, 3> ref1,
const std::array<double, 3> ref2) -> double {
// We need to rotate Stokes Parameters to (or from) the scattering plane from (or to)
// the meridian frame such that Q=1 is in the scattering plane and along ref1

Expand Down Expand Up @@ -256,7 +252,7 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
}

// Routine to compute the meridian frame axes ref1 and ref2
[[nodiscard]] [[gnu::const]] constexpr auto meridian(const std::array<double, 3> n)
[[nodiscard]] constexpr auto meridian(const std::array<double, 3> n)
-> std::tuple<std::array<double, 3>, std::array<double, 3>> {
// for ref_1 use (from triple product rule)
const double n_xylen = std::sqrt(n[0] * n[0] + n[1] * n[1]);
Expand All @@ -268,8 +264,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
return {ref1, ref2};
}

[[nodiscard]] [[gnu::const]] constexpr auto lorentz(const std::array<double, 3> e_rf, const std::array<double, 3> n_rf,
const std::array<double, 3> v) -> std::array<double, 3> {
[[nodiscard]] constexpr auto lorentz(const std::array<double, 3> e_rf, const std::array<double, 3> n_rf,
const std::array<double, 3> v) -> std::array<double, 3> {
// Use Lorentz transformations to get e_cmf from e_rf

const auto beta = std::array<double, 3>{v[0] / CLIGHT, v[1] / CLIGHT, v[2] / CLIGHT};
Expand Down

0 comments on commit a815dca

Please sign in to comment.