Merge 64845f4 into 9ed100e

galpy/df/__init__.py

@@ -5,6 +5,10 @@
 from . import streamdf
 from . import streamgapdf
 from . import jeans
+from . import eddingtondf
+from . import isotropicHernquistdf
+from . import constantbetaHernquistdf
+from . import kingdf
 #
 # Functions
 #
@@ -32,3 +36,7 @@
 quasiisothermaldf= quasiisothermaldf.quasiisothermaldf
 streamdf= streamdf.streamdf
 streamgapdf= streamgapdf.streamgapdf
+eddingtondf= eddingtondf.eddingtondf
+isotropicHernquistdf= isotropicHernquistdf.isotropicHernquistdf
+constantbetaHernquistdf= constantbetaHernquistdf.constantbetaHernquistdf
+kingdf= kingdf.kingdf

galpy/df/constantbetaHernquistdf.py

@@ -0,0 +1,99 @@
+# Class that implements the anisotropic spherical Hernquist DF with constant
+# beta parameter
+import numpy
+import scipy.special
+import scipy.integrate
+from ..util import conversion
+from ..potential import evaluatePotentials,HernquistPotential
+from .constantbetadf import constantbetadf
+
+class constantbetaHernquistdf(constantbetadf):
+    """Class that implements the anisotropic spherical Hernquist DF with constant beta parameter"""
+    def __init__(self,pot=None,beta=0,ro=None,vo=None):
+        """
+        NAME:
+
+            __init__
+
+        PURPOSE:
+
+            Initialize a Hernquist DF with constant anisotropy
+
+        INPUT:
+
+            pot - Hernquist potential which determines the DF
+
+            beta - anisotropy parameter
+
+        OUTPUT:
+
+            None
+
+        HISTORY:
+
+            2020-07-22 - Written - Lane (UofT)
+        """
+        assert isinstance(pot,HernquistPotential),'pot= must be potential.HernquistPotential'
+        constantbetadf.__init__(self,pot=pot,beta=beta,ro=ro,vo=vo)
+        self._psi0= -evaluatePotentials(self._pot,0,0,use_physical=False)
+        self._GMa = self._psi0*self._pot.a**2.
+        self._fEnorm= (2.**self._beta/(2.*numpy.pi)**2.5)\
+            *scipy.special.gamma(5.-2.*self._beta)\
+            /scipy.special.gamma(1.-self._beta)\
+            /scipy.special.gamma(3.5-self._beta)\
+            /self._GMa**(1.5-self._beta)
+
+    def fE(self,E):
+        """
+        NAME:
+
+            fE
+
+        PURPOSE
+
+            Calculate the energy portion of a Hernquist distribution function
+
+        INPUT:
+
+            E - The energy (can be Quantity)
+
+        OUTPUT:
+
+            fE - The value of the energy portion of the DF
+
+        HISTORY:
+
+            2020-07-22 - Written
+        """
+        Etilde= -conversion.parse_energy(E,vo=self._vo)/self._psi0
+        # Handle potential E outside of bounds
+        Etilde_out = numpy.where(numpy.logical_or(Etilde<0,Etilde>1))[0]
+        if len(Etilde_out)>0:
+            # Dummy variable now and 0 later, prevents numerical issues?
+            Etilde[Etilde_out]=0.5
+        # First check algebraic solutions
+        if self._beta == 0.: # isotropic case
+            sqrtEtilde= numpy.sqrt(Etilde)
+            fE= 1./numpy.sqrt(2.)/(2*numpy.pi)**3/self._GMa**1.5\
+                *sqrtEtilde/(1-Etilde)**2.\
+                *((1.-2.*Etilde)*(8.*Etilde**2.-8.*Etilde-3.)\
+                  +((3.*numpy.arcsin(sqrtEtilde))\
+                    /numpy.sqrt(Etilde*(1.-Etilde))))
+        elif self._beta == 0.5:
+            fE= (3.*Etilde**2.)/(4.*numpy.pi**3.*self._GMa)
+        elif self._beta == -0.5:
+            fE= ((20.*Etilde**3.-20.*Etilde**4.+6.*Etilde**5.)\
+                 /(1.-Etilde)**4)/(4.*numpy.pi**3.*self._GMa**2.)
+        else:
+            fE= self._fEnorm*numpy.power(Etilde,2.5-self._beta)*\
+                scipy.special.hyp2f1(5.-2.*self._beta,1.-2.*self._beta,
+                                     3.5-self._beta,Etilde)
+        if len(Etilde_out) > 0:
+            fE[Etilde_out]= 0.
+        return fE
+
+
+    def _icmf(self,ms):
+        '''Analytic expression for the normalized inverse cumulative mass 
+        function. The argument ms is normalized mass fraction [0,1]'''
+        return self._pot.a*numpy.sqrt(ms)/(1-numpy.sqrt(ms))

galpy/df/constantbetadf.py

@@ -0,0 +1,86 @@
+# Class that implements DFs of the form f(E,L) = L^{-2\beta} f(E) with constant
+# beta anisotropy parameter
+import numpy
+import scipy.interpolate
+from scipy import integrate, special
+from ..potential import evaluatePotentials
+from .sphericaldf import anisotropicsphericaldf
+
+class constantbetadf(anisotropicsphericaldf):
+    """Class that implements DFs of the form f(E,L) = L^{-2\beta} f(E) with constant beta anisotropy parameter"""
+    def __init__(self,pot=None,beta=None,scale=None,ro=None,vo=None):
+        """
+        NAME:
+
+            __init__
+
+        PURPOSE:
+
+            Initialize a spherical DF with constant anisotropy parameter
+
+        INPUT:
+
+            pot - Spherical potential which determines the DF
+
+            scale - Characteristic scale radius to aid sampling calculations. 
+                Not necessary, and will also be overridden by value from pot if 
+                available.
+
+        """
+        anisotropicsphericaldf.__init__(self,pot=pot,scale=scale,ro=ro,vo=vo)
+        self._beta= beta
+        self._potInf= evaluatePotentials(pot,10**12,0)
+
+    def _call_internal(self,*args):
+        """
+        NAME:
+
+            _call_internal
+
+        PURPOSE:
+
+            Evaluate the DF for a constant anisotropy Hernquist
+
+        INPUT:
+
+            E - The energy
+
+            L - The angular momentum
+
+        OUTPUT:
+
+            fH - The value of the DF
+
+        HISTORY:
+
+            2020-07-22 - Written - Lane (UofT)
+        """
+        E, L, _= args
+        return L**(-2*self._beta)*self.fE(E)
+
+    def _sample_eta(self,n=1):
+        """Sample the angle eta which defines radial vs tangential velocities"""
+        deta = 0.00005*numpy.pi
+        etas = (numpy.arange(0, numpy.pi, deta)+deta/2)
+        eta_pdf_cml = numpy.cumsum(numpy.power(numpy.sin(etas),1.-2.*self._beta))
+        eta_pdf_cml_norm = eta_pdf_cml / eta_pdf_cml[-1]
+        eta_icml_interp = scipy.interpolate.interp1d(eta_pdf_cml_norm, etas, 
+            bounds_error=False, fill_value='extrapolate')
+        eta_samples = eta_icml_interp(numpy.random.uniform(size=n))
+        return eta_samples
+
+    def _p_v_at_r(self,v,r):
+        return self.fE(evaluatePotentials(self._pot,r,0,use_physical=False)\
+                       +0.5*v**2.)*v**(2.-2.*self._beta)
+
+    def _vmomentdensity(self,r,n,m):
+         if m%2 == 1 or n%2 == 1:
+             return 0.
+         return 2.*numpy.pi*r**(-2.*self._beta)\
+             *integrate.quad(lambda v: v**(2.-2.*self._beta+m+n)
+                             *self.fE(evaluatePotentials(self._pot,r,0,
+                                                         use_physical=False)
+                                      +0.5*v**2.),
+                             0.,self._vmax_at_r(self._pot,r))[0]\
+            *special.gamma(m/2.-self._beta+1.)*special.gamma((n+1)/2.)/\
+            special.gamma(0.5*(m+n-2.*self._beta+3.))

galpy/df/eddingtondf.py

@@ -0,0 +1,15 @@
+# Class that implements isotropic spherical DFs computed using the Eddington
+# formula
+from ..potential import evaluatePotentials
+from .sphericaldf import isotropicsphericaldf
+
+class eddingtondf(isotropicsphericaldf):
+    """Class that implements isotropic spherical DFs computed using the Eddington formula"""
+    def __init__(self,pot=None,scale=None,ro=None,vo=None):
+        """
+            scale - Characteristic scale radius to aid sampling calculations. 
+                Not necessary, and will also be overridden by value from pot if 
+                available.
+        """
+        isotropicsphericaldf.__init__(self,pot=pot,scale=scale,ro=ro,vo=vo)
+        self._potInf= evaluatePotentials(pot,10**12,0)

galpy/df/isotropicHernquistdf.py

@@ -0,0 +1,58 @@
+# Class that implements isotropic spherical Hernquist DF
+# computed using the Eddington formula
+import numpy
+from ..util import conversion
+from ..potential import evaluatePotentials,HernquistPotential
+from .eddingtondf import eddingtondf
+
+class isotropicHernquistdf(eddingtondf):
+    """Class that implements isotropic spherical Hernquist DF computed using the Eddington formula"""
+    def __init__(self,pot=None,ro=None,vo=None):
+        assert isinstance(pot,HernquistPotential),'pot= must be potential.HernquistPotential'
+        eddingtondf.__init__(self,pot=pot,ro=ro,vo=vo)
+        self._psi0= -evaluatePotentials(self._pot,0,0,use_physical=False)
+        self._GMa = self._psi0*self._pot.a**2.
+        self._fEnorm= 1./numpy.sqrt(2.)/(2*numpy.pi)**3/self._GMa**1.5
+
+    def fE(self,E):
+        """
+        NAME:
+
+            fE
+
+        PURPOSE
+
+            Calculate the energy portion of an isotropic Hernquist distribution function
+
+        INPUT:
+
+            E - The energy (can be Quantity)
+
+        OUTPUT:
+
+            fE - The value of the energy portion of the DF
+
+        HISTORY:
+
+            2020-08-09 - Written - James Lane (UofT)
+        """
+        Etilde= -conversion.parse_energy(E,vo=self._vo)/self._psi0
+        # Handle E out of bounds
+        Etilde_out = numpy.where(numpy.logical_or(Etilde<0,Etilde>1))[0]
+        if len(Etilde_out)>0:
+            # Set to dummy and 0 later, prevents functions throwing errors?
+            Etilde[Etilde_out]=0.5
+        sqrtEtilde= numpy.sqrt(Etilde)
+        fE= self._fEnorm*sqrtEtilde/(1-Etilde)**2.\
+            *((1.-2.*Etilde)*(8.*Etilde**2.-8.*Etilde-3.)\
+              +((3.*numpy.arcsin(sqrtEtilde))\
+                /numpy.sqrt(Etilde*(1.-Etilde))))
+        # Fix out of bounds values
+        if len(Etilde_out) > 0:
+            fE[Etilde_out] = 0
+        return fE
+
+    def _icmf(self,ms):
+        '''Analytic expression for the normalized inverse cumulative mass 
+        function. The argument ms is normalized mass fraction [0,1]'''
+        return self._pot.a*numpy.sqrt(ms)/(1-numpy.sqrt(ms))