Source code for galpy.df.quasiisothermaldf

#A 'Binney' quasi-isothermal DF
import warnings
import hashlib
import numpy
from scipy import optimize, interpolate, integrate
from .. import potential
from .. import actionAngle
from ..actionAngle import actionAngleIsochrone
from ..potential import IsochronePotential
from ..potential import flatten as flatten_potential
from ..orbit import Orbit
from .df import df
from ..util import galpyWarning
from ..util.conversion import physical_conversion, \
    potential_physical_input, actionAngle_physical_input, _APY_UNITS, \
    physical_compatible, parse_length, parse_velocity, parse_angmom, \
    parse_length_kpc, parse_velocity_kms, _APY_LOADED
if _APY_LOADED:
    from astropy import units
_NSIGMA=4
_DEFAULTNGL=10
_DEFAULTNGL2=20
[docs]class quasiisothermaldf(df): """Class that represents a 'Binney' quasi-isothermal DF"""
[docs] def __init__(self,hr,sr,sz,hsr,hsz,pot=None,aA=None, cutcounter=False, _precomputerg=True,_precomputergrmax=None, _precomputergnLz=51, refr=1.,lo=10./220./8., ro=None,vo=None): """ NAME: __init__ PURPOSE: Initialize a quasi-isothermal DF INPUT: hr - radial scale length (can be Quantity) sr - radial velocity dispersion at the solar radius (can be Quantity) sz - vertical velocity dispersion at the solar radius (can be Quantity) hsr - radial-velocity-dispersion scale length (can be Quantity) hsz - vertial-velocity-dispersion scale length (can be Quantity) pot= Potential instance or list thereof aA= actionAngle instance used to convert (x,v) to actions [must be an instance of an actionAngle class that computes (J,Omega,angle) for a given (x,v)] cutcounter= if True, set counter-rotating stars' DF to zero refr= reference radius for dispersions (can be different from ro) (can be Quantity) lo= reference angular momentum below where there are significant numbers of retrograde stars (can be Quantity) ro= distance from vantage point to GC (kpc; can be Quantity) vo= circular velocity at ro (km/s; can be Quantity) OTHER INPUTS: _precomputerg= if True (default), pre-compute the rL(L) _precomputergrmax= if set, this is the maximum R for which to pre-compute rg (default: 5*hr) _precomputergnLz if set, number of Lz to pre-compute rg for (default: 51) OUTPUT: object HISTORY: 2012-07-25 - Started - Bovy (IAS@MPIA) """ df.__init__(self,ro=ro,vo=vo) self._hr= parse_length(hr,ro=self._ro) self._sr= parse_velocity(sr,vo=self._vo) self._sz= parse_velocity(sz,vo=self._vo) self._hsr= parse_length(hsr,ro=self._ro) self._hsz= parse_length(hsz,ro=self._ro) self._refr= parse_length(refr,ro=self._ro) self._lo= parse_angmom(lo,ro=self._ro,vo=self._vo) self._lnsr= numpy.log(self._sr) self._lnsz= numpy.log(self._sz) self._maxVT_hash= None self._maxVT_ip= None if pot is None: raise IOError("pot= must be set") self._pot= flatten_potential(pot) if aA is None: raise IOError("aA= must be set") self._aA= aA if not self._aA._pot == self._pot: if not isinstance(self._aA,actionAngleIsochrone): raise IOError("Potential in aA does not appear to be the same as given potential pot") elif isinstance(self._pot,IsochronePotential) and \ not self._aA.b == self._pot.b and \ not self._aA.amp == self._pot._amp: raise IOError("Potential in aA does not appear to be the same as given potential pot") self._check_consistent_units() self._cutcounter= cutcounter if _precomputerg: if _precomputergrmax is None: _precomputergrmax= 5*self._hr self._precomputergrmax= _precomputergrmax self._precomputergnLz= _precomputergnLz self._precomputergLzmin= 0.01 self._precomputergLzmax= self._precomputergrmax\ *potential.vcirc(self._pot,self._precomputergrmax) self._precomputergLzgrid= numpy.linspace(self._precomputergLzmin,self._precomputergLzmax,self._precomputergnLz) self._rls= numpy.array([potential.rl(self._pot,l) for l in self._precomputergLzgrid]) #Spline interpolate self._rgInterp= interpolate.InterpolatedUnivariateSpline(self._precomputergLzgrid,self._rls,k=3) else: self._precomputergrmax= 0. self._rgInterp= None self._rls= None self._precomputergnr= None self._precomputergLzgrid= None self._precomputergLzmin= \ numpy.finfo(numpy.dtype(numpy.float64)).max self._precomputergLzmax= \ numpy.finfo(numpy.dtype(numpy.float64)).min self._precomputerg= _precomputerg self._glxdef, self._glwdef= \ numpy.polynomial.legendre.leggauss(_DEFAULTNGL) self._glxdef2, self._glwdef2= \ numpy.polynomial.legendre.leggauss(_DEFAULTNGL2) self._glxdef12, self._glwdef12= \ numpy.polynomial.legendre.leggauss(_DEFAULTNGL//2) return None
[docs] @physical_conversion('phasespacedensity',pop=True) def __call__(self,*args,**kwargs): """ NAME: __call__ PURPOSE: return the DF INPUT: Either: a)(jr,lz,jz) tuple; each can be a Quantity where: jr - radial action lz - z-component of angular momentum jz - vertical action b) R,vR,vT,z,vz c) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well log= if True, return the natural log +scipy.integrate.quadrature kwargs func= function of (jr,lz,jz) to multiply f with (useful for moments) OUTPUT: value of DF HISTORY: 2012-07-25 - Written - Bovy (IAS@MPIA) NOTE: For Miyamoto-Nagai/adiabatic approximation this seems to take about 30 ms / evaluation in the extended Solar neighborhood For a MWPotential/adiabatic approximation this takes about 50 ms / evaluation in the extended Solar neighborhood For adiabatic-approximation grid this seems to take about 0.67 to 0.75 ms / evaluation in the extended Solar neighborhood (includes some out of the grid) up to 200x faster when called with vector R,vR,vT,z,vz """ #First parse log log= kwargs.pop('log',False) _return_actions= kwargs.pop('_return_actions',False) _return_freqs= kwargs.pop('_return_freqs',False) _func= kwargs.pop('func',None) if 'rg' in kwargs: thisrg= kwargs.pop('rg') kappa= kwargs.pop('kappa') nu= kwargs.pop('nu') Omega= kwargs.pop('Omega') else: thisrg= None kappa= None nu= None Omega= None #First parse args if len(args) == 1 and not isinstance(args[0],Orbit): #(jr,lz,jz) jr,lz,jz= args[0] jr= parse_angmom(jr,ro=self._ro,vo=self._vo) lz= parse_angmom(lz,ro=self._ro,vo=self._vo) jz= parse_angmom(jz,ro=self._ro,vo=self._vo) else: #Use self._aA to calculate the actions if isinstance(args[0],Orbit) and len(args[0].shape) > 1: raise RuntimeError("Evaluating quasiisothermaldf with Orbit instances with multi-dimensional shapes is not supported") #pragma: no cover try: jr,lz,jz= self._aA(*args,use_physical=False,**kwargs) except actionAngle.UnboundError: if log: return -numpy.finfo(numpy.dtype(numpy.float64)).max else: return 0. #if isinstance(jr,(list,numpy.ndarray)) and len(jr) > 1: jr= jr[0] #if isinstance(jz,(list,numpy.ndarray)) and len(jz) > 1: jz= jz[0] if not isinstance(lz,numpy.ndarray) and self._cutcounter and lz < 0.: if log: return -numpy.finfo(numpy.dtype(numpy.float64)).max else: return 0. #First calculate rg if thisrg is None: thisrg= self._rg(lz) #Then calculate the epicycle and vertical frequencies kappa, nu= self._calc_epifreq(thisrg), self._calc_verticalfreq(thisrg) Omega= numpy.fabs(lz)/thisrg/thisrg #calculate surface-densities and sigmas lnsurfmass= (self._refr-thisrg)/self._hr lnsr= self._lnsr+(self._refr-thisrg)/self._hsr lnsz= self._lnsz+(self._refr-thisrg)/self._hsz #Calculate func if not _func is None: if log: funcTerm= numpy.log(_func(jr,lz,jz)) else: funcFactor= _func(jr,lz,jz) #Calculate fsr else: if log: funcTerm= 0. else: funcFactor= 1. if log: lnfsr= numpy.log(Omega)+lnsurfmass-2.*lnsr-numpy.log(numpy.pi)\ -numpy.log(kappa)\ +numpy.log(1.+numpy.tanh(lz/self._lo))\ -kappa*jr*numpy.exp(-2.*lnsr) lnfsz= numpy.log(nu)-numpy.log(2.*numpy.pi)\ -2.*lnsz-nu*jz*numpy.exp(-2.*lnsz) out= lnfsr+lnfsz+funcTerm if isinstance(lz,numpy.ndarray): out[numpy.isnan(out)]= -numpy.finfo(numpy.dtype(numpy.float64)).max if self._cutcounter: out[(lz < 0.)]= -numpy.finfo(numpy.dtype(numpy.float64)).max elif numpy.isnan(out): out= -numpy.finfo(numpy.dtype(numpy.float64)).max else: srm2= numpy.exp(-2.*lnsr) fsr= Omega*numpy.exp(lnsurfmass)*srm2/numpy.pi/kappa\ *(1.+numpy.tanh(lz/self._lo))\ *numpy.exp(-kappa*jr*srm2) szm2= numpy.exp(-2.*lnsz) fsz= nu/2./numpy.pi*szm2*numpy.exp(-nu*jz*szm2) out= fsr*fsz*funcFactor if isinstance(lz,numpy.ndarray): out[numpy.isnan(out)]= 0. if self._cutcounter: out[(lz < 0.)]= 0. elif numpy.isnan(out): out= 0. if _return_actions and _return_freqs: return (out,jr,lz,jz,thisrg,kappa,nu,Omega) elif _return_actions: return (out,jr,lz,jz) elif _return_freqs: return (out,thisrg,kappa,nu,Omega) else: return out
[docs] @potential_physical_input @physical_conversion('position',pop=True) def estimate_hr(self,R,z=0.,dR=10.**-8.,**kwargs): """ NAME: estimate_hr PURPOSE: estimate the exponential scale length at R INPUT: R - Galactocentric radius (can be Quantity) z= height (default: 0 pc) (can be Quantity) dR- range in R to use (can be Quantity) density kwargs OUTPUT: estimated hR HISTORY: 2012-09-11 - Written - Bovy (IAS) 2013-01-28 - Re-written - Bovy """ Rs= [R-dR/2.,R+dR/2.] if z is None: sf= numpy.array([self.surfacemass_z(r,use_physical=False, **kwargs) for r in Rs]) else: sf= numpy.array([self.density(r,z,use_physical=False, **kwargs) for r in Rs]) lsf= numpy.log(sf) return -dR/(lsf[1]-lsf[0])
[docs] @potential_physical_input @physical_conversion('position',pop=True) def estimate_hz(self,R,z,dz=10.**-8.,**kwargs): """ NAME: estimate_hz PURPOSE: estimate the exponential scale height at R INPUT: R - Galactocentric radius (can be Quantity) dz - z range to use (can be Quantity) density kwargs OUTPUT: estimated hz HISTORY: 2012-08-30 - Written - Bovy (IAS) 2013-01-28 - Re-written - Bovy """ if z == 0.: zs= [z,z+dz] else: zs= [z-dz/2.,z+dz/2.] sf= numpy.array([self.density(R,zz,use_physical=False, **kwargs) for zz in zs]) lsf= numpy.log(sf) return -dz/(lsf[1]-lsf[0])
[docs] @potential_physical_input @physical_conversion('position',pop=True) def estimate_hsr(self,R,z=0.,dR=10.**-8.,**kwargs): """ NAME: estimate_hsr PURPOSE: estimate the exponential scale length of the radial dispersion at R INPUT: R - Galactocentric radius (can be Quantity) z= height (default: 0 pc) (can be Quantity) dR- range in R to use (can be Quantity) density kwargs OUTPUT: estimated hsR HISTORY: 2013-03-08 - Written - Bovy (IAS) """ Rs= [R-dR/2.,R+dR/2.] sf= numpy.array([self.sigmaR2(r,z,use_physical=False, **kwargs) for r in Rs]) lsf= numpy.log(sf)/2. return -dR/(lsf[1]-lsf[0])
[docs] @potential_physical_input @physical_conversion('position',pop=True) def estimate_hsz(self,R,z=0.,dR=10.**-8.,**kwargs): """ NAME: estimate_hsz PURPOSE: estimate the exponential scale length of the vertical dispersion at R INPUT: R - Galactocentric radius (can be Quantity) z= height (default: 0 pc) (can be Quantity) dR- range in R to use (can be Quantity) density kwargs OUTPUT: estimated hsz HISTORY: 2013-03-08 - Written - Bovy (IAS) """ Rs= [R-dR/2.,R+dR/2.] sf= numpy.array([self.sigmaz2(r,z,use_physical=False, **kwargs) for r in Rs]) lsf= numpy.log(sf)/2. return -dR/(lsf[1]-lsf[0])
[docs] @potential_physical_input @physical_conversion('numbersurfacedensity',pop=True) def surfacemass_z(self,R,nz=7,zmax=1.,fixed_quad=True,fixed_order=8, **kwargs): """ NAME: surfacemass_z PURPOSE: calculate the vertically-integrated surface density INPUT: R - Galactocentric radius (can be Quantity) fixed_quad= if True (default), use Gauss-Legendre integration fixed_order= (20), order of GL integration to use nz= number of zs to use to estimate zmax= maximum z to use (can be Quantity) density kwargs OUTPUT: \Sigma(R) HISTORY: 2012-08-30 - Written - Bovy (IAS) """ if fixed_quad: return 2.*integrate.fixed_quad(lambda x: self.density(R*numpy.ones(fixed_order),x,use_physical=False), 0.,.5,n=fixed_order)[0] zs= numpy.linspace(0.,zmax,nz) sf= numpy.array([self.density(R,z,use_physical=False, **kwargs) for z in zs]) lsf= numpy.log(sf) #Interpolate lsfInterp= interpolate.UnivariateSpline(zs, lsf, k=3) #Integrate return 2.*integrate.quad((lambda x: numpy.exp(lsfInterp(x))), 0.,1.)[0]
def vmomentdensity(self,*args,**kwargs): """ NAME: vmomentdensity PURPOSE: calculate the an arbitrary moment of the velocity distribution at R times the density INPUT: R - radius at which to calculate the moment(/ro) n - vR^n m - vT^m o - vz^o OPTIONAL INPUT: nsigma - number of sigma to integrate the vR and vz velocities over (when doing explicit numerical integral; default: 4) vTmax - upper limit for integration over vT (default: 1.5) mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= use Gauss-Legendre _returngl= if True, return the evaluated DF _return_actions= if True, return the evaluated actions (does not work with _returngl currently) _return_freqs= if True, return the evaluated frequencies and rg (does not work with _returngl currently) OUTPUT: <vR^n vT^m x density> at R,z (no support for units) HISTORY: 2012-08-06 - Written - Bovy (IAS@MPIA) """ use_physical= kwargs.pop('use_physical',True) ro= kwargs.pop('ro',None) if ro is None and hasattr(self,'_roSet') and self._roSet: ro= self._ro ro= parse_length_kpc(ro) vo= kwargs.pop('vo',None) if vo is None and hasattr(self,'_voSet') and self._voSet: vo= self._vo vo= parse_velocity_kms(vo) if use_physical and not vo is None and not ro is None: fac= vo**(args[2]+args[3]+args[4])/ro**3 if _APY_UNITS: u= 1/units.kpc**3*(units.km/units.s)**(args[2]+args[3]+args[4]) out= self._vmomentdensity(*args,**kwargs) if _APY_UNITS: return units.Quantity(out*fac,unit=u) else: return out*fac else: return self._vmomentdensity(*args,**kwargs)
[docs] def _vmomentdensity(self,R,z,n,m,o,nsigma=None,mc=False,nmc=10000, _returnmc=False,_vrs=None,_vts=None,_vzs=None, _rawgausssamples=False, gl=False,ngl=_DEFAULTNGL,_returngl=False,_glqeval=None, _return_actions=False,_jr=None,_lz=None,_jz=None, _return_freqs=False, _rg=None,_kappa=None,_nu=None,_Omega=None, _sigmaR1=None,_sigmaz1=None, **kwargs): """Non-physical version of vmomentdensity, otherwise the same""" if isinstance(R,numpy.ndarray): return numpy.array([self._vmomentdensity(r,zz,n,m,o,nsigma=nsigma, mc=mc,nmc=nmc, gl=gl,ngl=ngl,**kwargs) for r,zz in zip(R,z)]) if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): if n % 2 == 1. or o % 2 == 1.: return 0. #we know this must be the case if nsigma == None: nsigma= _NSIGMA if _sigmaR1 is None: sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr) else: sigmaR1= _sigmaR1 if _sigmaz1 is None: sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz) else: sigmaz1= _sigmaz1 thisvc= potential.vcirc(self._pot,R,use_physical=False) #Use the asymmetric drift equation to estimate va gamma= numpy.sqrt(0.5) va= sigmaR1**2./2./thisvc\ *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5 +R*(1./self._hr+2./self._hsr)) if numpy.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") if not _glqeval is None and ngl != _glqeval.shape[0]: _glqeval= None #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vRgl= nsigma*sigmaR1/2.*(glx+1.) vzgl= nsigma*sigmaz1/2.*(glx+1.) vRglw= glw vzglw= glw else: vRgl= nsigma*sigmaR1/2.*(glx12+1.) #vRgl= 1.5/2.*(glx12+1.) vRgl= list(vRgl) vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.)) #vRgl.extend(-1.5/2.*(glx12+1.)) vRgl= numpy.array(vRgl) vzgl= nsigma*sigmaz1/2.*(glx12+1.) #vzgl= 1.5/2.*(glx12+1.) vzgl= list(vzgl) vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.)) #vzgl.extend(-1.5/2.*(glx12+1.)) vzgl= numpy.array(vzgl) vRglw= glw12 vRglw= list(vRglw) vRglw.extend(glw12) vRglw= numpy.array(vRglw) vzglw= glw12 vzglw= list(vzglw) vzglw.extend(glw12) vzglw= numpy.array(vzglw) if 'vTmax' in kwargs: vTmax = kwargs['vTmax'] else: vTmax = 1.5 vTgl= vTmax/2.*(glx+1.) #Tile everything vTgl= numpy.tile(vTgl,(ngl,ngl,1)).T vRgl= numpy.tile(numpy.reshape(vRgl,(1,ngl)).T,(ngl,1,ngl)) vzgl= numpy.tile(vzgl,(ngl,ngl,1)) vTglw= numpy.tile(glw,(ngl,ngl,1)).T #also tile weights vRglw= numpy.tile(numpy.reshape(vRglw,(1,ngl)).T,(ngl,1,ngl)) vzglw= numpy.tile(vzglw,(ngl,ngl,1)) #evaluate if _glqeval is None and _jr is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R+numpy.zeros(ngl*ngl*ngl), vRgl.flatten(), vTgl.flatten(), z+numpy.zeros(ngl*ngl*ngl), vzgl.flatten(), log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl)) elif not _jr is None and _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz), log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl)) elif not _jr is None and not _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz), rg=_rg,kappa=_kappa,nu=_nu, Omega=_Omega, log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl)) else: logqeval= _glqeval if _returngl: return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2, logqeval) elif _return_actions and _return_freqs: return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2, jr,lz,jz, rg,kappa,nu,Omega) elif _return_actions: return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2, jr,lz,jz) else: return numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o *vTglw*vRglw*vzglw*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2) elif mc: mvT= (thisvc-va)/gamma/sigmaR1 if _vrs is None: vrs= numpy.random.normal(size=nmc) else: vrs= _vrs if _vts is None: vts= numpy.random.normal(size=nmc)+mvT else: if _rawgausssamples: vts= _vts+mvT else: vts= _vts if _vzs is None: vzs= numpy.random.normal(size=nmc) else: vzs= _vzs Is= _vmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R, numpy.ones(nmc)*z, self,sigmaR1,gamma,sigmaz1,mvT, n,m,o) if _returnmc: if _rawgausssamples: return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o), vrs,vts-mvT,vzs) else: return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o), vrs,vts,vzs) else: return numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o) else: #pragma: no cover because this is too slow; a warning is shown warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning) return integrate.tplquad(_vmomentsurfaceIntegrand, 1./gamma*(thisvc-va)/sigmaR1-nsigma, 1./gamma*(thisvc-va)/sigmaR1+nsigma, lambda x: 0., lambda x: nsigma, lambda x,y: 0., lambda x,y: nsigma, (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o), **kwargs)[0]*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
def jmomentdensity(self,*args,**kwargs): """ NAME: jmomentdensity PURPOSE: calculate the an arbitrary moment of an action of the velocity distribution at R times the surfacmass INPUT: R - radius at which to calculate the moment(/ro) n - jr^n m - lz^m o - jz^o OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over (when doing explicit numerical integral) mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples OUTPUT: <jr^n lz^m jz^o x density> at R (no support for units) HISTORY: 2012-08-09 - Written - Bovy (IAS@MPIA) """ use_physical= kwargs.pop('use_physical',True) ro= kwargs.pop('ro',None) if ro is None and hasattr(self,'_roSet') and self._roSet: ro= self._ro ro= parse_length_kpc(ro) vo= kwargs.pop('vo',None) if vo is None and hasattr(self,'_voSet') and self._voSet: vo= self._vo vo= parse_velocity_kms(vo) if use_physical and not vo is None and not ro is None: fac= (ro*vo)**(args[2]+args[3]+args[4])/ro**3 if _APY_UNITS: u= 1/units.kpc**3*(units.kpc*units.km/units.s)**(args[2]+args[3]+args[4]) out= self._jmomentdensity(*args,**kwargs) if _APY_UNITS: return units.Quantity(out*fac,unit=u) else: return out*fac else: return self._jmomentdensity(*args,**kwargs)
[docs] def _jmomentdensity(self,R,z,n,m,o,nsigma=None,mc=True,nmc=10000, _returnmc=False,_vrs=None,_vts=None,_vzs=None, **kwargs): """Non-physical version of jmomentdensity, otherwise the same""" if nsigma == None: nsigma= _NSIGMA sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr) sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz) thisvc= potential.vcirc(self._pot,R,use_physical=False) #Use the asymmetric drift equation to estimate va gamma= numpy.sqrt(0.5) va= sigmaR1**2./2./thisvc\ *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5 +R*(1./self._hr+2./self._hsr)) if numpy.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center if mc: mvT= (thisvc-va)/gamma/sigmaR1 if _vrs is None: vrs= numpy.random.normal(size=nmc) else: vrs= _vrs if _vts is None: vts= numpy.random.normal(size=nmc)+mvT else: vts= _vts if _vzs is None: vzs= numpy.random.normal(size=nmc) else: vzs= _vzs Is= _jmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R,numpy.ones(nmc)*z,self,sigmaR1,gamma,sigmaz1,mvT,n,m,o) if _returnmc: return (numpy.mean(Is)*sigmaR1**2.*gamma*sigmaz1, vrs,vts,vzs) else: return numpy.mean(Is)*sigmaR1**2.*gamma*sigmaz1 else: #pragma: no cover because this is too slow; a warning is shown warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning) return integrate.tplquad(_jmomentsurfaceIntegrand, 1./gamma*(thisvc-va)/sigmaR1-nsigma, 1./gamma*(thisvc-va)/sigmaR1+nsigma, lambda x: 0., lambda x: nsigma, lambda x,y: 0., lambda x,y: nsigma, (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o), **kwargs)[0]*sigmaR1**2.*gamma*sigmaz1
[docs] @potential_physical_input @physical_conversion('numberdensity',pop=True) def density(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: density PURPOSE: calculate the density at R,z by marginalizing over velocity INPUT: R - radius at which to calculate the density (can be Quantity) z - height at which to calculate the density (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: density at (R,z) HISTORY: 2012-07-26 - Written - Bovy (IAS@MPIA) """ return self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, gl=gl,ngl=ngl, **kwargs)
[docs] @potential_physical_input @physical_conversion('velocity2',pop=True) def sigmaR2(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: sigmaR2 PURPOSE: calculate sigma_R^2 by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: sigma_R^2 HISTORY: 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,2.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,2.,0.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,2.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('velocity2',pop=True) def sigmaRz(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: sigmaRz PURPOSE: calculate sigma_RZ^2 by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: sigma_Rz^2 HISTORY: 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,1.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,1.,0.,1., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,1.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('angle',pop=True) def tilt(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: tilt PURPOSE: calculate the tilt of the velocity ellipsoid by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: tilt in rad HISTORY: 2012-12-23 - Written - Bovy (IAS) 2017-10-28 - Changed return unit to rad - Bovy (UofT) """ warnings.warn("In versions >1.3, the output unit of quasiisothermaldf.tilt has been changed to radian (from degree before)",galpyWarning) if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) tsigmar2= self._vmomentdensity(R,z,2.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass tsigmaz2= self._vmomentdensity(R,z,0.,0.,2., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass tsigmarz= self._vmomentdensity(R,z,1.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass return 0.5*numpy.arctan(2.*tsigmarz/(tsigmar2-tsigmaz2)) elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) tsigmar2= self._vmomentdensity(R,z,2.,0.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass tsigmaz2= self._vmomentdensity(R,z,0.,0.,2., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass tsigmarz= self._vmomentdensity(R,z,1.,0.,1., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass return 0.5*numpy.arctan(2.*tsigmarz/(tsigmar2-tsigmaz2)) else: raise NotImplementedError("Use either mc=True or gl=True")
[docs] @potential_physical_input @physical_conversion('velocity2',pop=True) def sigmaz2(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: sigmaz2 PURPOSE: calculate sigma_z^2 by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: sigma_z^2 HISTORY: 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,0.,0.,2., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,0.,0.,2., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,0.,0.,2., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('velocity',pop=True) def meanvT(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: meanvT PURPOSE: calculate the mean rotational velocity by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: meanvT HISTORY: 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,0.,1.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,0.,1.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,0.,1.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('velocity',pop=True) def meanvR(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: meanvR PURPOSE: calculate the mean radial velocity by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: meanvR HISTORY: 2012-12-23 - Written - Bovy (IAS) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,1.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,1.,0.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,1.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('velocity',pop=True) def meanvz(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: meanvz PURPOSE: calculate the mean vertical velocity by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: meanvz HISTORY: 2012-12-23 - Written - Bovy (IAS) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._vmomentdensity(R,z,0.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) return self._vmomentdensity(R,z,0.,0.,1., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._vmomentdensity(R,z,0.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('velocity2',pop=True) def sigmaT2(self,R,z,nsigma=None,mc=False,nmc=10000, gl=True,ngl=_DEFAULTNGL,**kwargs): """ NAME: sigmaT2 PURPOSE: calculate sigma_T^2 by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples gl= if True, calculate using Gauss-Legendre integration ngl= if gl, use ngl-th order Gauss-Legendre integration for each dimension OUTPUT: sigma_T^2 HISTORY: 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) mvt= self._vmomentdensity(R,z,0.,1.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass return self._vmomentdensity(R,z,0.,2.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass\ -mvt**2. elif gl: surfmass, glqeval= self._vmomentdensity(R,z,0.,0.,0., gl=gl,ngl=ngl, _returngl=True, **kwargs) mvt= self._vmomentdensity(R,z,0.,1.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass return self._vmomentdensity(R,z,0.,2.,0., ngl=ngl,gl=gl, _glqeval=glqeval, **kwargs)/surfmass-mvt**2. else: #pragma: no cover because this is too slow; a warning is shown surfmass= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs) return (self._vmomentdensity(R,z,0.,2.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/surfmass\ -(self._vmomentdensity(R,z,0.,2.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/surfmass)**2.)
[docs] @potential_physical_input @physical_conversion('action',pop=True) def meanjr(self,R,z,nsigma=None,mc=True,nmc=10000,**kwargs): """ NAME: meanjr PURPOSE: calculate the mean radial action by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples OUTPUT: meanjr HISTORY: 2012-08-09 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._jmomentdensity(R,z,1.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._jmomentdensity(R,z,1.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('action',pop=True) def meanlz(self,R,z,nsigma=None,mc=True,nmc=10000,**kwargs): """ NAME: meanlz PURPOSE: calculate the mean angular momemtum by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples OUTPUT: meanlz HISTORY: 2012-08-09 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._jmomentdensity(R,z,0.,1.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._jmomentdensity(R,z,0.,1.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input @physical_conversion('action',pop=True) def meanjz(self,R,z,nsigma=None,mc=True,nmc=10000,**kwargs): """ NAME: meanjz PURPOSE: calculate the mean vertical action by marginalizing over velocity INPUT: R - radius at which to calculate this (can be Quantity) z - height at which to calculate this (can be Quantity) OPTIONAL INPUT: nsigma - number of sigma to integrate the velocities over scipy.integrate.tplquad kwargs epsabs and epsrel mc= if True, calculate using Monte Carlo integration nmc= if mc, use nmc samples OUTPUT: meanjz HISTORY: 2012-08-09 - Written - Bovy (IAS@MPIA) """ if mc: surfmass, vrs, vts, vzs= self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=True, **kwargs) return self._jmomentdensity(R,z,0.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc,_returnmc=False, _vrs=vrs,_vts=vts,_vzs=vzs, **kwargs)/surfmass else: #pragma: no cover because this is too slow; a warning is shown return (self._jmomentdensity(R,z,0.,0.,1., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs)/ self._vmomentdensity(R,z,0.,0.,0., nsigma=nsigma,mc=mc,nmc=nmc, **kwargs))
[docs] @potential_physical_input def sampleV(self,R,z,n=1,**kwargs): """ NAME: sampleV PURPOSE: sample a radial, azimuthal, and vertical velocity at R,z INPUT: R - Galactocentric distance (can be Quantity) z - height (can be Quantity) n= number of distances to sample OUTPUT: list of samples HISTORY: 2012-12-17 - Written - Bovy (IAS) """ use_physical= kwargs.pop('use_physical',True) vo= kwargs.pop('vo',None) if vo is None and hasattr(self,'_voSet') and self._voSet: vo= self._vo vo= parse_velocity_kms(vo) #Determine the maximum of the velocity distribution maxVR= 0. maxVz= 0. # scipy 1.5.0: issue scipy#12298: fmin_powell now returns multiD array, # so squeeze out single dimensions by hand maxVT= numpy.squeeze(\ optimize.fmin_powell((lambda x: -self(R,0.,x,z,0.,log=True, use_physical=False)), 1.)) logmaxVD= self(R,maxVR,maxVT,z,maxVz,log=True,use_physical=False) #Now rejection-sample vRs= [] vTs= [] vzs= [] while len(vRs) < n: nmore= n-len(vRs)+1 #sample propvR= numpy.random.normal(size=nmore)*2.*self._sr propvT= numpy.random.normal(size=nmore)*2.*self._sr+maxVT propvz= numpy.random.normal(size=nmore)*2.*self._sz VDatprop= self(R+numpy.zeros(nmore), propvR,propvT,z+numpy.zeros(nmore), propvz,log=True,use_physical=False)-logmaxVD VDatprop-= -0.5*(propvR**2./4./self._sr**2.+propvz**2./4./self._sz**2.\ +(propvT-maxVT)**2./4./self._sr**2.) VDatprop= numpy.reshape(VDatprop,(nmore)) indx= (VDatprop > numpy.log(numpy.random.random(size=nmore))) #accept vRs.extend(list(propvR[indx])) vTs.extend(list(propvT[indx])) vzs.extend(list(propvz[indx])) out= numpy.empty((n,3)) out[:,0]= vRs[0:n] out[:,1]= vTs[0:n] out[:,2]= vzs[0:n] if use_physical and not vo is None: if _APY_UNITS: return units.Quantity(out*vo,unit=units.km/units.s) else: return out*vo else: return out
[docs] @potential_physical_input def sampleV_interpolate(self,R,z,R_pixel,z_pixel,num_std=3,R_min=None, R_max=None,z_max=None,**kwargs): """ NAME: sampleV_interpolate PURPOSE: Given an array of R and z coordinates of stars, return the positions and their radial, azimuthal, and vertical velocity. INPUT: R - array of Galactocentric distance (can be Quantity) z - array of height (can be Quantity) R_pixel, z_pixel= the pixel size for creating the grid for interpolation (in natural unit) num_std= number of standard deviation to be considered outliers sampled separately from interpolation R_min, R_max, z_max= optional edges of the grid OUTPUT: coord_v= a numpy array containing the sampled velocity, (vR, vT, vz), where each row correspond to the row of (R,z) HISTORY: 2018-08-10 - Written - Samuel Wong (University of Toronto) """ use_physical= kwargs.pop('use_physical',True) vo= kwargs.pop('vo',None) if vo is None and hasattr(self,'_voSet') and self._voSet: vo= self._vo vo= parse_velocity_kms(vo) #Initialize output array coord_v= numpy.empty((numpy.size(R), 3)) #Since the sign of z doesn't matter, work with absolute value of z z= numpy.abs(z) # Grid edges if R_min is None: R_min= numpy.amax([numpy.mean(R)-num_std*numpy.std(R), numpy.amin(R)]) if R_max is None: R_max= numpy.amin([numpy.mean(R)+num_std*numpy.std(R), numpy.amax(R)]) if z_max is None: z_max= numpy.amin([numpy.mean(z)+num_std*numpy.std(z), numpy.amax(z)]) z_min= 0. #Always start grid at z=0 for stars close to plane #Separate the coodinates into outliers and normal points #Define outliers as points outside of grid mask= numpy.any([R < R_min, R > R_max, z > z_max],axis = 0) outliers_R= R[mask] outliers_z= z[mask] normal_R= R[~mask] normal_z= z[~mask] #Sample the velocity of outliers directly (without interpolation) outlier_coord_v= numpy.empty((outliers_R.size, 3)) for i in range(outliers_R.size): outlier_coord_v[i]= self.sampleV(outliers_R[i], outliers_z[i], use_physical=False)[0] #Prepare for optimizing maxVT on a grid #Get the new hash of the parameters of grid new_hash= hashlib.md5(numpy.array([R_min,R_max,z_max,R_pixel,z_pixel])).hexdigest() #Reuse old interpolated object if new hash matches the old one if new_hash == self._maxVT_hash: ip_max_vT= self._maxVT_ip #Generate a new interpolation object if different from before else: R_number= int((R_max - R_min)/R_pixel) z_number= int((z_max - z_min)/z_pixel) R_linspace= numpy.linspace(R_min, R_max, R_number) z_linspace= numpy.linspace(z_min, z_max, z_number) Rv, zv= numpy.meshgrid(R_linspace, z_linspace) grid= numpy.dstack((Rv, zv)) #This grid stores (R,z) coordinate #Grid is a 3 dimensional array since it stores pairs of values, but #grid max vT is a 2 dimensinal array grid_max_vT= numpy.empty((grid.shape[0], grid.shape[1])) #Optimize max_vT on the grid for i in range(z_number): for j in range(R_number): R, z= grid[i][j] grid_max_vT[i][j]= numpy.squeeze(\ optimize.fmin_powell((lambda x: -self( R,0.,x,z,0.,log=True, use_physical=False)),1.)) #Determine degree of interpolation ky= numpy.min([R_number-1,3]) kx= numpy.min([z_number-1,3]) #Generate interpolation object ip_max_vT= interpolate.RectBivariateSpline(z_linspace,R_linspace, grid_max_vT,kx=kx,ky=ky) #Store interpolation object self._maxVT_ip= ip_max_vT #Update hash of parameters self._maxVT_hash= new_hash #Evaluate interpolation object to get maxVT at the normal coordinates normal_max_vT= ip_max_vT.ev(normal_z, normal_R) #Sample all 3 velocities at a normal point and use interpolated vT normal_coord_v= \ self._sampleV_preoptimized(normal_R,normal_z,normal_max_vT) #Combine normal and outlier result, preserving original order coord_v[mask]= outlier_coord_v coord_v[~mask]= normal_coord_v if use_physical and not vo is None: if _APY_UNITS: return units.Quantity(coord_v*vo,unit=units.km/units.s) else: return coord_v*vo else: return coord_v
def _sampleV_preoptimized(self,R,z,maxVT): """ NAME: _sampleV_preoptimized PURPOSE: sample a radial, azimuthal, and vertical velocity at R,z; R,z can be an array of positions maxVT is already optimized INPUT: R - Galactocentric distance (can be Quantity) z - height (can be Quantity) maxVT - an array of pre-optimized maximum vT at corresponding R,z OUTPUT: a numpy array containing the sampled velocity, (vR, vT, vz), where each row correspond to the row of (R,z) HISTORY: 2018-08-09 - Written - Samuel Wong (University of Toronto) """ length = numpy.size(R) out= numpy.empty((length,3)) #Initialize output #Determine the maximum of the velocity distribution maxVR= numpy.zeros(length) maxVz= numpy.zeros(length) logmaxVD= self(R,maxVR,maxVT,z,maxVz,log=True,use_physical=False) #Now rejection-sample #Intiialize boolean index of position remaining to be sampled remain_indx = numpy.full(length,True) while numpy.any(remain_indx): nmore= numpy.sum(remain_indx) propvR= numpy.random.normal(size=nmore)*2.*self._sr propvT= numpy.random.normal(size=nmore)*2.*self._sr+maxVT[remain_indx] propvz= numpy.random.normal(size=nmore)*2.*self._sz VDatprop= self(R[remain_indx],propvR,propvT,z[remain_indx],propvz, log=True, use_physical=False)-logmaxVD[remain_indx] VDatprop-= -0.5*(propvR**2./4./self._sr**2.+ propvz**2./4./self._sz**2.+ (propvT-maxVT[remain_indx])**2./4./self._sr**2.) accept_indx= (VDatprop > numpy.log(numpy.random.random(size=nmore))) vR_accept= propvR[accept_indx] vT_accept= propvT[accept_indx] vz_accept= propvz[accept_indx] #Get the indexing of rows of output array that need to be updated #with newly accepted velocity to_change= numpy.copy(remain_indx) to_change[remain_indx]= accept_indx out[to_change]= numpy.stack((vR_accept,vT_accept,vz_accept), axis = 1) #Removing accepted sampled from remain index remain_indx[remain_indx]= ~accept_indx return out
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity2',pop=True) def pvR(self,vR,R,z,gl=True,ngl=_DEFAULTNGL2,nsigma=4.,vTmax=1.5): """ NAME: pvR PURPOSE: calculate the marginalized vR probability at this location (NOT normalized by the density) INPUT: vR - radial velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration nsigma - sets integration limits to [-1,+1]*nsigma*sigma_z(R) for integration over vz (default: 4) vTmax - sets integration limits to [0,vTmax] for integration over vT (default: 1.5) OUTPUT: p(vR,R,z) HISTORY: 2012-12-22 - Written - Bovy (IAS) """ sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz) if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vzgl= nsigma*sigmaz1/2.*(glx+1.) vzglw= glw vzfac= nsigma*sigmaz1 #2 x integration over [0,nsigma*sigmaz1] else: vzgl= nsigma*sigmaz1/2.*(glx12+1.) vzgl= list(vzgl) vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.)) vzgl= numpy.array(vzgl) vzglw= glw12 vzglw= list(vzglw) vzglw.extend(glw12) vzglw= numpy.array(vzglw) vzfac = 0.5*nsigma*sigmaz1 #integration over [-nsigma*sigmaz1,0] and [0,nsigma*sigmaz1] vTgl= vTmax/2.*(glx+1.) vTfac= 0.5 * vTmax #integration over [0.,vTmax] #Tile everything vTgl= numpy.tile(vTgl,(ngl,1)).T vzgl= numpy.tile(vzgl,(ngl,1)) vTglw= numpy.tile(glw,(ngl,1)).T #also tile weights vzglw= numpy.tile(vzglw,(ngl,1)) #evaluate logqeval= numpy.reshape(self(R+numpy.zeros(ngl*ngl), vR+numpy.zeros(ngl*ngl), vTgl.flatten(), z+numpy.zeros(ngl*ngl), vzgl.flatten(), log=True, use_physical=False), (ngl,ngl)) return numpy.sum(numpy.exp(logqeval)*vTglw*vzglw*vzfac)*vTfac
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity2',pop=True) def pvT(self,vT,R,z,gl=True,ngl=_DEFAULTNGL2,nsigma=4.): """ NAME: pvT PURPOSE: calculate the marginalized vT probability at this location (NOT normalized by the density) INPUT: vT - tangential velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration nsigma - sets integration limits to [-1,+1]*nsigma*sigma(R) for integration over vz and vR (default: 4) OUTPUT: p(vT,R,z) HISTORY: 2012-12-22 - Written - Bovy (IAS) 2018-01-12 - Added Gauss-Legendre integration prefactor nsigma^2/4 - Trick (MPA) """ sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr) sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz) if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vRgl= nsigma*sigmaR1/2.*(glx+1.) vzgl= nsigma*sigmaz1/2.*(glx+1.) vRglw= glw vzglw= glw vRfac= nsigma*sigmaR1 #2 x integration over [0,nsigma*sigmaR1] vzfac= nsigma*sigmaz1 #2 x integration over [0,nsigma*sigmaz1] else: vRgl= nsigma*sigmaR1/2.*(glx12+1.) vRgl= list(vRgl) vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.)) vRgl= numpy.array(vRgl) vzgl= nsigma*sigmaz1/2.*(glx12+1.) vzgl= list(vzgl) vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.)) vzgl= numpy.array(vzgl) vRglw= glw12 vRglw= list(vRglw) vRglw.extend(glw12) vRglw= numpy.array(vRglw) vzglw= glw12 vzglw= list(vzglw) vzglw.extend(glw12) vzglw= numpy.array(vzglw) vRfac = 0.5*nsigma*sigmaR1 #integration over [-nsigma*sigmaR1,0] and [0,nsigma*sigmaR1] vzfac = 0.5*nsigma*sigmaz1 #integration over [-nsigma*sigmaz1,0] and [0,nsigma*sigmaz1] #Tile everything vRgl= numpy.tile(vRgl,(ngl,1)).T vzgl= numpy.tile(vzgl,(ngl,1)) vRglw= numpy.tile(vRglw,(ngl,1)).T #also tile weights vzglw= numpy.tile(vzglw,(ngl,1)) #evaluate logqeval= numpy.reshape(self(R+numpy.zeros(ngl*ngl), vRgl.flatten(), vT+numpy.zeros(ngl*ngl), z+numpy.zeros(ngl*ngl), vzgl.flatten(), log=True, use_physical=False), (ngl,ngl)) return numpy.sum(numpy.exp(logqeval)*vRglw*vzglw*vRfac*vzfac)
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity2',pop=True) def pvz(self,vz,R,z,gl=True,ngl=_DEFAULTNGL2, nsigma=4.,vTmax=1.5, _return_actions=False,_jr=None,_lz=None,_jz=None, _return_freqs=False, _rg=None,_kappa=None,_nu=None,_Omega=None, _sigmaR1=None): """ NAME: pvz PURPOSE: calculate the marginalized vz probability at this location (NOT normalized by the density) INPUT: vz - vertical velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration nsigma - sets integration limits to [-1,+1]*nsigma*sigma_R(R) for integration over vR (default: 4) vTmax - sets integration limits to [0,vTmax] for integration over vT (default: 1.5) OUTPUT: p(vz,R,z) HISTORY: 2012-12-22 - Written - Bovy (IAS) """ if _sigmaR1 is None: sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr) else: sigmaR1= _sigmaR1 if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vRgl= (glx+1.) vRglw= glw vRfac= nsigma*sigmaR1 #2 x integration over [0,nsigma*sigmaR1] else: vRgl= (glx12+1.) vRgl= list(vRgl) vRgl.extend(-(glx12+1.)) vRgl= numpy.array(vRgl) vRglw= glw12 vRglw= list(vRglw) vRglw.extend(glw12) vRglw= numpy.array(vRglw) vRfac = 0.5*nsigma*sigmaR1 #integration over [-nsigma*sigmaR1,0] and [0,nsigma*sigmaR1] vTgl= vTmax/2.*(glx+1.) vTfac= 0.5 * vTmax #integration over [0.,vTmax] #Tile everything vTgl= numpy.tile(vTgl,(ngl,1)).T vRgl= numpy.tile(vRgl,(ngl,1)) vTglw= numpy.tile(glw,(ngl,1)).T #also tile weights vRglw= numpy.tile(vRglw,(ngl,1)) #If inputs are arrays, tile if isinstance(R,numpy.ndarray): nR= len(R) R= numpy.tile(R,(ngl,ngl,1)).T.flatten() z= numpy.tile(z,(ngl,ngl,1)).T.flatten() vz= numpy.tile(vz,(ngl,ngl,1)).T.flatten() vTgl= numpy.tile(vTgl,(nR,1,1)).flatten() vRgl= numpy.tile(vRgl,(nR,1,1)).flatten() vTglw= numpy.tile(vTglw,(nR,1,1)) vRglw= numpy.tile(vRglw,(nR,1,1)) scalarOut= False else: R= R+numpy.zeros(ngl*ngl) z= z+numpy.zeros(ngl*ngl) vz= vz+numpy.zeros(ngl*ngl) nR= 1 scalarOut= True vRgl= vRgl.flatten() vRgl*= numpy.tile(nsigma*sigmaR1/2.,(ngl,ngl,1)).T.flatten() #evaluate if _jr is None and _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R, vRgl.flatten(), vTgl.flatten(), z, vz, log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(nR,ngl*ngl)) elif not _jr is None and not _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz), rg=_rg,kappa=_kappa,nu=_nu, Omega=_Omega, log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(nR,ngl*ngl)) elif not _jr is None and _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz), log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(nR,ngl*ngl)) elif _jr is None and not _rg is None: logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R, vRgl.flatten(), vTgl.flatten(), z, vz, rg=_rg,kappa=_kappa,nu=_nu, Omega=_Omega, log=True, _return_actions=True, _return_freqs=True, use_physical=False) logqeval= numpy.reshape(logqeval,(nR,ngl*ngl)) vRglw= numpy.reshape(vRglw,(nR,ngl*ngl)) vTglw= numpy.reshape(vTglw,(nR,ngl*ngl)) if scalarOut: result= numpy.sum(numpy.exp(logqeval)*vTglw*vRglw,axis=1)[0]*vRfac*vTfac else: result= numpy.sum(numpy.exp(logqeval)*vTglw*vRglw,axis=1)*vRfac*vTfac if _return_actions and _return_freqs: return (result, jr,lz,jz, rg, kappa, nu, Omega) elif _return_freqs: return (result, rg, kappa, nu, Omega) elif _return_actions: return (result, jr,lz,jz) else: return result
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity',pop=True) def pvRvT(self,vR,vT,R,z,gl=True,ngl=_DEFAULTNGL2,nsigma=4.): """ NAME: pvRvT PURPOSE: calculate the marginalized (vR,vT) probability at this location (NOT normalized by the density) INPUT: vR - radial velocity (can be Quantity) vT - tangential velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration nsigma - sets integration limits to [-1,+1]*nsigma*sigma_z(R) for integration over vz (default: 4) OUTPUT: p(vR,vT,R,z) HISTORY: 2013-01-02 - Written - Bovy (IAS) 2018-01-12 - Added Gauss-Legendre integration prefactor nsigma/2 - Trick (MPA) """ sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz) if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vzgl= nsigma*sigmaz1/2.*(glx+1.) vzglw= glw vzfac= nsigma*sigmaz1 #2 x integration over [0,nsigma*sigmaz1] else: vzgl= nsigma*sigmaz1/2.*(glx12+1.) vzgl= list(vzgl) vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.)) vzgl= numpy.array(vzgl) vzglw= glw12 vzglw= list(vzglw) vzglw.extend(glw12) vzglw= numpy.array(vzglw) vzfac = 0.5*nsigma*sigmaz1 #integration over [-nsigma*sigmaz1,0] and [0,nsigma*sigmaz1] #evaluate logqeval= self(R+numpy.zeros(ngl), vR+numpy.zeros(ngl), vT+numpy.zeros(ngl), z+numpy.zeros(ngl), vzgl, log=True,use_physical=False) return numpy.sum(numpy.exp(logqeval)*vzglw*vzfac)
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity',pop=True) def pvTvz(self,vT,vz,R,z,gl=True,ngl=_DEFAULTNGL2,nsigma=4.): """ NAME: pvTvz PURPOSE: calculate the marginalized (vT,vz) probability at this location (NOT normalized by the density) INPUT: vT - tangential velocity (can be Quantity) vz - vertical velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration nsigma - sets integration limits to [-1,+1]*nsigma*sigma_R(R) for integration over vR (default: 4) OUTPUT: p(vT,vz,R,z) HISTORY: 2012-12-22 - Written - Bovy (IAS) 2018-01-12 - Added Gauss-Legendre integration prefactor nsigma/2 - Trick (MPA) """ sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr) if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere if isinstance(self._aA,(actionAngle.actionAngleAdiabatic, actionAngle.actionAngleAdiabaticGrid)): vRgl= nsigma*sigmaR1/2.*(glx+1.) vRglw= glw vRfac= nsigma*sigmaR1 #2 x integration over [0,nsigma*sigmaR1] else: vRgl= nsigma*sigmaR1/2.*(glx12+1.) vRgl= list(vRgl) vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.)) vRgl= numpy.array(vRgl) vRglw= glw12 vRglw= list(vRglw) vRglw.extend(glw12) vRglw= numpy.array(vRglw) vRfac = 0.5*nsigma*sigmaR1 #integration over [-nsigma*sigmaR1,0] and [0,nsigma*sigmaR1] #evaluate logqeval= self(R+numpy.zeros(ngl), vRgl, vT+numpy.zeros(ngl), z+numpy.zeros(ngl), vz+numpy.zeros(ngl), log=True,use_physical=False) return numpy.sum(numpy.exp(logqeval)*vRglw*vRfac)
[docs] @actionAngle_physical_input @physical_conversion('phasespacedensityvelocity',pop=True) def pvRvz(self,vR,vz,R,z,gl=True,ngl=_DEFAULTNGL2,vTmax=1.5): """ NAME: pvR PURPOSE: calculate the marginalized (vR,vz) probability at this location (NOT normalized by the density) INPUT: vR - radial velocity (can be Quantity) vz - vertical velocity (can be Quantity) R - radius (can be Quantity) z - height (can be Quantity) gl - use Gauss-Legendre integration (True, currently the only option) ngl - order of Gauss-Legendre integration vTmax - sets integration limits to [0,vTmax] for integration over vT (default: 1.5) OUTPUT: p(vR,vz,R,z) HISTORY: 2013-01-02 - Written - Bovy (IAS) 2018-01-12 - Added Gauss-Legendre integration prefactor vTmax/2 - Trick (MPA) """ if gl: if ngl % 2 == 1: raise ValueError("ngl must be even") #Use Gauss-Legendre integration for all if ngl == _DEFAULTNGL: glx, glw= self._glxdef, self._glwdef glx12, glw12= self._glxdef12, self._glwdef12 elif ngl == _DEFAULTNGL2: glx, glw= self._glxdef2, self._glwdef2 glx12, glw12= self._glxdef, self._glwdef else: glx, glw= numpy.polynomial.legendre.leggauss(ngl) glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2) #Evaluate everywhere vTgl= vTmax/2.*(glx+1.) vTglw= glw vTfac= 0.5 * vTmax #integration over [0.,vTmax] #If inputs are arrays, tile if isinstance(R,numpy.ndarray): nR= len(R) R= numpy.tile(R,(ngl,1)).T.flatten() z= numpy.tile(z,(ngl,1)).T.flatten() vR= numpy.tile(vR,(ngl,1)).T.flatten() vz= numpy.tile(vz,(ngl,1)).T.flatten() vTgl= numpy.tile(vTgl,(nR,1)).flatten() vTglw= numpy.tile(vTglw,(nR,1)) scalarOut= False else: R= R+numpy.zeros(ngl) vR= vR+numpy.zeros(ngl) z= z+numpy.zeros(ngl) vz= vz+numpy.zeros(ngl) nR= 1 scalarOut= True #evaluate logqeval= numpy.reshape(self(R, vR, vTgl, z, vz, log=True, use_physical=False), (nR,ngl)) out= numpy.sum(numpy.exp(logqeval)*vTglw*vTfac,axis=1) if scalarOut: return out[0] else: return out
def _calc_epifreq(self,r): """ NAME: _calc_epifreq PURPOSE: calculate the epicycle frequency at r INPUT: r - radius OUTPUT: kappa HISTORY: 2012-07-25 - Written - Bovy (IAS@MPIA) NOTE: takes about 0.1 ms for a Miyamoto-Nagai potential """ return potential.epifreq(self._pot,r) def _calc_verticalfreq(self,r): """ NAME: _calc_verticalfreq PURPOSE: calculate the vertical frequency at r INPUT: r - radius OUTPUT: nu HISTORY: 2012-07-25 - Written - Bovy (IAS@MPIA) NOTE: takes about 0.05 ms for a Miyamoto-Nagai potential """ return potential.verticalfreq(self._pot,r) def _rg(self,lz): """ NAME: _rg PURPOSE: calculate the radius of a circular orbit of Lz INPUT: lz - Angular momentum OUTPUT: radius HISTORY: 2012-07-25 - Written - Bovy (IAS@MPIA) NOTE: seems to take about ~0.5 ms for a Miyamoto-Nagai potential; ~0.75 ms for a MWPotential about the same with or without interpolation of the rotation curve Not sure what to do about negative lz... """ if isinstance(lz,numpy.ndarray): indx= (lz > self._precomputergLzmax)*(lz < self._precomputergLzmin) indxc= True^indx out= numpy.empty(lz.shape) out[indxc]= self._rgInterp(lz[indxc]) out[indx]= numpy.array([potential.rl(self._pot,lz[indx][ii]) for ii in range(numpy.sum(indx))]) return out else: if lz > self._precomputergLzmax or lz < self._precomputergLzmin: return potential.rl(self._pot,lz) return numpy.atleast_1d(self._rgInterp(lz))
def _vmomentsurfaceIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,n,m,o): #pragma: no cover because this is too slow; a warning is shown """Internal function that is the integrand for the vmomentsurface mass integration""" return vR**n*vT**m*vz**o*df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1, use_physical=False) def _vmomentsurfaceMCIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,mvT,n,m,o): """Internal function that is the integrand for the vmomentsurface mass integration""" return vR**n*vT**m*vz**o*df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False)*numpy.exp(vR**2./2.+(vT-mvT)**2./2.+vz**2./2.) def _jmomentsurfaceIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,n,m,o): #pragma: no cover because this is too slow; a warning is shown """Internal function that is the integrand for the vmomentsurface mass integration""" return df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False, func= (lambda x,y,z: x**n*y**m*z**o)) def _jmomentsurfaceMCIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,mvT,n,m,o): """Internal function that is the integrand for the vmomentsurface mass integration""" return df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False, func=(lambda x,y,z: x**n*y**m*z**o))\ *numpy.exp(vR**2./2.+(vT-mvT)**2./2.+vz**2./2.)