# Source code for galpy.potential.NumericalPotentialDerivativesMixin

###############################################################################
#   NumericalPotentialDerivativesMixin: helper class to add numerical derivs
###############################################################################
[docs]class NumericalPotentialDerivativesMixin(object):
"""Mixin to add numerical derivatives to a Potential class, use as, e.g.,

.. highlight:: python
.. code-block:: python

class PotWithNumericalDerivs(Potential,NumericalPotentialDerivativesMixin):
def __init__(self,*args,**kwargs):
NumericalPotentialDerivativesMixin.__init__(self,kwargs) # *not* **kwargs!
# Remainder of initialization
...

def _evaluate(self,R,z,phi=0.,t=0.):
# Evaluate the potential

# All forces and second derivatives then computed by NumericalPotentialDerivativesMixin

to add numerical derivatives to a new potential class PotWithNumericalDerivs that only implements the potential itself, but not the forces. The class may implement any of the forces or second derivatives, all non-implemented forces/second-derivatives will be computed numerically by adding this Mixin

The step used to compute the first (force) and second derivatives can be controlled at object instantiation by the keyword arguments dR, dz, dphi (for the forces; 1e-8 default) and dR2, dz2, and dphi2 (for the second derivaives; 1e-4 default)"""
[docs]    def __init__(self,kwargs): # no **kwargs to get a reference, not a copy!
# For first derivatives
self._dR= kwargs.pop('dR',1e-8)
self._dphi= kwargs.pop('dphi',1e-8)
self._dz= kwargs.pop('dz',1e-8)
# For second derivatives
self._dR2= kwargs.pop('dR2',1e-4)
self._dphi2= kwargs.pop('dphi2',1e-4)
self._dz2= kwargs.pop('dz2',1e-4)

def _Rforce(self,R,z,phi=0.,t=0.):
# Do forward difference because R cannot be negative
RplusdR= R+self._dR
Rplus2dR= R+2.*self._dR
dR= (Rplus2dR-R)/2.
return (1.5*self._evaluate(R,z,phi=phi,t=t)
-2.*self._evaluate(RplusdR,z,phi=phi,t=t)
+0.5*self._evaluate(Rplus2dR,z,phi=phi,t=t))/dR

def _zforce(self,R,z,phi=0.,t=0.):
# Central difference to get derivative at z=0 right
zplusdz= z+self._dz
zminusdz= z-self._dz
dz= zplusdz-zminusdz
return (self._evaluate(R,zminusdz,phi=phi,t=t)
-self._evaluate(R,zplusdz,phi=phi,t=t))/dz

def _phiforce(self,R,z,phi=0.,t=0.):
if not self.isNonAxi: return 0.
# Central difference
phiplusdphi= phi+self._dphi
phiminusdphi= phi-self._dphi
dphi= phiplusdphi-phiminusdphi
return (self._evaluate(R,z,phi=phiminusdphi,t=t)
-self._evaluate(R,z,phi=phiplusdphi,t=t))/dphi

def _R2deriv(self,R,z,phi=0.,t=0.):
# Do forward difference because R cannot be negative
RplusdR= R+self._dR2
Rplus2dR= R+2.*self._dR2
Rplus3dR= R+3.*self._dR2
dR= (Rplus3dR-R)/3.
return (2.*self._evaluate(R,z,phi=phi,t=t)
-5.*self._evaluate(RplusdR,z,phi=phi,t=t)
+4.*self._evaluate(Rplus2dR,z,phi=phi,t=t)
-1.*self._evaluate(Rplus3dR,z,phi=phi,t=t))/dR**2.

def _z2deriv(self,R,z,phi=0.,t=0.):
# Central derivative
zplusdz= z+self._dz2
zminusdz= z-self._dz2
dz= (zplusdz-zminusdz)/2.
return (self._evaluate(R,zplusdz,phi=phi,t=t)
+self._evaluate(R,zminusdz,phi=phi,t=t)
-2.*self._evaluate(R,z,phi=phi,t=t))/dz**2.

def _phi2deriv(self,R,z,phi=0.,t=0.):
if not self.isNonAxi: return 0.
# Central derivative
phiplusdphi= phi+self._dphi2
phiminusdphi= phi-self._dphi2
dphi= (phiplusdphi-phiminusdphi)/2.
return (self._evaluate(R,z,phi=phiplusdphi,t=t)
+self._evaluate(R,z,phi=phiminusdphi,t=t)
-2.*self._evaluate(R,z,phi=phi,t=t))/dphi**2.

def _Rzderiv(self,R,z,phi=0.,t=0.):
# Do forward difference in R because R cannot be negative
RplusdR= R+self._dR2
Rplus2dR= R+2.*self._dR2
dR= (Rplus2dR-R)/2.
zplusdz= z+self._dz2
zminusdz= z-self._dz2
dz= zplusdz-zminusdz
return (-1.5*self._evaluate(R,zplusdz,phi=phi,t=t)
+2.*self._evaluate(RplusdR,zplusdz,phi=phi,t=t)
-0.5*self._evaluate(Rplus2dR,zplusdz,phi=phi,t=t)
+1.5*self._evaluate(R,zminusdz,phi=phi,t=t)
-2.*self._evaluate(RplusdR,zminusdz,phi=phi,t=t)
+0.5*self._evaluate(Rplus2dR,zminusdz,phi=phi,t=t))/dR/dz

def _Rphideriv(self,R,z,phi=0.,t=0.):
if not self.isNonAxi: return 0.
# Do forward difference in R because R cannot be negative
RplusdR= R+self._dR2
Rplus2dR= R+2.*self._dR2
dR= (Rplus2dR-R)/2.
phiplusdphi= phi+self._dphi2
phiminusdphi= phi-self._dphi2
dphi= phiplusdphi-phiminusdphi
return (-1.5*self._evaluate(R,z,phi=phiplusdphi,t=t)
+2.*self._evaluate(RplusdR,z,phi=phiplusdphi,t=t)
-0.5*self._evaluate(Rplus2dR,z,phi=phiplusdphi,t=t)
+1.5*self._evaluate(R,z,phi=phiminusdphi,t=t)
-2.*self._evaluate(RplusdR,z,phi=phiminusdphi,t=t)
+0.5*self._evaluate(Rplus2dR,z,phi=phiminusdphi,t=t))\
/dR/dphi

def _phizderiv(self,R,z,phi=0.,t=0.):
if not self.isNonAxi: return 0.
# Central derivative
phiplusdphi= phi+self._dphi2
phiminusdphi= phi-self._dphi2
dphi= (phiplusdphi-phiminusdphi)/2.
zplusdz= z+self._dz2
zminusdz= z-self._dz2
dz= zplusdz-zminusdz
return (self._evaluate(R,zplusdz,phi=phiplusdphi,t=t)
-self._evaluate(R,zplusdz,phi=phiminusdphi,t=t)
-self._evaluate(R,zminusdz,phi=phiplusdphi,t=t)
+self._evaluate(R,zminusdz,phi=phiminusdphi,t=t))\
/dz/dphi/2.