Source code for galpy.df.kingdf

# Class that represents a King DF
import numpy
from scipy import special, integrate, interpolate
from ..util import conversion
from .df import df
from .sphericaldf import isotropicsphericaldf

_FOURPI= 4.*numpy.pi
_TWOOVERSQRTPI= 2./numpy.sqrt(numpy.pi)

[docs]class kingdf(isotropicsphericaldf): """Class that represents a King DF"""
[docs] def __init__(self,W0,M=1.,rt=1.,npt=1001,ro=None,vo=None): """ NAME: __init__ PURPOSE: Initialize a King DF INPUT: W0 - dimensionless central potential W0 = Psi(0)/sigma^2 (in practice, needs to be <~ 200, where the DF is essentially isothermal) M= (1.) total mass (can be a Quantity) rt= (1.) tidal radius (can be a Quantity) npt= (1001) number of points to use to solve for Psi(r) ro=, vo= standard galpy unit scaling parameters OUTPUT: (none; sets up instance) HISTORY: 2020-07-09 - Written - Bovy (UofT) """ # Just run df init to set up unit-conversion parameters df.__init__(self,ro=ro,vo=vo) self.W0= W0 self.M= conversion.parse_mass(M,ro=self._ro,vo=self._vo) self.rt= conversion.parse_length(rt,ro=self._ro) # Solve (mass,rtidal)-scale-free model, which is the basis for # the full solution self._scalefree_kdf= _scalefreekingdf(self.W0) self._scalefree_kdf.solve(npt) # Set up scaling factors self._radius_scale= self.rt/self._scalefree_kdf.rt self._mass_scale= self.M/self._scalefree_kdf.mass self._velocity_scale= numpy.sqrt(self._mass_scale/self._radius_scale) self._density_scale= self._mass_scale/self._radius_scale**3. # Store central density, r0... self.rho0= self._scalefree_kdf.rho0*self._density_scale self.r0= self._scalefree_kdf.r0*self._radius_scale self.c= self._scalefree_kdf.c # invariant self.sigma= self._velocity_scale self._sigma2= self.sigma**2. self.rho1= self._density_scale # Setup the potential, use original params in case they had units # because then the initialization will turn on units for this object from ..potential import KingPotential pot= KingPotential(W0=self.W0,M=M,rt=rt,_sfkdf=self._scalefree_kdf, ro=ro,vo=vo) # Now initialize the isotropic DF isotropicsphericaldf.__init__(self,pot=pot,scale=self.r0, rmax=self.rt,ro=ro,vo=vo) self._potInf= self._pot(self.rt,0.,use_physical=False) # Setup inverse cumulative mass function for radius sampling self._icmf= interpolate.InterpolatedUnivariateSpline(\ self._mass_scale*self._scalefree_kdf._cumul_mass/self.M, self._radius_scale*self._scalefree_kdf._r, k=3) # Setup velocity DF interpolator for velocity sampling here self._rmin_sampling= 0. self._v_vesc_pvr_interpolator= self._make_pvr_interpolator(\ r_a_end=numpy.log10(self.rt/self._scale))
def dens(self,r): return self._scalefree_kdf.dens(r/self._radius_scale)\ *self._density_scale def fE(self,E): out= numpy.zeros(numpy.atleast_1d(E).shape) varE= self._potInf-E if numpy.sum(varE > 0.) > 0: out[varE > 0.]= (numpy.exp(varE[varE > 0.]/self._sigma2)-1.)\ *(2.*numpy.pi*self._sigma2)**-1.5*self.rho1 return out# mass density, not /self.M as for number density
class _scalefreekingdf(object): """Internal helper class to solve the scale-free King DF model, that is, the one that only depends on W = Psi/sigma^2""" def __init__(self,W0): self.W0= W0 def solve(self,npt=1001): """Solve the model W(r) at npt points (note: not equally spaced in either r or W, because combination of two ODEs for different r ranges)""" # Set up arrays for outputs r= numpy.zeros(npt) W= numpy.zeros(npt) dWdr= numpy.zeros(npt) # Initialize (r[0]=0 already) W[0]= self.W0 # Determine central density and r0 self.rho0= self._dens_W(self.W0) self.r0= numpy.sqrt(9./4./numpy.pi/self.rho0) # First solve Poisson equation ODE from r=0 to r0 using form # d^2 Psi / dr^2 = ... (d psi / dr = v, r^2 dv / dr = RHS-2*r*v) if self.W0 < 2.: rbreak= self.r0/100. else: rbreak= self.r0 #Using linspace focuses on what happens ~rbreak rather than on <<rbreak # which is what you want, because W ~ constant at r <~ r0 r[:npt//2]= numpy.linspace(0.,rbreak,npt//2) sol= integrate.solve_ivp(\ lambda t,y: [y[1],-_FOURPI*self._dens_W(y[0]) -(2.*y[1]/t if t > 0. else 0.)], [0.,rbreak],[self.W0,0.],method='LSODA',t_eval=r[:npt//2]) W[:npt//2]= sol.y[0] dWdr[:npt//2]= sol.y[1] # Then solve Poisson equation ODE from Psi(r0) to Psi=0 using form # d^2 r / d Psi^2 = ... (d r / d psi = 1/v, dv / dpsi = 1/v(RHS-2*r*v)) # Added advantage that this becomes ~log-spaced in r, which is what # you want W[npt//2-1:]= numpy.linspace(sol.y[0,-1],0.,npt-npt//2+1) sol= integrate.solve_ivp(\ lambda t,y: [1./y[1], -1./y[1]*(_FOURPI*self._dens_W(t) +2.*y[1]/y[0])], [sol.y[0,-1],0.],[rbreak,sol.y[1,-1]], method='LSODA',t_eval=W[npt//2-1:]) r[npt//2-1:]= sol.y[0] dWdr[npt//2-1:]= sol.y[1] # Store solution self._r= r self._W= W self._dWdr= dWdr # Also store density at these points, and the tidal radius self._rho= self._dens_W(self._W) self.rt= r[-1] self.c= numpy.log10(self.rt/self.r0) # Interpolate solution self._W_from_r=\ interpolate.InterpolatedUnivariateSpline(self._r,self._W,k=3) # Compute the cumulative mass and store the total mass mass_shells= numpy.array([\ integrate.quad(lambda r: _FOURPI*r**2*self.dens(r), rlo,rhi)[0] for rlo,rhi in zip(r[:-1],r[1:])]) self._cumul_mass= numpy.hstack((\ integrate.quad(lambda r: _FOURPI*r**2*self.dens(r),0.,r[0])[0], numpy.cumsum(mass_shells))) self.mass= self._cumul_mass[-1] return None def _dens_W(self,W): """Density as a function of W""" sqW= numpy.sqrt(W) return numpy.exp(W)*special.erf(sqW)-_TWOOVERSQRTPI*sqW*(1.+2./3.*W) def dens(self,r): return self._dens_W(self._W_from_r(r))