Questions on the issue with time-segmented integration in galpy.

[MandLforever]

Dear Jo and the community,
I am trying to integrate the Sgr orbit backwards by time-segmented integration, considering the dynamical friction. Surprisingly, if the integration is conducted by intervals, the generated orbits are not continuous and have weird shape. A piece of the code is below:

M_sgr1 = 4e8units.Msun
M_sgr2 = 4e9units.Msun
M_sgr3 = 4e10*units.Msun
o = copy.deepcopy(o_sgr)

cdf1= ChandrasekharDynamicalFrictionForce(GMs= M_sgr1, rhm = rhm,dens=mw14,sigmar = sigmar, gamma = np.sqrt(2))
cdf2= ChandrasekharDynamicalFrictionForce(GMs= M_sgr2, rhm = rhm*(40)(1/3),dens=mw14,sigmar = sigmar, gamma = np.sqrt(2))
cdf3= ChandrasekharDynamicalFrictionForce(GMs= M_sgr3, rhm = rhm*(400)(1/3),dens=mw14,sigmar = sigmar, gamma = np.sqrt(2))

t = np.linspace(0,-6,6001)*units.Gyr

t1 = t[:3000]
t2 = t[3000:]

o.integrate(t1,mw14 + cdf2)
r_2 = list(o.r(t1))
o_ = copy.deepcopy(o(t1[-1]))
o_.integrate(t2,mw14 + cdf3) ##It is also wrong if set cdf3 to cdf2
r_2 = list(o_.r(t2)) + r_2

plt.plot(t,r_2)
plt.xlabel("t (Gyr)")
plt.ylabel("r (kpc)")

And according to my limited understandings, with larger dynamical friction, the orbit should have larger apocenter but the plot also did not give this feature .

And I am also wondering how I should deal with mass-(time)dependent dynamical friction when this method is wrong.

[jobovy]

Hi!

I think the issue is simply with how you get the final r_2 array. Doing

r_2 = list(o_.r(t2)) + r_2

puts the second orbit integration first (-3 Gyr to -6 Gyr) and then attaches the first orbit integration (0 to -3 Gyr). You can see that the first and the last point in your plot line up, because that's where they should connect. So I think you simply need to do

r_2 = r2 + list(o_.r(t2))

I should also point out that in the main branch of galpy, it is now possible to continue an orbit integration (see #781), so you wouldn't need to instantiate o_ and would get the full r_2 array in the end. To install this version, you can do

python -m pip install --pre --extra-index-url https://www.galpy.org/wheelhouse/simple --only-binary galpy galpy

and then you can do

o.integrate(t1,mw14 + cdf2)
o.integrate(t2,mw14 + cdf3)
r_2 = o(t)

or at the end just do

o.plotr()