jobovy / galpy

Galactic Dynamics in python
https://www.galpy.org
BSD 3-Clause "New" or "Revised" License
226 stars 98 forks source link

Wrapper potentials silently fail when applied to DissipativeForces #627

Closed jobovy closed 7 months ago

jobovy commented 8 months ago

🐛 Bug

Wrapper potentials only work for potentials, not general forces (because they inherit from Potential) and when applied to forces, they often silently fail (e.g., in C).

Reproducible example

import numpy
from galpy import potential
from galpy.orbit import Orbit
def M(t):
    return 1.

MW= potential.MWPotential2014
past= numpy.linspace(0.,-4.,1001)*u.Gyr
#Initialize potentials and time-varying potentials
df = potential.ChandrasekharDynamicalFrictionForce(GMs = 3.0*10**10*u.Msun)
df_wrap = potential.TimeDependentAmplitudeWrapperPotential(A=M, amp=1, pot=df)

#Initialize and integrate orbits.
o = Orbit.from_name("Fornax")
oweird = Orbit.from_name("Fornax")
o.integrate(past, pot=MW+df)
oweird.integrate(past, pot=MW+df_wrap,method='dop853')
print(o.x(past[5]), oweird.x(past[5])) # They are not equivalent

(Example adapted from one sent by Urmila Chadayammuri). Note that applying TimeDependentAmplitudeWrapperPotential isn't quite right, because it doesn't consistently change the mass)

Expected behavior

No silent failure, solution could be:

System Details

Emscripten-3.1.14-wasm32-32bit
Python 3.10.2 (main, Aug 29 2022, 18:09:32) [Clang 15.0.0 (https://github.com/llvm/llvm-project 7effcbda49ba32991b8955821b8f
numpy 1.22.4
scipy 1.8.1
matplotlib 3.5.2
galpy 1.9.2.dev0
astropy 5.1

Additional context

Originally reported by Urmila Chadayammuri