Hi,
First of all, thank you for maintaining gala! I'm using version 1.11.0 and encountered two issues when integrating orbits with specific time arrays.
Code:
mw = gp.MilkyWayPotential(version="latest")
w0 = gd.PhaseSpacePosition(
pos=np.array([-8.122, 0.0, 0.0208]) * u.kpc,
vel=np.array([12.9, 245.6, 7.78]) * u.km / u.s
)
orbit = mw.integrate_orbit(w0, t=ts * u.Myr, Integrator='dop853', cython_if_possible=True)
Issue 1: Integration fails with extremely small time spacing
When ts is defined as:
ts = np.array([1e-3, 1e-3 + 1e-18])
I get the following error:
It seems that when two consecutive time points are extremely close (difference ~1e-18), the integrator fails.
Issue 2: Significant deviation with certain time array
When ts is defined as:
ts = np.array([0.0, 0.0022884973945761, 0.0973861319484548, 1.653381907704304, 1.729685002780781])
The integration result shows significant deviation from the expected orbit.
Interestingly, I found two workarounds that produce correct results:
1. ts = np.round(ts, 10)
2. ts = ts + 1
It seems that certain floating-point values in the time array can cause accuracy issues, while rounding or offsetting resolves them.
Thank you for looking into this! Please let me know if you need any additional information.
Best regards
Hi,
First of all, thank you for maintaining gala! I'm using version 1.11.0 and encountered two issues when integrating orbits with specific time arrays.
Code:
mw = gp.MilkyWayPotential(version="latest")
w0 = gd.PhaseSpacePosition(
pos=np.array([-8.122, 0.0, 0.0208]) * u.kpc,
vel=np.array([12.9, 245.6, 7.78]) * u.km / u.s
)
orbit = mw.integrate_orbit(w0, t=ts * u.Myr, Integrator='dop853', cython_if_possible=True)
Issue 1: Integration fails with extremely small time spacing
When ts is defined as:
ts = np.array([1e-3, 1e-3 + 1e-18])
I get the following error:
It seems that when two consecutive time points are extremely close (difference ~1e-18), the integrator fails.
Issue 2: Significant deviation with certain time array
When ts is defined as:
ts = np.array([0.0, 0.0022884973945761, 0.0973861319484548, 1.653381907704304, 1.729685002780781])
The integration result shows significant deviation from the expected orbit.
Interestingly, I found two workarounds that produce correct results:
1. ts = np.round(ts, 10)
2. ts = ts + 1
It seems that certain floating-point values in the time array can cause accuracy issues, while rounding or offsetting resolves them.
Thank you for looking into this! Please let me know if you need any additional information.
Best regards