-
Notifications
You must be signed in to change notification settings - Fork 47
integrate_backwards=True in EMRIInspiral returns inverted physical phases #183
Copy link
Copy link
Open
Labels
bugSomething isn't workingSomething isn't working
Description
From FEWv2.0.0 tag:
When generating a backward trajectory using the integrate_backwards=True flag, the initial phase signs are flipped to use the same integrator module (see few.trajectory.inspiral.EMRIInspiral.get_inspiral Line 242).
However, when returning the trajectory in Line 269, the phases are not flipped back, which results in
Reproduce the bug:
import numpy as np
from few.trajectory.inspiral import EMRIInspiral
from few.trajectory.ode import KerrEccEqFlux
from few.utils.utility import get_p_at_t
kerr_traj = EMRIInspiral(func=KerrEccEqFlux)
T = 1.0
m1, m2, a, e0, xI0 = 1e6, 1e1, 0.9, 0.2, 1.0
Phi_phi0, Phi_theta0, Phi_r0 = 1.0, 0.0, 1.0
p0 = 12.0
# 1. Forward integration
traj = kerr_traj(m1, m2, a, p0, e0, xI0,
Phi_phi0=Phi_phi0, Phi_theta0=Phi_theta0, Phi_r0=Phi_r0,
T=T, err_tol=1e-11)
# 2. Backward integration from the final state
traj_back = kerr_traj(m1, m2, a, traj[1][-1], traj[2][-1], traj[3][-1],
Phi_phi0=-traj[4][-1], # <= user must supply the negative of the final phases
Phi_theta0=-traj[5][-1],
Phi_r0=-traj[6][-1],
T=T, upsample=False, integrate_backwards=True, err_tol=1e-11)
# The raw returned phase is the negative of the true physical phase
print(f"Input Initial Phase: {Phi_phi0}")
print(f"Raw Returned t=0 Phase: {traj_back[4][0]}")
print(f"Negated Returned t=0 Phase: {-traj_back[4][0]}")Output:
Input Initial Phase: 1.0
Raw Returned t=0 Phase: -0.9999967547482811
Negated Returned t=0 Phase: 0.9999967547482811
Proposed fix:
Add a condition to negate the phase columns of the out array right before return:
if temp_kwargs["integrate_backwards"]:
out[:, 4] *= -1
out[:, 5] *= -1
out[:, 6] *= -1
t, p, e, x, Phi_phi, Phi_theta, Phi_r = out.T.copy()
return t, p, e, x, Phi_phi, Phi_theta, Phi_rReactions are currently unavailable
Metadata
Metadata
Assignees
Labels
bugSomething isn't workingSomething isn't working