From 8be16e011b21db9bdc6cc2ab23b5a454e4bd994f Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Sun, 1 Oct 2023 23:04:51 -0400 Subject: [PATCH] Numpy-style docstrings for DF functions and methods --- galpy/df/constantbetaHernquistdf.py | 64 +- galpy/df/constantbetadf.py | 152 +- galpy/df/df.py | 81 +- galpy/df/diskdf.py | 2188 +++++++++++++------------ galpy/df/eddingtondf.py | 75 +- galpy/df/evolveddiskdf.py | 1252 +++++++------- galpy/df/isotropicHernquistdf.py | 65 +- galpy/df/isotropicNFWdf.py | 71 +- galpy/df/isotropicPlummerdf.py | 63 +- galpy/df/jeans.py | 102 +- galpy/df/kingdf.py | 48 +- galpy/df/osipkovmerrittHernquistdf.py | 68 +- galpy/df/osipkovmerrittNFWdf.py | 71 +- galpy/df/osipkovmerrittdf.py | 165 +- galpy/df/quasiisothermaldf.py | 1856 ++++++++++----------- galpy/df/sphericaldf.py | 505 +++--- galpy/df/streamdf.py | 1502 ++++++++--------- galpy/df/streamgapdf.py | 914 +++++------ galpy/df/streamspraydf.py | 127 +- galpy/df/surfaceSigmaProfile.py | 233 +-- 20 files changed, 4619 insertions(+), 4983 deletions(-) diff --git a/galpy/df/constantbetaHernquistdf.py b/galpy/df/constantbetaHernquistdf.py index fd7b159e9..e0cf88bca 100644 --- a/galpy/df/constantbetaHernquistdf.py +++ b/galpy/df/constantbetaHernquistdf.py @@ -14,27 +14,22 @@ class constantbetaHernquistdf(_constantbetadf): 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) + Initialize a Hernquist DF with constant anisotropy. + + Parameters + ---------- + pot : HernquistPotential + Hernquist potential which determines the DF. + beta : float + Anisotropy parameter. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-07-22 - Written - Lane (UofT) """ assert isinstance( pot, HernquistPotential @@ -56,25 +51,22 @@ def __init__(self, pot=None, beta=0, ro=None, vo=None): def fE(self, E): """ - NAME: - - fE - - PURPOSE - - Calculate the energy portion of a Hernquist distribution function - - INPUT: - - E - The energy (can be Quantity) + Calculate the energy portion of a Hernquist distribution function - OUTPUT: + Parameters + ---------- + E : float, numpy.ndarray, or Quantity + The energy. - fE - The value of the energy portion of the DF + Returns + ------- + float or numpy.ndarray + The value of the energy portion of the DF - HISTORY: + Notes + ----- + - 2020-07-22 - Written - 2020-07-22 - Written """ Etilde = -numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo) / self._psi0) # Handle potential E outside of bounds diff --git a/galpy/df/constantbetadf.py b/galpy/df/constantbetadf.py index b97ceddf9..f955a8aa8 100644 --- a/galpy/df/constantbetadf.py +++ b/galpy/df/constantbetadf.py @@ -23,26 +23,24 @@ def __init__( self, pot=None, denspot=None, beta=None, rmax=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 - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale - Characteristic scale radius to aid sampling calculations. - Not necessary, and will also be overridden by value from pot if - available. - + Initialize a spherical DF with constant anisotropy parameter. + + Parameters + ---------- + pot : Potential or list of Potential instances, optional + Spherical potential which determines the DF. + denspot : Potential or list of Potential instances, optional + Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot). + beta : float, optional + Anisotropy parameter. Default is None. + rmax : float or Quantity, optional + Maximum radius to consider; DF is cut off at E = Phi(rmax). Default is None. + scale : float or Quantity, optional + Characteristic scale radius to aid sampling calculations. Not necessary, and will also be overridden by value from pot if available. Default is None. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). """ anisotropicsphericaldf.__init__( self, pot=pot, denspot=denspot, rmax=rmax, scale=scale, ro=ro, vo=vo @@ -51,27 +49,24 @@ def __init__( def _call_internal(self, *args): """ - NAME: + Evaluate the DF for a constant anisotropy Hernquist. - _call_internal + Parameters + ---------- + E : float + The energy. + L : float + The angular momentum. - PURPOSE: + Returns + ------- + float + The value of the DF. - Evaluate the DF for a constant anisotropy Hernquist + Notes + ----- + - 2020-07-22 - Written - Lane (UofT) - 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) @@ -166,37 +161,30 @@ def __init__( vo=None, ): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a spherical DF with constant anisotropy parameter - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - beta= (0.) anisotropy parameter - - twobeta= (None) twice the anisotropy parameter (useful for \beta = half-integer, which is a special case); has priority over beta - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale - Characteristic scale radius to aid sampling calculations. Optionaland will also be overridden by value from pot if available. - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2021-02-14 - Written - Bovy (UofT) + Initialize a spherical DF with constant anisotropy parameter + + Parameters + ---------- + pot : Potential instance or list thereof, optional + Potential instance or list thereof + denspot : Potential instance or list thereof, optional + Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) + beta : float, optional + anisotropy parameter + twobeta : float, optional + twice the anisotropy parameter (useful for \beta = half-integer, which is a special case); has priority over beta + rmax : float or Quantity, optional + maximum radius to consider; DF is cut off at E = Phi(rmax) + scale : float or Quantity, optional + Characteristic scale radius to aid sampling calculations. Optional and will also be overridden by value from pot if available. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2021-02-14 - Written - Bovy (UofT) """ if not _JAX_LOADED: # pragma: no cover @@ -317,25 +305,21 @@ def sample(self, R=None, z=None, phi=None, n=1, return_orbit=True, rmin=0.0): def fE(self, E): """ - NAME: - - fE - - PURPOSE - - Calculate the energy portion of a constant-beta distribution function - - INPUT: - - E - The energy (can be Quantity) - - OUTPUT: + Calculate the energy portion of a constant-beta distribution function - fE - The value of the energy portion of the DF + Parameters + ---------- + E : float, numpy.ndarray, or Quantity + The energy. - HISTORY: + Returns + ------- + numpy.ndarray + The value of the energy portion of the DF - 2021-02-14 - Written - Bovy (UofT) + Notes + ----- + - 2021-02-14 - Written - Bovy (UofT) """ Eint = numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo)) out = numpy.zeros_like(Eint) diff --git a/galpy/df/df.py b/galpy/df/df.py index 1f5ce370a..7196db2f7 100644 --- a/galpy/df/df.py +++ b/galpy/df/df.py @@ -7,16 +7,19 @@ class df: def __init__(self, ro=None, vo=None): """ - NAME: - __init__ - PURPOSE: - initialize a DF object - INPUT: - ro= (None) distance scale - vo= (None) velocity scale - OUTPUT: - HISTORY: - 2016-02-28 - Written - Bovy (UofT) + Initialize a DF object. + + Parameters + ---------- + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2016-02-28 - Written - Bovy (UofT) + """ # Parse ro and vo if ro is None: @@ -41,25 +44,11 @@ def _check_consistent_units(self): def turn_physical_off(self): """ - NAME: - - turn_physical_off - - PURPOSE: - - turn off automatic returning of outputs in physical units - - INPUT: - - (none) - - OUTPUT: - - (none) - - HISTORY: + Turn off automatic returning of outputs in physical units. - 2017-06-05 - Written - Bovy (UofT) + Notes + ----- + - 2017-06-05 - Written - Bovy (UofT) """ self._roSet = False @@ -68,29 +57,19 @@ def turn_physical_off(self): def turn_physical_on(self, ro=None, vo=None): """ - NAME: - - turn_physical_on - - PURPOSE: - - turn on automatic returning of outputs in physical units - - INPUT: - - ro= reference distance (kpc; can be Quantity) - - vo= reference velocity (km/s; can be Quantity) - - OUTPUT: - - (none) - - HISTORY: - - 2016-06-05 - Written - Bovy (UofT) - - 2020-04-22 - Don't turn on a parameter when it is False - Bovy (UofT) + Turn on automatic returning of outputs in physical units. + + Parameters + ---------- + ro : float or Quantity, optional + Reference distance (kpc). If False, don't turn it on. + vo : float or Quantity, optional + Reference velocity (km/s). If False, don't turn it on. + + Notes + ----- + - 2016-06-05 - Written - Bovy (UofT) + - 2020-04-22 - Don't turn on a parameter when it is False - Bovy (UofT) """ if not ro is False: diff --git a/galpy/df/diskdf.py b/galpy/df/diskdf.py index 6d7ee7ccb..e0438a227 100644 --- a/galpy/df/diskdf.py +++ b/galpy/df/diskdf.py @@ -70,28 +70,34 @@ def __init__( **kwargs ): """ - NAME: - __init__ - PURPOSE: - Initialize a DF - INPUT: - dftype= 'dehnen' or 'corrected-dehnen', 'shu' or 'corrected-shu' - surfaceSigma - instance or class name of the target - surface density and sigma_R profile - (default: both exponential) - profileParams - parameters of the surface and sigma_R profile: - (xD,xS,Sro) where - xD - disk surface mass scalelength / Ro - xS - disk velocity dispersion scalelength / Ro - Sro - disk velocity dispersion at Ro (/vo) - Directly given to the 'surfaceSigmaProfile class, so - could be anything that class takes - beta - power-law index of the rotation curve - correct - correct the DF (i.e., DFcorrection kwargs are also given) - + DFcorrection kwargs (except for those already specified) - OUTPUT: - HISTORY: - 2010-03-10 - Written - Bovy (NYU) + Initialize a DF + + Parameters + ---------- + dftype : str, optional + 'dehnen' or 'corrected-dehnen', 'shu' or 'corrected-shu' + surfaceSigma : instance or class name of the target surface density and sigma_R profile, optional + (default: both exponential) + profileParams : tuple, optional + parameters of the surface and sigma_R profile: (xD,xS,Sro) where + * xD - disk surface mass scalelength / Ro + * xS - disk velocity dispersion scalelength / Ro + * Sro - disk velocity dispersion at Ro (/vo) + Directly given to the 'surfaceSigmaProfile class, so could be anything that class takes + beta : float, optional + power-law index of the rotation curve + correct : bool, optional + correct the DF (i.e., DFcorrection kwargs are also given) + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + **kwargs : dict, optional + DFcorrection kwargs (except for those already specified) + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) """ df.__init__(self, ro=ro, vo=vo) self._dftype = dftype @@ -137,47 +143,34 @@ def __init__( @physical_conversion("phasespacedensity2d", pop=True) def __call__(self, *args, **kwargs): """ - NAME: - - __call__ - - PURPOSE: - - evaluate the distribution function - - INPUT: - - either an orbit instance, a list of such instances, or E,Lz - - 1) Orbit instance or list: - a) Orbit instance alone: use initial condition - b) Orbit instance + t: call the Orbit instance (for list, each instance is called at t) - - 2) - E - energy (/vo^2; or can be Quantity) - L - angular momentun (/ro/vo; or can be Quantity) - - 3) array vxvv [3/4,nt] [must be in natural units /vo,/ro; use Orbit interface for physical-unit input) - - KWARGS: - - marginalizeVperp - marginalize over perpendicular velocity (only supported with 1a) for single orbits above) - - - marginalizeVlos - marginalize over line-of-sight velocity (only supported with 1a) for single orbits above) - - nsigma= number of sigma to integrate over when marginalizing - - +scipy.integrate.quad keywords - - OUTPUT: - - DF(orbit/E,L) - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - + Evaluate the distribution function + + Parameters + ---------- + *args : tuple + Either: + 1) Orbit instance or list: + a) Orbit instance alone: use initial condition + b) Orbit instance + t: call the Orbit instance (for list, each instance is called at t) + 2) E,L - energy (/vo^2; or can be Quantity) and angular momentun (/ro/vo; or can be Quantity) + 3) array vxvv [3/4,nt] [must be in natural units /vo,/ro; use Orbit interface for physical-unit input] + marginalizeVperp : bool, optional + marginalize over perpendicular velocity (only supported with 1a) for single orbits above) + marginalizeVlos : bool, optional + marginalize over line-of-sight velocity (only supported with 1a) for single orbits above) + nsigma : float, optional + number of sigma to integrate over when marginalizing + **kwargs: dict, optional + scipy.integrate.quad keywords + + Returns + ------- + float or numpy.ndarray + value of DF + + Notes + ----- + - 2010-07-10 - Written - Bovy (NYU) """ if isinstance(args[0], Orbit): if len(args[0]) > 1: @@ -404,28 +397,23 @@ def _call_marginalizevlos(self, o, **kwargs): @physical_conversion("velocity2", pop=True) def targetSigma2(self, R, log=False): """ - NAME: - - targetSigma2 - - PURPOSE: - - evaluate the target Sigma_R^2(R) - - INPUT: - - R - radius at which to evaluate (can be Quantity) - - OUTPUT: - - target Sigma_R^2(R) + Evaluate the target Sigma_R^2(R) - log - if True, return the log (default: False) + Parameters + ---------- + R : float or Quantity + Radius at which to evaluate. + log : bool, optional + If True, return the log (default: False). - HISTORY: - - 2010-03-28 - Written - Bovy (NYU) + Returns + ------- + float + Target Sigma_R^2(R). + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) """ return self._surfaceSigmaProfile.sigma2(R, log=log) @@ -433,27 +421,23 @@ def targetSigma2(self, R, log=False): @physical_conversion("surfacedensity", pop=True) def targetSurfacemass(self, R, log=False): """ - NAME: - - targetSurfacemass - - PURPOSE: - - evaluate the target surface mass at R - - INPUT: - - R - radius at which to evaluate (can be Quantity) - - log - if True, return the log (default: False) + Evaluate the target surface mass at R. - OUTPUT: + Parameters + ---------- + R : float or Quantity + Radius at which to evaluate. + log : bool, optional + If True, return the log (default: False). - Sigma(R) + Returns + ------- + float or Quantity + Target surface mass at R. - HISTORY: - - 2010-03-28 - Written - Bovy (NYU) + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) """ return self._surfaceSigmaProfile.surfacemass(R, log=log) @@ -461,32 +445,27 @@ def targetSurfacemass(self, R, log=False): @physical_conversion("surfacedensitydistance", pop=True) def targetSurfacemassLOS(self, d, l, log=False, deg=True): """ - NAME: - - targetSurfacemassLOS - - PURPOSE: - - evaluate the target surface mass along the LOS given l and d - - INPUT: - - d - distance along the line of sight (can be Quantity) + Evaluate the target surface mass along the line of sight given Galactic longitude and distance. - l - Galactic longitude (in deg, unless deg=False; can be Quantity) + Parameters + ---------- + d : float or Quantity + Distance along the line of sight. + l : float or Quantity + Galactic longitude in degrees, unless deg=False. + deg : bool, optional + If False, l is in radians. Default is True. + log : bool, optional + If True, return the logarithm of the surface mass. Default is False. - deg= if False, l is in radians - - log - if True, return the log (default: False) - - OUTPUT: - - Sigma(d,l) x d - - HISTORY: - - 2011-03-23 - Written - Bovy (NYU) + Returns + ------- + float or Quantity + Surface mass times distance. + Notes + ----- + - 2011-03-23 - Written - Bovy (NYU) """ # Calculate R and phi if _APY_LOADED and isinstance(l, units.Quantity): @@ -507,40 +486,33 @@ def surfacemassLOS( self, d, l, deg=True, target=True, romberg=False, nsigma=None, relative=None ): """ - NAME: - - surfacemassLOS - - PURPOSE: - - evaluate the surface mass along the LOS given l and d - - INPUT: - - d - distance along the line of sight (can be Quantity) - - l - Galactic longitude (in deg, unless deg=False; can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - target= if True, use target surfacemass (default) - - romberg - if True, use a romberg integrator (default: False) - - deg= if False, l is in radians - - OUTPUT: - - Sigma(d,l) x d - - HISTORY: - - 2011-03-24 - Written - Bovy (NYU) + Evaluate the surface mass along the line of sight (LOS) given Galactic longitude and distance. + + Parameters + ---------- + d : float or Quantity + Distance along the line of sight. + l : float or Quantity + Galactic longitude (in deg, unless deg=False). + nsigma : float, optional + Number of sigma to integrate the velocities over. + target : bool, optional + If True, use target surfacemass (default). + romberg : bool, optional + If True, use a romberg integrator (default: False). + deg : bool, optional + If False, l is in radians. + relative : bool, optional + If True, return d. + + Returns + ------- + float + Sigma(d,l) x d + Notes + ----- + - 2011-03-24 - Written - Bovy (NYU) """ # Calculate R and phi if _APY_LOADED and isinstance(l, units.Quantity): @@ -571,31 +543,27 @@ def surfacemassLOS( @physical_conversion("position", pop=True) def sampledSurfacemassLOS(self, l, n=1, maxd=None, target=True): """ - NAME: - - sampledSurfacemassLOS - - PURPOSE: + Sample a distance along the line of sight - sample a distance along the line of sight + Parameters + ---------- + l : float or Quantity + Galactic longitude. + n : int, optional + Number of distances to sample. + maxd : float or Quantity, optional + Maximum distance to consider (for the rejection sampling). + target : bool, optional + If True, sample from the 'target' surface mass density, rather than the actual surface mass density (default=True). - INPUT: + Returns + ------- + list + List of samples. - l - Galactic longitude (in rad; can be Quantity) - - n= number of distances to sample - - maxd= maximum distance to consider (for the rejection sampling) (can be Quantity) - - target= if True, sample from the 'target' surface mass density, rather than the actual surface mass density (default=True) - - OUTPUT: - - list of samples - - HISTORY: - - 2011-03-24 - Written - Bovy (NYU) + Notes + ----- + - 2011-03-24 - Written - Bovy (NYU) """ # First calculate where the maximum is @@ -640,36 +608,27 @@ def sampledSurfacemassLOS(self, l, n=1, maxd=None, target=True): @physical_conversion("velocity", pop=True) def sampleVRVT(self, R, n=1, nsigma=None, target=True): """ - NAME: - - sampleVRVT - - PURPOSE: - - sample a radial and azimuthal velocity at R - - INPUT: - - R - Galactocentric distance (can be Quantity) - - n= number of distances to sample + Sample a radial and azimuthal velocity at R - nsigma= number of sigma to rejection-sample on + Parameters + ---------- + R : float or Quantity + Galactocentric distance. + n : int, optional + Number of distances to sample. + nsigma : float, optional + Number of sigma to rejection-sample on. + target : bool, optional + If True, sample using the 'target' sigma_R rather than the actual sigma_R (default=True). - target= if True, sample using the 'target' sigma_R rather than the actual sigma_R (default=True) - - OUTPUT: - - list of samples - - BUGS: - - should use the fact that vR and vT separate - - HISTORY: - - 2011-03-24 - Written - Bovy (NYU) + Returns + ------- + list + List of samples. + Notes + ----- + - 2011-03-24 - Written - Bovy (NYU) """ # Determine where the max of the v-distribution is using asymmetric drift maxVR = 0.0 @@ -706,36 +665,34 @@ def sampleLOS( targetSigma2=True, ): """ - NAME: - - sampleLOS - - PURPOSE: - - sample along a given LOS - - INPUT: - - los - line of sight (in deg, unless deg=False; can be Quantity) - - n= number of desired samples - - deg= los in degrees? (default=True) - - targetSurfmass, targetSigma2= if True, use target surface mass and sigma2 profiles, respectively (there is not much point to doing the latter) - (default=True) - - OUTPUT: - - returns list of Orbits - - BUGS: - target=False uses target distribution for derivatives (this is a detail) - - HISTORY: - - 2011-03-24 - Started - Bovy (NYU) - + Sample along a given LOS + + Parameters + ---------- + los : float or Quantity + Line of sight Galactic longitude. + n : int, optional + Number of distances to sample. + deg : bool, optional + If False, los is in radians. + maxd : float or Quantity, optional + Maximum distance to consider (for the rejection sampling). + nsigma : float, optional + Number of sigma to integrate the velocities over. + targetSurfmass : bool, optional + If True, use target surface mass (default=True). + targetSigma2 : bool, optional + If True, use target sigma_R^2 (default=True). + + Returns + ------- + list + List of Orbits sampled. + + Notes + ----- + - target=False uses target distribution for derivatives (this is a detail) + - 2011-03-24 - Written - Bovy (NYU) """ if _APY_LOADED and isinstance(los, units.Quantity): l = conversion.parse_angle(los) @@ -767,26 +724,21 @@ def sampleLOS( @physical_conversion("velocity", pop=True) def asymmetricdrift(self, R): """ - NAME: - - asymmetricdrift - - PURPOSE: - - estimate the asymmetric drift (vc-mean-vphi) from an approximation to the Jeans equation - - INPUT: - - R - radius at which to calculate the asymmetric drift (can be Quantity) + Estimate the asymmetric drift (vc-mean-vphi) from an approximation to the Jeans equation. - OUTPUT: + Parameters + ---------- + R : float or Quantity + Radius at which to calculate the asymmetric drift. - asymmetric drift at R - - HISTORY: - - 2011-04-02 - Written - Bovy (NYU) + Returns + ------- + float + Asymmetric drift at R. + Notes + ----- + - 2011-04-02 - Written - Bovy (NYU). """ sigmaR2 = self.targetSigma2(R, use_physical=False) return ( @@ -805,34 +757,27 @@ def asymmetricdrift(self, R): @physical_conversion("surfacedensity", pop=True) def surfacemass(self, R, romberg=False, nsigma=None, relative=False): """ - NAME: - - surfacemass - - PURPOSE: - - calculate the surface-mass at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate the surfacemass density (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - surface mass at R + Calculate the surface-mass at R by marginalizing over velocity - HISTORY: + Parameters + ---------- + R : float or Quantity + Radius at which to calculate the surfacemass density. + romberg : bool, optional + If True, use a romberg integrator (default: False) + nsigma : float, optional + Number of sigma to integrate the velocities over + relative : bool, optional + If True, return the relative surface mass at R (default: False) - 2010-03-XX - Written - Bovy (NYU) + Returns + ------- + float + Surface mass at R + Notes + ----- + - 2011-03-XX - Bovy (NYU) """ if nsigma == None: nsigma = _NSIGMA @@ -891,34 +836,27 @@ def surfacemass(self, R, romberg=False, nsigma=None, relative=False): @physical_conversion("velocity2surfacedensity", pop=True) def sigma2surfacemass(self, R, romberg=False, nsigma=None, relative=False): """ + Calculate the product sigma_R^2 x surface-mass at R by marginalizing over velocity. - NAME: + Parameters + ---------- + R : float or Quantity + Radius at which to calculate the sigma_R^2 x surfacemass density. + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : float, optional + Number of sigma to integrate the velocities over. + relative : bool, optional + If True, return the relative density (default: False). - sigma2surfacemass + Returns + ------- + float + Sigma_R^2 x surface-mass at R. - PURPOSE: - - calculate the product sigma_R^2 x surface-mass at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate the sigma_R^2 x surfacemass density (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - sigma_R^2 x surface-mass at R - - HISTORY: - - 2010-03-XX - Written - Bovy (NYU) + Notes + ----- + - 2010-03-XX - Written - Bovy (NYU). """ if nsigma == None: @@ -976,41 +914,31 @@ def sigma2surfacemass(self, R, romberg=False, nsigma=None, relative=False): def vmomentsurfacemass(self, *args, **kwargs): """ - NAME: - - vmomentsurfacemass - - PURPOSE: - - calculate the an arbitrary moment of the velocity distribution - at R times the surfacmass - - INPUT: - - R - radius at which to calculate the moment (in natural units) - - n - vR^n - - m - vT^m - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - deriv= None, 'R', or 'phi': calculates derivative of the moment wrt R or phi - - OUTPUT: - - at R (no support for units) - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) - + Calculate the an arbitrary moment of the velocity distribution at R times the surfacmass + + Parameters + ---------- + R: float or Quantity + Galactocentric radius at which to calculate the moment. + n: int + vR^n in the moment + m: int + vT^m in the moment + nsigma : int, optional + number of sigma to integrate the velocities over + romberg : bool, optional + If True, use a romberg integrator (default: False) + deriv : str, optional + Calculates derivative of the moment wrt R or phi (default: None) + + Returns + ------- + float or Quantity + at R (no support for units) + + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU) """ use_physical = kwargs.pop("use_physical", True) ro = kwargs.pop("ro", None) @@ -1155,40 +1083,29 @@ def _vmomentsurfacemass( @physical_conversion("frequency-kmskpc", pop=True) def oortA(self, R, romberg=False, nsigma=None, phi=0.0): """ - - NAME: - - oortA - - PURPOSE: - - calculate the Oort function A - - INPUT: - - R - radius at which to calculate A (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - Oort A at R - - HISTORY: - - 2011-04-19 - Written - Bovy (NYU) - - BUGS: - - could be made more efficient, e.g., surfacemass is calculated multiple times - - """ + Calculate the Oort function A. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate A. + phi : float, optional + Azimuth (default: 0.0). + nsigma : int, optional + Number of sigma to integrate the velocities over. + romberg : bool, optional + If True, use a romberg integrator (default: False). + + Returns + ------- + float or Quantity + Oort A at R. + + Notes + ----- + - 2011-04-19 - Written - Bovy (NYU) + """ + # Could be made more efficient, e.g., surfacemass is calculated multiple times. # 2A= meanvphi/R-dmeanvR/R/dphi-dmeanvphi/dR meanvphi = self.meanvT( R, romberg=romberg, nsigma=nsigma, phi=phi, use_physical=False @@ -1210,39 +1127,29 @@ def oortA(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("frequency-kmskpc", pop=True) def oortB(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - oortB - - PURPOSE: - - calculate the Oort function B - - INPUT: - - R - radius at which to calculate B (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - Oort B at R - - HISTORY: - - 2011-04-19 - Written - Bovy (NYU) - - BUGS: - - could be made more efficient, e.g., surfacemass is calculated multiple times - - """ + Calculate the Oort function B. + + Parameters + ---------- + R : float + Radius at which to calculate B (can be Quantity). + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : float, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth angle (in radians) at which to calculate B. + + Returns + ------- + float or Quantity + Oort B at R. + + Notes + ----- + - 2011-04-19 - Written - Bovy (NYU). + """ + # Could be made more efficient, e.g., surfacemass is calculated multiple times. # 2B= -meanvphi/R+dmeanvR/R/dphi-dmeanvphi/dR meanvphi = self.meanvT( R, romberg=romberg, nsigma=nsigma, phi=phi, use_physical=False @@ -1264,40 +1171,30 @@ def oortB(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("frequency-kmskpc", pop=True) def oortC(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - oortC - - PURPOSE: - - calculate the Oort function C - - INPUT: - - R - radius at which to calculate C (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - Oort C at R - - HISTORY: - - 2011-04-19 - Written - Bovy (NYU) - - BUGS: - - could be made more efficient, e.g., surfacemass is calculated multiple times - we know this is zero, but it is calculated anyway (bug or feature?) - - """ + Calculate the Oort function C. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate C (can be Quantity). + nsigma : int, optional + Number of sigma to integrate the velocities over. + romberg : bool, optional + If True, use a romberg integrator (default: False). + phi : float, optional + Azimuth (default: 0.0). + + Returns + ------- + float or Quantity + Oort C at R. + + Notes + ----- + - 2011-04-19 - Written - Bovy (NYU) + """ + # - Could be made more efficient, e.g., surfacemass is calculated multiple times. + # - We know this is zero, but it is calculated anyway (bug or feature?). # 2C= -meanvR/R-dmeanvphi/R/dphi+dmeanvR/dR meanvr = self.meanvR( R, romberg=romberg, nsigma=nsigma, phi=phi, use_physical=False @@ -1318,40 +1215,30 @@ def oortC(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("frequency-kmskpc", pop=True) def oortK(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - oortK - - PURPOSE: - - calculate the Oort function K - - INPUT: - - R - radius at which to calculate K (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - Oort K at R - - HISTORY: - - 2011-04-19 - Written - Bovy (NYU) - - BUGS: - - could be made more efficient, e.g., surfacemass is calculated multiple times - we know this is zero, but it is calculated anyway (bug or feature?) - - """ + Calculate the Oort function K. + + Parameters + ---------- + R : float + Radius at which to calculate K (can be Quantity). + phi : float, optional + Azimuth angle (in radians) at which to calculate K. + nsigma : int, optional + Number of sigma to integrate the velocities over. + romberg : bool, optional + If True, use a romberg integrator (default: False). + + Returns + ------- + float or Quantity + Oort K at R. + + Notes + ----- + - 2011-04-19 - Written - Bovy (NYU) + """ + # - Could be made more efficient, e.g., surfacemass is calculated multiple times. + # - We know this is zero, but it is calculated anyway (bug or feature?). # 2K= meanvR/R+dmeanvphi/R/dphi+dmeanvR/dR meanvr = self.meanvR( R, romberg=romberg, nsigma=nsigma, phi=phi, use_physical=False @@ -1372,35 +1259,29 @@ def oortK(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("velocity2", pop=True) def sigma2(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - sigma2 - - PURPOSE: - - calculate sigma_R^2 at R by marginalizing over velocity - - INPUT: + Calculate sigma_R^2 at R by marginalizing over velocity. - R - radius at which to calculate sigma_R^2 density (can be Quantity) + Parameters + ---------- + R : float + Radius at which to calculate sigma_R^2 density. + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : int, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth angle at which to calculate sigma_R^2 density. - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - sigma_R^2 at R - - HISTORY: - - 2010-03-XX - Written - Bovy (NYU) + Returns + ------- + float or Quantity + Sigma_R^2 at R. + Notes + ----- + - 2010-03-XX - Written - Bovy (NYU) """ + return self.sigma2surfacemass( R, romberg, nsigma, use_physical=False ) / self.surfacemass(R, romberg, nsigma, use_physical=False) @@ -1409,34 +1290,27 @@ def sigma2(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("velocity2", pop=True) def sigmaT2(self, R, romberg=False, nsigma=None, phi=0.0): """ + Calculate sigma_T^2 at R by marginalizing over velocity - NAME: - - sigmaT2 - - PURPOSE: - - calculate sigma_T^2 at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate sigma_T^2 (can be Quantity) + Parameters + ---------- + R : float + Radius at which to calculate sigma_T^2 (can be Quantity) + romberg : bool, optional + If True, use a romberg integrator (default: False) + nsigma : int, optional + Number of sigma to integrate the velocities over + phi : float, optional + Azimuth (default: 0.0) - OPTIONAL INPUT: + Returns + ------- + float or Quantity + Sigma_T^2 at R - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - sigma_T^2 at R - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU) """ surfmass = self.surfacemass( @@ -1452,33 +1326,27 @@ def sigmaT2(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("velocity2", pop=True) def sigmaR2(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - sigmaR2 (duplicate of sigma2 for consistency) - - PURPOSE: - - calculate sigma_R^2 at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate sigma_R^2 (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over + Calculate sigma_R^2 at R by marginalizing over velocity. - KEYWORDS: + Parameters + ---------- + R : float + Radius at which to calculate sigma_R^2. + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : int, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth (default: 0.0). - romberg - if True, use a romberg integrator (default: False) + Returns + ------- + float or Quantity + Sigma_R^2 at R. - OUTPUT: - - sigma_R^2 at R - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU). """ return self.sigma2(R, romberg=romberg, nsigma=nsigma, use_physical=False) @@ -1487,34 +1355,27 @@ def sigmaR2(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("velocity", pop=True) def meanvT(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - meanvT - - PURPOSE: - - calculate at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) + Calculate the mean tangential velocity at a given radius by marginalizing over velocity. - OUTPUT: + Parameters + ---------- + R : float + Radius at which to calculate the mean tangential velocity. + romberg : bool, optional + If True, use a Romberg integrator. Default is False. + nsigma : float, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth angle at which to calculate the mean tangential velocity. - at R - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) + Returns + ------- + float or Quantity + The mean tangential velocity at the given radius. + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU) """ return self._vmomentsurfacemass( R, 0, 1, romberg=romberg, nsigma=nsigma @@ -1524,35 +1385,29 @@ def meanvT(self, R, romberg=False, nsigma=None, phi=0.0): @physical_conversion("velocity", pop=True) def meanvR(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - meanvR - - PURPOSE: - - calculate at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - at R + Calculate at R by marginalizing over velocity. - HISTORY: + Parameters + ---------- + R : float + Radius at which to calculate . + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : float, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth angle at which to calculate . - 2011-03-30 - Written - Bovy (NYU) + Returns + ------- + float or Quantity + at R. + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU). """ + return self._vmomentsurfacemass( R, 1, 0, romberg=romberg, nsigma=nsigma ) / self.surfacemass(R, romberg=romberg, nsigma=nsigma, use_physical=False) @@ -1560,34 +1415,27 @@ def meanvR(self, R, romberg=False, nsigma=None, phi=0.0): @potential_physical_input def skewvT(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - skewvT - - PURPOSE: - - calculate skew in vT at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - skewvT + Calculate skew in vT at R by marginalizing over velocity - HISTORY: + Parameters + ---------- + R : float + Radius at which to calculate + romberg : bool, optional + If True, use a romberg integrator (default: False) + nsigma : float, optional + Number of sigma to integrate the velocities over + phi : float, optional + Azimuth (default: 0.0) - 2011-12-07 - Written - Bovy (NYU) + Returns + ------- + float + Skew in vT + Notes + ----- + - 2011-12-07 - Written - Bovy (NYU) """ surfmass = self.surfacemass( R, romberg=romberg, nsigma=nsigma, use_physical=False @@ -1607,35 +1455,29 @@ def skewvT(self, R, romberg=False, nsigma=None, phi=0.0): @potential_physical_input def skewvR(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: + Calculate skew in vR at R by marginalizing over velocity. - skewvR + Parameters + ---------- + R : float or Quantity + Radius at which to calculate . + romberg : bool, optional + If True, use a romberg integrator (default: False). + nsigma : float, optional + Number of sigma to integrate the velocities over. + phi : float, optional + Azimuth (in radians) at which to calculate the skew in vR. - PURPOSE: - - calculate skew in vR at R by marginalizing over velocity - - INPUT: - - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - skewvR - - HISTORY: - - 2011-12-07 - Written - Bovy (NYU) + Returns + ------- + float + Skew in vR. + Notes + ----- + - 2011-12-07 - Written - Bovy (NYU). """ + surfmass = self.surfacemass( R, romberg=romberg, nsigma=nsigma, use_physical=False ) @@ -1654,33 +1496,27 @@ def skewvR(self, R, romberg=False, nsigma=None, phi=0.0): @potential_physical_input def kurtosisvT(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - kurtosisvT + Calculate excess kurtosis in vT at R by marginalizing over velocity - PURPOSE: + Parameters + ---------- + R : float or Quantity + Radius at which to calculate + romberg : bool, optional + If True, use a romberg integrator (default: False) + nsigma : float, optional + Number of sigma to integrate the velocities over + phi : float, optional + (default: 0.0) - calculate excess kurtosis in vT at R by marginalizing over velocity + Returns + ------- + float + kurtosisvT - INPUT: - - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - kurtosisvT - - HISTORY: - - 2011-12-07 - Written - Bovy (NYU) + Notes + ----- + - 2011-12-07 - Written - Bovy (NYU) """ surfmass = self.surfacemass( @@ -1706,34 +1542,27 @@ def kurtosisvT(self, R, romberg=False, nsigma=None, phi=0.0): @potential_physical_input def kurtosisvR(self, R, romberg=False, nsigma=None, phi=0.0): """ - NAME: - - kurtosisvR - - PURPOSE: - - calculate excess kurtosis in vR at R by marginalizing over velocity + Calculate excess kurtosis in vR at R by marginalizing over velocity - INPUT: + Parameters + ---------- + R : float or Quantity + Radius at which to calculate + romberg : bool, optional + If True, use a romberg integrator (default: False) + nsigma : float, optional + Number of sigma to integrate the velocities over + phi : float or Quantity, optional + Azimuth (default: 0.0) - R - radius at which to calculate (can be Quantity) - - OPTIONAL INPUT: - - nsigma - number of sigma to integrate the velocities over - - KEYWORDS: - - romberg - if True, use a romberg integrator (default: False) - - OUTPUT: - - kurtosisvR - - HISTORY: - - 2011-12-07 - Written - Bovy (NYU) + Returns + ------- + float + KurtosisvR + Notes + ----- + - 2011-12-07 - Written - Bovy (NYU) """ surfmass = self.surfacemass( R, romberg=romberg, nsigma=nsigma, use_physical=False @@ -1757,18 +1586,29 @@ def kurtosisvR(self, R, romberg=False, nsigma=None, phi=0.0): def _ELtowRRapRperi(self, E, L): """ - NAME: - _ELtowRRapRperi - PURPOSE: - calculate the radial frequency based on E,L, also return rap and - rperi - INPUT: - E - energy - L - angular momentum - OUTPUT: - wR(E.L) - HISTORY: - 2010-07-11 - Written - Bovy (NYU) + Calculate the radial frequency based on energy and angular momentum and return the pericenter and apocenter radii. + + Parameters + ---------- + E : float + Energy. + L : float + Angular momentum. + + Returns + ------- + tuple + Tuple containing: + - wR(E.L) : float + Radial frequency. + - rap : float + Apocenter radius. + - rperi : float + Pericenter radius. + + Notes + ----- + - 2010-07-11 - Written - Bovy (NYU) """ if self._beta == 0.0: xE = numpy.exp(E - 0.5) @@ -1800,124 +1640,241 @@ def sample( maxd=None, target=True, ): - r""" - NAME: + """ + Sample n*nphi points from this disk DF. + + Parameters + ---------- + n : int, optional + Number of desired samples. Default is 1. + rrange : list, optional + If you only want samples in this rrange, set this keyword (only works when asking for an (RZ)Orbit). + returnROrbit : bool, optional + If True, return a planarROrbit instance: [R,vR,vT] (default). + returnOrbit : bool, optional + If True, return a planarOrbit instance (including phi). + nphi : float, optional + Number of azimuths to sample for each E,L. + los : float, optional + Line of sight sampling along this line of sight. + losdeg : bool, optional + If True, los is in degrees (default). + nsigma : int, optional + Number of sigma to rejection-sample on. + maxd : float, optional + Maximum distance to consider (for the rejection sampling). + target : bool, optional + If True, use target surface mass and sigma2 profiles (default). + + Returns + ------- + list + n*nphi list of [[E,Lz],...] or list of planar(R)Orbits. + CAUTION: lists of EL need to be post-processed to account for the + \\kappa/\\omega_R discrepancy + + Notes + ----- + - 2010-07-10 - Started - Bovy (NYU) + + """ + raise NotImplementedError("'sample' method for this disk df is not implemented") - sample + def _estimatemeanvR(self, R, phi=0.0, log=False): + """ + Quickly estimate the mean radial velocity at a given radius R. - PURPOSE: + Parameters + ---------- + R : float + Radius at which to evaluate (/ro). + phi : float, optional + Azimuth angle (not used). + log : bool, optional + If True, return the logarithm of the target Sigma_R^2(R). - sample n*nphi points from this DF + Returns + ------- + float + The target Sigma_R^2(R). - INPUT: + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) - n - number of desired sample (specifying this rather than calling this routine n times is more efficient) + """ + return 0.0 - rrange - if you only want samples in this rrange, set this keyword (only works when asking for an (RZ)Orbit) (can be Quantity) + def _estimatemeanvR(self, R, phi=0.0, log=False): + """ + Quickly estimate the mean radial velocity at a given radius. - returnROrbit - if True, return a planarROrbit instance: - [R,vR,vT] (default) + Parameters + ---------- + R : float + Radius at which to evaluate the mean radial velocity (/ro). + phi : float, optional + Azimuth angle (not used). + log : bool, optional + If True, return the logarithm of the target Sigma_R^2(R) (default: False). - returnOrbit - if True, return a planarOrbit instance (including phi) + Returns + ------- + float + The estimated mean radial velocity. - nphi - number of azimuths to sample for each E,L + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) - los= line of sight sampling along this line of sight (can be Quantity) + """ + return 0.0 - losdeg= los in degrees? (default=True) + def _estimatemeanvR(self, R, phi=0.0, log=False): + """ + Quickly estimate the mean radial velocity at a given radius. - target= if True, use target surface mass and sigma2 profiles (default=True) + Parameters + ---------- + R : float + Radius at which to evaluate the mean radial velocity. Given in natural units (/ro). + phi : float, optional + Azimuth angle. Not used. + log : bool, optional + If True, return the logarithm of the target Sigma_R^2(R). Default is False. - nsigma= number of sigma to rejection-sample on + Returns + ------- + float + The estimated mean radial velocity. - maxd= maximum distance to consider (for the rejection sampling) + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) - OUTPUT: + """ + return 0.0 - n*nphi list of [[E,Lz],...] or list of planar(R)Orbits + def _estimatemeanvR(self, R, phi=0.0, log=False): + """ + Quickly estimate the mean radial velocity at a given radius R. - CAUTION: lists of EL need to be post-processed to account for the - \kappa/\omega_R discrepancy + Parameters + ---------- + R : float + Radius at which to evaluate the mean radial velocity. [galpy units] + phi : float, optional + Azimuth angle. Not used. (default: 0.0) + log : bool, optional + If True, return the logarithm of the target Sigma_R^2(R). (default: False) - HISTORY: + Returns + ------- + float + The estimated mean radial velocity at the given radius R. - 2010-07-10 - Started - Bovy (NYU) + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) """ - raise NotImplementedError("'sample' method for this disk df is not implemented") + return 0.0 def _estimatemeanvR(self, R, phi=0.0, log=False): """ - NAME: - _estimatemeanvR - PURPOSE: - quickly estimate meanvR (useful in evolveddiskdf where we - need an estimate of this but we do not want to spend too - much time on it) - INPUT: - R - radius at which to evaluate (/ro) - phi= azimuth (not used) - OUTPUT: - target Sigma_R^2(R) - log - if True, return the log (default: False) - HISTORY: - 2010-03-28 - Written - Bovy (NYU) + Quickly estimate the mean radial velocity at a given radius. + + Parameters + ---------- + R : float + Radius at which to evaluate the mean radial velocity (/ro). + phi : float, optional + Azimuth angle (not used). + log : bool, optional + If True, return the logarithm of the target Sigma_R^2(R) (default: False). + + Returns + ------- + float + The estimated mean radial velocity. + + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) """ return 0.0 def _estimatemeanvT(self, R, phi=0.0, log=False): """ - NAME: - _estimatemeanvT - PURPOSE: - quickly estimate meanvR (useful in evolveddiskdf where we - need an estimate of this but we do not want to spend too - much time on it) - INPUT: - R - radius at which to evaluate (/ro) - phi= azimuth (not used) - OUTPUT: - target Sigma_R^2(R) - HISTORY: - 2010-03-28 - Written - Bovy (NYU) + Quickly estimate the mean tangential velocity at a given radius. + + Parameters + ---------- + R : float + Radius at which to evaluate (/ro). + phi : float, optional + Azimuth angle (not used). + log : bool, optional + If True, return the logarithm of the estimate. + + Returns + ------- + float + The estimated mean tangential velocity. + + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) + """ return R**self._beta - self.asymmetricdrift(R, use_physical=False) def _estimateSigmaR2(self, R, phi=0.0, log=False): """ - NAME: - _estimateSigmaR2 - PURPOSE: - quickly estimate SigmaR2 (useful in evolveddiskdf where we - need an estimate of this but we do not want to spend too - much time on it) - INPUT: - R - radius at which to evaluate (/ro) - phi= azimuth (not used) - OUTPUT: - target Sigma_R^2(R) - log - if True, return the log (default: False) - HISTORY: - 2010-03-28 - Written - Bovy (NYU) + Quickly estimate SigmaR2. + + Parameters + ---------- + R : float + Radius at which to evaluate (/ro). + phi : float, optional + Azimuth (not used). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + float + Target Sigma_R^2(R). + + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) """ return self.targetSigma2(R, log=log, use_physical=False) def _estimateSigmaT2(self, R, phi=0.0, log=False): """ - NAME: - _estimateSigmaT2 - PURPOSE: - quickly estimate SigmaT2 (useful in evolveddiskdf where we - need an estimate of this but we do not want to spend too - much time on it) - INPUT: - R - radius at which to evaluate (/ro) - phi= azimuth (not used) - OUTPUT: - target Sigma_R^2(R) - log - if True, return the log (default: False) - HISTORY: - 2010-03-28 - Written - Bovy (NYU) + Quickly estimate SigmaT2. + + Parameters + ---------- + R : float + Radius at which to evaluate (/ro). + phi : float, optional + Azimuth (not used). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + float + Target Sigma_R^2(R). + + Notes + ----- + - 2010-03-28 - Written - Bovy (NYU) + """ if log: return self.targetSigma2(R, log=log, use_physical=False) - 2.0 * numpylog( @@ -1941,43 +1898,32 @@ def __init__( **kwargs ): """ - NAME: - __init__ - PURPOSE: - Initialize a Dehnen 'new' DF - INPUT: - surfaceSigma - instance or class name of the target - surface density and sigma_R profile - (default: both exponential) - profileParams - parameters of the surface and sigma_R profile: - (xD,xS,Sro) where - - xD - disk surface mass scalelength (can be Quantity) - - xS - disk velocity dispersion scalelength (can be Quantity) - - Sro - disk velocity dispersion at Ro (can be Quantity) - - Directly given to the 'surfaceSigmaProfile class, so - could be anything that class takes - - beta - power-law index of the rotation curve - - correct - if True, correct the DF - - ro= distance from vantage point to GC (kpc; can be Quantity) - - vo= circular velocity at ro (km/s; can be Quantity) - - +DFcorrection kwargs (except for those already specified) - - OUTPUT: - - instance - - HISTORY: - - 2010-03-10 - Written - Bovy (NYU) + Initialize a Dehnen 'new' DF. + + Parameters + ---------- + surfaceSigma : instance or class name of the target surface density and sigma_R profile, optional + Default: both exponential. + profileParams : tuple, optional + Parameters of the surface and sigma_R profile: (xD,xS,Sro) where: + * xD - disk surface mass scalelength (can be Quantity) + * xS - disk velocity dispersion scalelength (can be Quantity) + * Sro - disk velocity dispersion at Ro (can be Quantity) + Directly given to the 'surfaceSigmaProfile class, so could be anything that class takes. + beta : float, optional + Power-law index of the rotation curve. + correct : bool, optional + If True, correct the DF. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + **kwargs: dict, optional + DFcorrection kwargs (except for those already specified). + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) """ return diskdf.__init__( @@ -1992,18 +1938,28 @@ def __init__( def eval(self, E, L, logSigmaR=0.0, logsigmaR2=0.0): """ - NAME: - eval - PURPOSE: - evaluate the distribution function - INPUT: - E - energy (can be Quantity) - L - angular momentum (can be Quantity) - OUTPUT: + Evaluate the distribution function. + + Parameters + ---------- + E : float or Quantity + Energy. + L : float or Quantity + Angular momentum. + logSigmaR : float, optional + Logarithm of the radial velocity dispersion. + logsigmaR2 : float, optional + Logarithm of the square of the radial velocity dispersion. + + Returns + ------- + float DF(E,L) - HISTORY: - 2010-03-10 - Written - Bovy (NYU) - 2010-03-28 - Moved to dehnenDF - Bovy (NYU) + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU). + - 2010-03-28 - Moved to dehnenDF - Bovy (NYU). """ if _PROFILE: # pragma: no cover import time @@ -2085,32 +2041,51 @@ def sample( maxd=None, **kwargs ): - r""" - NAME: - sample - PURPOSE: - sample n*nphi points from this DF - INPUT: - n - number of desired sample (specifying this rather than calling - this routine n times is more efficient) - rrange - if you only want samples in this rrange, set this keyword - (only works when asking for an (RZ)Orbit - returnROrbit - if True, return a planarROrbit instance: - [R,vR,vT] (default) - returnOrbit - if True, return a planarOrbit instance (including phi) - nphi - number of azimuths to sample for each E,L - los= if set, sample along this line of sight (deg) (assumes that the Sun is located at R=1,phi=0) - losdeg= if False, los is in radians (default=True) - targetSurfmass, targetSigma2= if True, use target surface mass and sigma2 profiles, respectively (there is not much point to doing the latter) - (default=True) - nsigma= number of sigma to rejection-sample on - maxd= maximum distance to consider (for the rejection sampling) - OUTPUT: - n*nphi list of [[E,Lz],...] or list of planar(R)Orbits - CAUTION: lists of EL need to be post-processed to account for the - \kappa/\omega_R discrepancy; EL not returned in physical units - HISTORY: - 2010-07-10 - Started - Bovy (NYU) + """ + Sample n*nphi points from this DF. + + Parameters + ---------- + n : int, optional + Number of desired samples (specifying this rather than calling + this routine n times is more efficient). Default is 1. + rrange : list or tuple, optional + If you only want samples in this rrange, set this keyword + (only works when asking for an (RZ)Orbit). Default is None. + returnROrbit : bool, optional + If True, return a planarROrbit instance: [R,vR,vT] (default). + Default is True. + returnOrbit : bool, optional + If True, return a planarOrbit instance (including phi). + Default is False. + nphi : float, optional + Number of azimuths to sample for each E,L. Default is 1.0. + los : float or Quantity, optional + If set, sample along this line of sight (deg) (assumes that the Sun is located at R=1,phi=0). + Default is None. + losdeg : bool, optional + If False, los is in radians (default=True). Default is True. + nsigma : int, optional + Number of sigma to rejection-sample on. Default is None. + targetSurfmass : bool, optional + If True, use target surface mass profile. Default is True. + targetSigma2 : bool, optional + If True, use target sigma2 profile. Default is True. + maxd : float or Quantity, optional + Maximum distance to consider (for the rejection sampling). Default is None. + **kwargs : dict, optional + Additional keyword arguments. + + Returns + ------- + out : list + n*nphi list of [[E,Lz],...] or list of planar(R)Orbits. + CAUTION: lists of EL need to be post-processed to account for the + \\kappa/\\omega_R discrepancy; EL not returned in physical units. + + Notes + ----- + - 2010-07-10 - Started - Bovy (NYU) """ if not los is None: return self.sampleLOS( @@ -2374,43 +2349,32 @@ def __init__( **kwargs ): """ - NAME: - __init__ - PURPOSE: - Initialize a Shu DF - INPUT: - surfaceSigma - instance or class name of the target - surface density and sigma_R profile - (default: both exponential) - profileParams - parameters of the surface and sigma_R profile: - (xD,xS,Sro) where - - xD - disk surface mass scalelength (can be Quantity) - - xS - disk velocity dispersion scalelength (can be Quantity) - - Sro - disk velocity dispersion at Ro (can be Quantity) - - Directly given to the 'surfaceSigmaProfile class, so - could be anything that class takes - - beta - power-law index of the rotation curve - - correct - if True, correct the DF - - ro= distance from vantage point to GC (kpc; can be Quantity) - - vo= circular velocity at ro (km/s; can be Quantity) - - +DFcorrection kwargs (except for those already specified) - - OUTPUT: - - instance - - HISTORY: - - 2010-05-09 - Written - Bovy (NYU) + Initialize a Shu DF. + + Parameters + ---------- + surfaceSigma : instance or class name of the target surface density and sigma_R profile, optional + Default: both exponential. + profileParams : tuple, optional + Parameters of the surface and sigma_R profile: (xD,xS,Sro) where + * xD - disk surface mass scalelength (can be Quantity) + * xS - disk velocity dispersion scalelength (can be Quantity) + * Sro - disk velocity dispersion at Ro (can be Quantity) + Directly given to the 'surfaceSigmaProfile class, so could be anything that class takes. + beta : float, optional + Power-law index of the rotation curve. + correct : bool, optional + If True, correct the DF. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + **kwargs: dict, optional + DFcorrection kwargs (except for those already specified). + + Notes + ----- + - 2010-05-09 - Written - Bovy (NYU) """ return diskdf.__init__( @@ -2425,17 +2389,27 @@ def __init__( def eval(self, E, L, logSigmaR=0.0, logsigmaR2=0.0): """ - NAME: - eval - PURPOSE: - evaluate the distribution function - INPUT: - E - energy (/vo^2) - L - angular momentun (/ro/vo) - OUTPUT: - DF(E,L) - HISTORY: - 2010-05-09 - Written - Bovy (NYU) + Evaluate the distribution function. + + Parameters + ---------- + E : float + Energy (/vo^2). + L : float + Angular momentum (/ro/vo). + logSigmaR : float, optional + Logarithm of the radial velocity dispersion squared. + logsigmaR2 : float, optional + Logarithm of the radial velocity dispersion squared. + + Returns + ------- + float + DF(E,L). + + Notes + ----- + - 2010-05-09 - Written - Bovy (NYU) """ E = conversion.parse_energy(E, vo=self._vo) L = conversion.parse_angmom(L, ro=self._ro, vo=self._vo) @@ -2484,32 +2458,51 @@ def sample( targetSigma2=True, **kwargs ): - r""" - NAME: - sample - PURPOSE: - sample n*nphi points from this DF - INPUT: - n - number of desired sample (specifying this rather than calling - this routine n times is more efficient) - rrange - if you only want samples in this rrange, set this keyword - (only works when asking for an (RZ)Orbit - returnROrbit - if True, return a planarROrbit instance: - [R,vR,vT] (default) - returnOrbit - if True, return a planarOrbit instance (including phi) - nphi - number of azimuths to sample for each E,L - los= if set, sample along this line of sight (deg) (assumes that the Sun is located at R=1,phi=0) - losdeg= if False, los is in radians (default=True) - targetSurfmass, targetSigma2= if True, use target surface mass and sigma2 profiles, respectively (there is not much point to doing the latter) - (default=True) - nsigma= number of sigma to rejection-sample on - maxd= maximum distance to consider (for the rejection sampling) - OUTPUT: - n*nphi list of [[E,Lz],...] or list of planar(R)Orbits - CAUTION: lists of EL need to be post-processed to account for the - \kappa/\omega_R discrepancy - HISTORY: - 2010-07-10 - Started - Bovy (NYU) + """ + Sample n*nphi points from this DF. + + Parameters + ---------- + n : int, optional + Number of desired samples (specifying this rather than calling + this routine n times is more efficient). Default is 1. + rrange : list or tuple, optional + If you only want samples in this rrange, set this keyword + (only works when asking for an (RZ)Orbit). Default is None. + returnROrbit : bool, optional + If True, return a planarROrbit instance: [R,vR,vT] (default). + Default is True. + returnOrbit : bool, optional + If True, return a planarOrbit instance (including phi). + Default is False. + nphi : float, optional + Number of azimuths to sample for each E,L. Default is 1.0. + los : float or Quantity, optional + If set, sample along this line of sight (deg) (assumes that the Sun is located at R=1,phi=0). + Default is None. + losdeg : bool, optional + If False, los is in radians (default=True). Default is True. + nsigma : int, optional + Number of sigma to rejection-sample on. Default is None. + targetSurfmass : bool, optional + If True, use target surface mass profile. Default is True. + targetSigma2 : bool, optional + If True, use target sigma2 profile. Default is True. + maxd : float or Quantity, optional + Maximum distance to consider (for the rejection sampling). Default is None. + **kwargs : dict, optional + Additional keyword arguments. + + Returns + ------- + out : list + n*nphi list of [[E,Lz],...] or list of planar(R)Orbits. + CAUTION: lists of EL need to be post-processed to account for the + \\kappa/\\omega_R discrepancy; EL not returned in physical units. + + Notes + ----- + - 2010-07-10 - Started - Bovy (NYU) """ if not los is None: return self.sampleLOS( @@ -2704,43 +2697,32 @@ def __init__( **kwargs ): """ - NAME: - __init__ - PURPOSE: - Initialize a Schwarzschild DF - INPUT: - surfaceSigma - instance or class name of the target - surface density and sigma_R profile - (default: both exponential) - profileParams - parameters of the surface and sigma_R profile: - (xD,xS,Sro) where - - xD - disk surface mass scalelength (can be Quantity) - - xS - disk velocity dispersion scalelength (can be Quantity) - - Sro - disk velocity dispersion at Ro (can be Quantity) - - Directly given to the 'surfaceSigmaProfile class, so - could be anything that class takes - - beta - power-law index of the rotation curve - - correct - if True, correct the DF - - ro= distance from vantage point to GC (kpc; can be Quantity) - - vo= circular velocity at ro (km/s; can be Quantity) - - +DFcorrection kwargs (except for those already specified) - - OUTPUT: - - instance - - HISTORY: - - 2017-09-17 - Written - Bovy (UofT) + Initialize a Schwarzschild DF. + + Parameters + ---------- + surfaceSigma : instance or class name of the target surface density and sigma_R profile, optional + (default: both exponential) + profileParams : tuple, optional + Parameters of the surface and sigma_R profile: (xD,xS,Sro) where + * xD - disk surface mass scalelength (can be Quantity) + * xS - disk velocity dispersion scalelength (can be Quantity) + * Sro - disk velocity dispersion at Ro (can be Quantity) + Directly given to the 'surfaceSigmaProfile class, so could be anything that class takes + beta : float, optional + Power-law index of the rotation curve + correct : bool, optional + If True, correct the DF + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + **kwargs : dict, optional + DFcorrection kwargs (except for those already specified) + + Notes + ----- + - 2017-09-17 - Written - Bovy (UofT) """ # Schwarzschild == Shu w/ energy computed in epicycle approx. @@ -2824,16 +2806,33 @@ def _oned_intFunc(x, twodfunc, gfun, hfun, tol, args): def bovy_dblquad(func, a, b, gfun, hfun, args=(), tol=1.48e-08): """ - NAME: - bovy_dblquad - PURPOSE: - like scipy.integrate's dblquad, but using Romberg integration for the one-d integrals and using tol - INPUT: - same as scipy.integrate.dblquad except for tol and epsrel,epsabs - OUTPUT: - value - HISTORY: - 2010-03-11 - Written - Bpvy (NYU) + Compute a double integral using Romberg integration for the one-dimensional integrals and a specified tolerance. + + Parameters + ---------- + func : callable + Function of two variables to integrate. + a : float + Lower limit of integration in the outer integral. + b : float + Upper limit of integration in the outer integral. + gfun : callable + Function of one variable that returns the lower limit of integration in the inner integral for a given value of the outer variable. + hfun : callable + Function of one variable that returns the upper limit of integration in the inner integral for a given value of the outer variable. + args : tuple, optional + Extra arguments to pass to the integrand function. + tol : float, optional + Desired absolute tolerance. + + Returns + ------- + float + The value of the double integral. + + Notes + ----- + - 2010-03-11 - Written - Bpvy (NYU) """ return integrate.romberg( _oned_intFunc, a, b, args=(func, gfun, hfun, tol, args), tol=tol @@ -2846,25 +2845,33 @@ class DFcorrection: def __init__(self, **kwargs): """ - NAME: - __init__ - PURPOSE: - initialize the corrections: set them, load them, or calculate - and save them - OPTIONAL INPUTS: - corrections - if Set, these are the corrections and they should - be used as such - npoints - number of points from 0 to Rmax - rmax - correct up to this radius (/ro) (default: 5) - savedir - save the corrections in this directory - surfaceSigmaProfile - target surfacemass and sigma_R^2 instance - beta - power-law index of the rotation curve (when calculating) - dftype - classname of the DF - niter - number of iterations to perform to calculate the corrections - interp_k - 'k' keyword to give to InterpolatedUnivariateSpline - OUTPUT: - HISTORY: - 2010-03-10 - Written - Bovy (NYU) + Initialize the corrections: set them, load them, or calculate and save them. + + Parameters + ---------- + corrections : numpy.ndarray, optional + If set, these are the corrections and they should be used as such. + npoints : int, optional + Number of points from 0 to Rmax. + rmax : float, optional + Correct up to this radius (/ro) (default: 5). + savedir : str, optional + Save the corrections in this directory. + surfaceSigmaProfile : object + Target surfacemass and sigma_R^2 instance. + beta : float, optional + Power-law index of the rotation curve (when calculating). + dftype : class, optional + Classname of the DF. + niter : int, optional + Number of iterations to perform to calculate the corrections. + interp_k : str, optional + 'k' keyword to give to InterpolatedUnivariateSpline. + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) + """ if not "surfaceSigmaProfile" in kwargs: raise DFcorrectionError("surfaceSigmaProfile not given") @@ -2944,17 +2951,23 @@ def _createSavefilename(self, niter): def correct(self, R, log=False): """ - NAME: - correct - PURPOSE: - calculate the correction in Sigma and sigma2 at R - INPUT: - R - Galactocentric radius(/ro) - log - if True, return the log of the correction - OUTPUT: - [Sigma correction, sigma2 correction] - HISTORY: - 2010-03-10 - Written - Bovy (NYU) + Calculate the correction in Sigma and sigma2 at R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the log of the correction. + + Returns + ------- + tuple + (Sigma correction, sigma2 correction). + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) """ if isinstance(R, numpy.ndarray): out = numpy.empty((2, len(R))) @@ -3007,17 +3020,21 @@ def correct(self, R, log=False): def derivLogcorrect(self, R): """ - NAME: - derivLogcorrect - PURPOSE: - calculate the derivative of the log of the correction in Sigma - and sigma2 at R - INPUT: - R - Galactocentric radius(/ro) - OUTPUT: - [d log(Sigma correction)/dR, d log(sigma2 correction)/dR] - HISTORY: - 2010-03-10 - Written - Bovy (NYU) + Calculate the derivative of the log of the correction in Sigma and sigma2 at R. + + Parameters + ---------- + R : float + Galactocentric radius(/ro) + + Returns + ------- + numpy.ndarray + [d log(Sigma correction)/dR, d log(sigma2 correction)/dR] + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) """ if R < _RMIN: out = numpy.array([self._surfaceDerivSmallR, self._sigma2DerivSmallR]) @@ -3103,17 +3120,30 @@ def __str__(self): def vRvTRToEL(vR, vT, R, beta, dftype="dehnen"): """ - NAME: - vRvTRToEL - PURPOSE: - calculate the energy and angular momentum - INPUT: - vR - radial velocity - vT - rotational velocity - R - Galactocentric radius - OUTPUT: - HISTORY: - 2010-03-10 - Written - Bovy (NYU) + Calculate the energy and angular momentum. + + Parameters + ---------- + vR : float + Radial velocity. + vT : float + Rotational velocity. + R : float + Galactocentric radius. + beta : float + Parameter that determines the shape of the rotation curve. + dftype : str, optional + Type of disk distribution function. Default is "dehnen". + + Returns + ------- + tuple + Energy and angular momentum. + + Notes + ----- + - 2010-03-10 - Written - Bovy (NYU) + """ if dftype == "schwarzschild": # Compute E in the epicycle approximation @@ -3136,17 +3166,24 @@ def vRvTRToEL(vR, vT, R, beta, dftype="dehnen"): def axipotential(R, beta=0.0): """ - NAME: - axipotential - PURPOSE: - return the axisymmetric potential at R/Ro - INPUT: - R - Galactocentric radius - beta - rotation curve power-law - OUTPUT: - Pot(R)/vo**2. - HISTORY: - 2010-03-01 - Written - Bovy (NYU) + Return the axisymmetric potential at R/Ro. + + Parameters + ---------- + R : float + Galactocentric radius. + beta : float, optional + Rotation curve power-law. + + Returns + ------- + float + Pot(R)/vo**2. + + Notes + ----- + - 2010-03-01 - Written - Bovy (NYU) + """ if beta == 0.0: if numpy.any(R == 0.0): @@ -3162,20 +3199,26 @@ def axipotential(R, beta=0.0): def _ars_hx(x, args): """ - NAME: - _ars_hx - PURPOSE: - h(x) for ARS sampling of the input surfacemass profile - INPUT: - x - R(/ro) - args= (surfaceSigma, dfcorr) - surfaceSigma - surfaceSigmaProfile instance - dfcorr - DFcorrection instance - OUTPUT: - log(x)+log surface(x) + log(correction) - HISTORY: - 2010-07-11 - Written - Bovy (NYU) + h(x) for ARS sampling of the input surfacemass profile + + Parameters + ---------- + x : float + R(/ro) + args : tuple + surfaceSigma - surfaceSigmaProfile instance + dfcorr - DFcorrection instance + + Returns + ------- + float + log(x)+log surface(x) + log(correction) + + Notes + ----- + - 2010-07-11 - Written - Bovy (NYU) """ + surfaceSigma, dfcorr = args if dfcorr is None: return numpylog(x) + surfaceSigma.surfacemass(x, log=True) @@ -3187,19 +3230,24 @@ def _ars_hx(x, args): def _ars_hpx(x, args): """ - NAME: - _ars_hpx - PURPOSE: - h'(x) for ARS sampling of the input surfacemass profile - INPUT: - x - R(/ro) - args= (surfaceSigma, dfcorr) - surfaceSigma - surfaceSigmaProfile instance - dfcorr - DFcorrection instance - OUTPUT: - derivative of log(x)+log surface(x) + log(correction) wrt x - HISTORY: - 2010-07-11 - Written - Bovy (NYU) + h'(x) for ARS sampling of the input surfacemass profile + + Parameters + ---------- + x : float + R(/ro) + args : tuple + surfaceSigma - surfaceSigmaProfile instance + dfcorr - DFcorrection instance + + Returns + ------- + float + derivative of log(x)+log surface(x) + log(correction) wrt x + + Notes + ----- + - 2010-07-11 - Written - Bovy (NYU) """ surfaceSigma, dfcorr = args if dfcorr is None: diff --git a/galpy/df/eddingtondf.py b/galpy/df/eddingtondf.py index 079b9e8a8..097d26368 100644 --- a/galpy/df/eddingtondf.py +++ b/galpy/df/eddingtondf.py @@ -21,33 +21,26 @@ class eddingtondf(isotropicsphericaldf): def __init__(self, pot=None, denspot=None, rmax=1e4, scale=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an isotropic distribution function computed using the Eddington inversion - - INPUT: - - pot= (None) Potential instance or list thereof that represents the gravitational potential (assumed to be spherical) - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale= Characteristic scale radius to aid sampling calculations. Optionaland will also be overridden by value from pot if available. - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2021-02-04 - Written - Bovy (UofT) + Initialize an isotropic distribution function computed using the Eddington inversion. + + Parameters + ---------- + pot : Potential instance or list thereof + Represents the gravitational potential (assumed to be spherical). + denspot : Potential instance or list thereof, optional + Represents the density of the tracers (assumed to be spherical; if None, set equal to pot). + rmax : float or Quantity, optional + Maximum radius to consider. DF is cut off at E = Phi(rmax). + scale : float or Quantity, optional + Characteristic scale radius to aid sampling calculations. Optional and will also be overridden by value from pot if available. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2021-02-04 - Written - Bovy (UofT) """ isotropicsphericaldf.__init__( @@ -90,25 +83,21 @@ def sample(self, R=None, z=None, phi=None, n=1, return_orbit=True, rmin=0.0): def fE(self, E): """ - NAME: - - fE - - PURPOSE - - Calculate the energy portion of a DF computed using the Eddington inversion - - INPUT: - - E - The energy (can be Quantity) - - OUTPUT: + Calculate the energy portion of a DF computed using the Eddington inversion - fE - The value of the energy portion of the DF + Parameters + ---------- + E : float or Quantity + The energy. - HISTORY: + Returns + ------- + fE : ndarray + The value of the energy portion of the DF. - 2021-02-04 - Written - Bovy (UofT) + Notes + ----- + - 2021-02-04 - Written - Bovy (UofT) """ Eint = conversion.parse_energy(E, vo=self._vo) out = numpy.zeros_like(Eint) diff --git a/galpy/df/evolveddiskdf.py b/galpy/df/evolveddiskdf.py index ec6e9923f..ede7eff85 100644 --- a/galpy/df/evolveddiskdf.py +++ b/galpy/df/evolveddiskdf.py @@ -35,30 +35,20 @@ class evolveddiskdf(df): def __init__(self, initdf, pot, to=0.0): """ - NAME: - - __init__ - - PURPOSE: - - initialize - - INPUT: - - initdf - the df at the start of the evolution (at to) (units are transferred) - - pot - potential to integrate orbits in - - to= initial time (time at which initdf is evaluated; orbits are integrated from current t back to to) (can be Quantity) - - OUTPUT: - - instance - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) - + Initialize the evolved disk distribution function. + + Parameters + ---------- + initdf : galpy.df.df instance + The distribution function at the start of the evolution (at to) (units are transferred). + pot : galpy.potential.Potential instance + Potential to integrate orbits in. + to : float or Quantity, optional + Initial time (time at which initdf is evaluated; orbits are integrated from current t back to to). + + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU) """ if initdf._roSet: ro = initdf._ro @@ -76,44 +66,37 @@ def __init__(self, initdf, pot, to=0.0): @physical_conversion("phasespacedensity2d", pop=True) def __call__(self, *args, **kwargs): """ - NAME: - - __call__ - - PURPOSE: - - evaluate the distribution function - - INPUT: - - Orbit instance: - - a) Orbit instance alone: use initial state and t=0 - - b) Orbit instance + t: Orbit instance *NOT* called (i.e., Orbit's initial condition is used, call Orbit yourself), t can be Quantity - - If t is a list of t, DF is returned for each t, times must be in descending order and equally spaced (does not work with marginalize...) - - marginalizeVperp - marginalize over perpendicular velocity (only supported with 1a) above) + nsigma, +scipy.integrate.quad keywords - - marginalizeVlos - marginalize over line-of-sight velocity (only supported with 1a) above) + nsigma, +scipy.integrate.quad keywords - - log= if True, return the log (not for deriv, bc that can be negative) - - integrate_method= method argument of orbit.integrate - - deriv= None, 'R', or 'phi': calculates derivative of the moment wrt R or phi **not with the marginalize options** - - OUTPUT: - - DF(orbit,t) - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) - - 2011-04-15 - Added list of times option - Bovy (NYU) - + Evaluate the distribution function + + Parameters + ---------- + *args : tuple + Either: + 1) Orbit instance alone: use initial state and t=0 + 2) Orbit instance + t: Orbit instance *NOT* called (i.e., Orbit's initial condition is used, call Orbit yourself), t can be Quantity + If t is a list of t, DF is returned for each t, times must be in descending order and equally spaced (does not work with marginalize...) + marginalizeVperp : bool, optional + marginalize over perpendicular velocity (only supported with 1a) for single orbits above) + marginalizeVlos : bool, optional + marginalize over line-of-sight velocity (only supported with 1a) for single orbits above) + integrate_method : str, optional + orbit.integrate method argument + log : bool, optional + if True, return the log (not for deriv, bc that can be negative) + deriv : str, optional + None, 'R', or 'phi': calculates derivative of the moment wrt R or phi **not with the marginalize options** + **kwargs: dict, optional + scipy.integrate.quad keywords + + Returns + ------- + float or numpy.ndarray + value of DF + + Notes + ----- + - 2011-03-30 - Written - Bovy (NYU) + - 2011-04-15 - Added list of times option - Bovy (NYU) """ integrate_method = kwargs.pop("integrate_method", "dopr54_c") # Must match Python fallback for non-C potentials here, bc odeint needs @@ -460,60 +443,54 @@ def vmomentsurfacemass( deriv=None, ): """ - NAME: - - vmomentsurfacemass - - PURPOSE: - - calculate the an arbitrary moment of the velocity distribution at (R,phi) times the surfacmass - - INPUT: - - R - radius at which to calculate the moment (in natural units) - - phi= azimuth (rad unless deg=True) - - n - vR^n - - m - vT^m - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous, but not too generous) - - deg= azimuth is in degree (default=False) - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid; if this was created for a list of times, moments are calculated for each time - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - print_progress= if True, print progress updates - - integrate_method= orbit.integrate method argument - - deriv= None, 'R', or 'phi': calculates derivative of the moment wrt R or phi **onnly with grid options** - - OUTPUT: - - at R,phi (no support for units) - - COMMENT: - - grid-based calculation is the only one that is heavily tested (although the test suite also tests the direct calculation) - - HISTORY: - - 2011-03-30 - Written - Bovy (NYU) - + Calculate the an arbitrary moment of the velocity distribution at (R,phi) times the surfacmass + + Parameters + ---------- + R : float + Radius at which to calculate the moment (in natural units). + phi : float, optional + Azimuth (rad unless deg=True). + n : int + vR^n. + m : int + vT^m. + t : float, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced). + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous, but not too generous). + deg : bool, optional + Azimuth is in degree (default=False). + epsrel : float, optional + scipy.integrate keyword (the integration calculates the ratio of this vmoment to that of the initial DF). + epsabs : float, optional + scipy.integrate keyword (the integration calculates the ratio of this vmoment to that of the initial DF). + grid : bool, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid; if this was created for a list of times, moments are calculated for each time. + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101). + returnGrid : bool, optional + If True, return the grid object (default=False). + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False). + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid. + print_progress : bool, optional + If True, print progress updates. + integrate_method : str, optional + orbit.integrate method argument. + deriv : str, optional + None, 'R', or 'phi': calculates derivative of the moment wrt R or phi **onnly with grid options**. + + Returns + ------- + float or numpy.ndarray + at R,phi (no support for units). + + Notes + ----- + - grid-based calculation is the only one that is heavily tested (although the test suite also tests the direct calculation) + - 2011-03-30 - Written - Bovy (NYU) """ # if we have already precalculated a grid, use that if not grid is None and isinstance(grid, evolveddiskdfGrid): @@ -654,50 +631,53 @@ def vertexdev( integrate_method="dopr54_c", ): """ - NAME: - - vertexdev - - PURPOSE: - - calculate the vertex deviation of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - sigmaR2, sigmaT2, sigmaRT= if set the vertex deviation is simply calculated using these - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - vertex deviation in rad - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the vertex deviation of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + radius at which to calculate the moment (can be Quantity) + t : float, optional + time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + nsigma : float, optional + number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + phi : float, optional + azimuth (rad unless deg=True; can be Quantity), by default 0.0 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + number of points to use for the grid in 1D (default=101), by default 101 + returnGrid : bool, optional + if True, return the grid object (default=False), by default False + sigmaR2 : float, optional + if set the vertex deviation is simply calculated using these, by default None + sigmaT2 : float, optional + if set the vertex deviation is simply calculated using these, by default None + sigmaRT : float, optional + if set the vertex deviation is simply calculated using these, by default None + surfacemass : None, optional + Not used, by default None + hierarchgrid : bool, optional + if True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + vertex deviation in rad or tuple of vertex deviation in rad and grid object + + Notes + ----- + - 2011-03-31 - Written - Bovy (NYU) """ # The following aren't actually the moments, but they are the moments # times the surface-mass density; that drops out @@ -814,50 +794,47 @@ def meanvR( integrate_method="dopr54_c", ): """ - NAME: - - meanvR - - PURPOSE: - - calculate the mean vR of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment(/ro) (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - surfacemass= if set use this pre-calculated surfacemass - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - mean vR - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the mean vR of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + radius at which to calculate the moment(/ro) (can be Quantity) + phi : float, optional + azimuth (rad unless deg=True; can be Quantity) (default=0.0) + t : float or array, optional + time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) (default=0.0) + surfacemass : float, optional + if set use this pre-calculated surfacemass (default=None) + nsigma : float, optional + number of sigma to integrate the velocities over (based on an estimate, so be generous) (default=None) + deg : bool, optional + azimuth is in degree (default=False); do not set this when giving phi as a Quantity + epsrel : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) (default=1.0e-02) + epsabs : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) (default=1.0e-05) + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid (default=None) + gridpoints : int, optional + number of points to use for the grid in 1D (default=101) + returnGrid : bool, optional + if True, return the grid object (default=False) + hierarchgrid : bool, optional + if True, use a hierarchical grid (default=False) + nlevels : int, optional + number of hierarchical levels for the hierarchical grid (default=2) + integrate_method : str, optional + orbit.integrate method argument (default="dopr54_c") + + Returns + ------- + float or tuple + mean vR + + Notes + ----- + -2011-03-31 - Written - Bovy (NYU) """ if isinstance(grid, evolveddiskdfGrid) or isinstance( grid, evolveddiskdfHierarchicalGrid @@ -966,50 +943,47 @@ def meanvT( integrate_method="dopr54_c", ): """ - NAME: - - meanvT - - PURPOSE: - - calculate the mean vT of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - surfacemass= if set use this pre-calculated surfacemass - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - mean vT - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the mean vT of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + radius at which to calculate the moment(/ro) (can be Quantity) + phi : float, optional + azimuth (rad unless deg=True; can be Quantity) (default=0.0) + t : float or array, optional + time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) (default=0.0) + surfacemass : float, optional + if set use this pre-calculated surfacemass (default=None) + nsigma : float, optional + number of sigma to integrate the velocities over (based on an estimate, so be generous) (default=None) + deg : bool, optional + azimuth is in degree (default=False); do not set this when giving phi as a Quantity + epsrel : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) (default=1.0e-02) + epsabs : float, optional + scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) (default=1.0e-05) + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid (default=None) + gridpoints : int, optional + number of points to use for the grid in 1D (default=101) + returnGrid : bool, optional + if True, return the grid object (default=False) + hierarchgrid : bool, optional + if True, use a hierarchical grid (default=False) + nlevels : int, optional + number of hierarchical levels for the hierarchical grid (default=2) + integrate_method : str, optional + orbit.integrate method argument (default="dopr54_c") + + Returns + ------- + float or tuple + mean vT + + Notes + ----- + -2011-03-31 - Written - Bovy (NYU) """ if isinstance(grid, evolveddiskdfGrid) or isinstance( grid, evolveddiskdfHierarchicalGrid @@ -1119,50 +1093,49 @@ def sigmaR2( integrate_method="dopr54_c", ): """ - NAME: - - sigmaR2 - - PURPOSE: - - calculate the radial variance of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - surfacemass, meanvR= if set use this pre-calculated surfacemass and mean vR - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - variance of vR - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the radial variance of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + t : float, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous) + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity) + epsrel : float, optional + Scipy.integrate keyword (the integration calculates the ratio of this vmoment to that of the initial DF) + epsabs : float, optional + Scipy.integrate keyword (the integration calculates the ratio of this vmoment to that of the initial DF) + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101) + returnGrid : bool, optional + If True, return the grid object (default=False) + surfacemass : float, optional + If set use this pre-calculated surfacemass + meanvR : float, optional + If set use this pre-calculated mean vR + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False) + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid + integrate_method : str, optional + orbit.integrate method argument + + Returns + ------- + float or tuple + Variance of vR. If returnGrid is True, return a tuple with the variance of vR and the grid object. + + Notes + ----- + - 2011-03-31 - Written - Bovy (NYU) """ # The following aren't actually the moments, but they are the moments # times the surface-mass density @@ -1297,50 +1270,49 @@ def sigmaT2( integrate_method="dopr54_c", ): """ - NAME: - - sigmaT2 - - PURPOSE: - - calculate the rotational-velocity variance of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - surfacemass, meanvT= if set use this pre-calculated surfacemass and mean rotational velocity - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - variance of vT - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the rotational-velocity variance of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrid : bool, optional + If True, return the grid object (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Variance of vT or tuple of variance of vT and grid object + + Notes + ----- + - 2011-03-31 - Written - Bovy (NYU) """ if isinstance(grid, evolveddiskdfGrid) or isinstance( grid, evolveddiskdfHierarchicalGrid @@ -1474,50 +1446,49 @@ def sigmaRT( integrate_method="dopr54_c", ): """ - NAME: - - sigmaRT - - PURPOSE: - - calculate the radial-rotational co-variance of the velocity distribution at (R,phi) - - INPUT: - - R - radius at which to calculate the moment (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - surfacemass, meanvR, meavT= if set use this pre-calculated surfacemass and mean vR and vT - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords (the integration calculates the ration of this vmoment to that of the initial DF) - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid object (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - covariance of vR and vT - - HISTORY: - - 2011-03-31 - Written - Bovy (NYU) - + Calculate the radial-rotational co-variance of the velocity distribution at (R,phi) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrid : bool, optional + If True, return the grid object (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Covariance of vR and vT or tuple of covariance of vR and vT and grid object + + Notes + ----- + - 2011-03-31 - Written - Bovy (NYU) """ # The following aren't actually the moments, but they are the moments # times the surface-mass density @@ -1675,54 +1646,57 @@ def oortA( integrate_method="dopr54_c", ): """ - NAME: - - oortA - - PURPOSE: - - calculate the Oort function A at (R,phi,t) - - INPUT: - - R - radius at which to calculate A (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - derivRGrid, derivphiGrid= if set to True, build a grid and use that to evaluate integrals of the derivatives of the DF;if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - derivGridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid objects (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - derivHierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - Oort A at R,phi,t - - HISTORY: - - 2011-10-16 - Written - Bovy (NYU) - + Calculate the Oort function A at (R,phi,t) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrids : bool, optional + If True, return the grid objects (default=False), by default False + derivRGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivphiGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivGridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + derivHierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Oort A at R,phi,t or tuple of Oort A at R,phi,t and grid objects + + Notes + ----- + - 2011-10-16 - Written - Bovy (NYU) """ # First calculate the grids if they are not given if isinstance(grid, bool) and grid: @@ -1975,54 +1949,57 @@ def oortB( integrate_method="dopr54_c", ): """ - NAME: - - oortB - - PURPOSE: - - calculate the Oort function B at (R,phi,t) - - INPUT: - - R - radius at which to calculate B (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - derivRGrid, derivphiGrid= if set to True, build a grid and use that to evaluat integrals of the derivatives of the DF: if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - derivGridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid objects (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - derivHierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - Oort B at R,phi,t - - HISTORY: - - 2011-10-16 - Written - Bovy (NYU) - + Calculate the Oort function B at (R,phi,t) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrids : bool, optional + If True, return the grid objects (default=False), by default False + derivRGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivphiGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivGridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + derivHierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Oort B at R,phi,t or tuple of Oort A at R,phi,t and grid objects + + Notes + ----- + - 2011-10-16 - Written - Bovy (NYU) """ # First calculate the grids if they are not given if isinstance(grid, bool) and grid: @@ -2275,54 +2252,57 @@ def oortC( integrate_method="dopr54_c", ): """ - NAME: - - oortC - - PURPOSE: - - calculate the Oort function C at (R,phi,t) - - INPUT: - - R - radius at which to calculate C (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - derivRGrid, derivphiGrid= if set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - derivGridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid objects (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - derivHierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - Oort C at R,phi,t - - HISTORY: - - 2011-10-16 - Written - Bovy (NYU) - + Calculate the Oort function C at (R,phi,t) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrids : bool, optional + If True, return the grid objects (default=False), by default False + derivRGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivphiGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivGridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + derivHierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Oort C at R,phi,t or tuple of Oort A at R,phi,t and grid objects + + Notes + ----- + - 2011-10-16 - Written - Bovy (NYU) """ # First calculate the grids if they are not given if isinstance(grid, bool) and grid: @@ -2575,54 +2555,57 @@ def oortK( integrate_method="dopr54_c", ): """ - NAME: - - oortK - - PURPOSE: - - calculate the Oort function K at (R,phi,t) - - INPUT: - - R - radius at which to calculate K (can be Quantity) - - phi= azimuth (rad unless deg=True; can be Quantity) - - t= time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity) - - nsigma - number of sigma to integrate the velocities over (based on an estimate, so be generous) - - deg= azimuth is in degree (default=False); do not set this when giving phi as a Quantity - - epsrel, epsabs - scipy.integrate keywords - - grid= if set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid - - derivRGrid, derivphiGrid= if set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid - - gridpoints= number of points to use for the grid in 1D (default=101) - - derivGridpoints= number of points to use for the grid in 1D (default=101) - - returnGrid= if True, return the grid objects (default=False) - - hierarchgrid= if True, use a hierarchical grid (default=False) - - derivHierarchgrid= if True, use a hierarchical grid (default=False) - - nlevels= number of hierarchical levels for the hierarchical grid - - integrate_method= orbit.integrate method argument - - OUTPUT: - - Oort K at R,phi,t - - HISTORY: - - 2011-10-16 - Written - Bovy (NYU) - + Calculate the Oort function K at (R,phi,t) + + Parameters + ---------- + R : float + Radius at which to calculate the moment (can be Quantity) + phi : float, optional + Azimuth (rad unless deg=True; can be Quantity), by default 0.0 + t : float or array, optional + Time at which to evaluate the DF (can be a list or ndarray; if this is the case, list needs to be in descending order and equally spaced) (can be Quantity), by default 0.0 + surfacemass : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + meanvT : float, optional + If set use this pre-calculated surfacemass and mean rotational velocity, by default None + nsigma : float, optional + Number of sigma to integrate the velocities over (based on an estimate, so be generous), by default None + deg : bool, optional + Azimuth is in degree (default=False); do not set this when giving phi as a Quantity, by default False + epsrel : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-02 + epsabs : float, optional + Scipy.integrate keywords (the integration calculates the ratio of this vmoment to that of the initial DF), by default 1.0e-05 + grid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + gridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + returnGrids : bool, optional + If True, return the grid objects (default=False), by default False + derivRGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivphiGrid : bool or evolveddiskdfGrid or evolveddiskdfHierarchicalGrid, optional + If set to True, build a grid and use that to evaluate integrals of the derivatives of the DF; if set to a grid-objects (such as returned by this procedure), use this grid, by default None + derivGridpoints : int, optional + Number of points to use for the grid in 1D (default=101), by default 101 + derivHierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + hierarchgrid : bool, optional + If True, use a hierarchical grid (default=False), by default False + nlevels : int, optional + Number of hierarchical levels for the hierarchical grid, by default 2 + integrate_method : str, optional + orbit.integrate method argument, by default "dopr54_c" + + Returns + ------- + float or tuple + Oort K at R,phi,t or tuple of Oort A at R,phi,t and grid objects + + Notes + ----- + - 2011-10-16 - Written - Bovy (NYU) """ # First calculate the grids if they are not given if isinstance(grid, bool) and grid: @@ -3140,16 +3123,30 @@ def __init__(self): def plot(self, tt=0): """ - NAME: - plot - PURPOSE: - plot the velocity distribution - INPUT: - t= optional time index - OUTPUT: - plot of velocity distribution to output device - HISTORY: - 2011-06-27 - Written - Bovy (NYU) + Plot the velocity distribution. + + Parameters + ---------- + tt : int, optional + Time index. Default is 0. + + Returns + ------- + None + + Other Parameters + ---------------- + xrange : list + Range of vRgrid. + yrange : list + Range of vTgrid. + plotthis : array + Array to be plotted. + + Notes + ----- + - 2011-06-27 - Written - Bovy (NYU) + """ xrange = [self.vRgrid[0], self.vRgrid[len(self.vRgrid) - 1]] yrange = [self.vTgrid[0], self.vTgrid[len(self.vTgrid) - 1]] @@ -3190,30 +3187,43 @@ def __init__( nlevelsTotal=None, ): """ - NAME: - __init__ - PURPOSE: - Initialize a hierarchical grid - INPUT: - edf - evolveddiskdf instance - R - Radius - phi- azimuth - nsigma - number of sigma to integrate over - t- time - sigmaR1 - radial dispersion - sigmaT1 - rotational-velocity dispersion - meanvR - mean of radial velocity - meanvT - mean of rotational velocity - gridpoints- number of gridpoints - nlevels- number of levels to build - deriv- None, 'R', or 'phi': calculates derivative of the moment wrt - R or phi - upperdxdy= area element of previous hierarchical level - print_progress= if True, print progress on building the grid - OUTPUT: - object - HISTORY: - 2011-04-21 - Written - Bovy (NYU) + Initialize a hierarchical grid + + Parameters + ---------- + edf : evolveddiskdf instance + R : float + Radius. + phi : float + Azimuth. + nsigma : int + Number of sigma to integrate over. + t : float + Time. + sigmaR1 : float + Radial dispersion. + sigmaT1 : float + Rotational-velocity dispersion. + meanvR : float + Mean of radial velocity. + meanvT : float + Mean of rotational velocity. + gridpoints : int + Number of gridpoints. + nlevels : int + Number of levels to build. + deriv : {None,'R','phi'} + Calculates derivative of the moment wrt R or phi. + upperdxdy : float, optional + Area element of previous hierarchical level. Default is None. + print_progress : bool, optional + If True, print progress on building the grid. Default is False. + nlevelsTotal : int, optional + Number of levels to build. Default is nlevels. + + Notes + ----- + - 2011-04-21 - Written - Bovy (NYU). """ self.sigmaR1 = sigmaR1 self.sigmaT1 = sigmaT1 @@ -3401,16 +3411,26 @@ def __call__(self, n, m): def plot(self, tt=0, vmax=None, aspect=None, extent=None): """ - NAME: - plot - PURPOSE: - plot the velocity distribution - INPUT: - t= optional time index - OUTPUT: - plot of velocity distribution to output device - HISTORY: - 2011-06-27 - Written - Bovy (NYU) + Plot the velocity distribution. + + Parameters + ---------- + tt : int, optional + Time index. Default is 0. + vmax : float, optional + Maximum value for the plot. Default is None. + aspect : float, optional + Aspect ratio of the plot. Default is None. + extent : list, optional + Extent of the plot. Default is None. + + Returns + ------- + None + + Notes + ----- + - 2011-06-27 - Written - Bovy (NYU). """ if vmax is None: vmax = self.max(tt=tt) * 2.0 diff --git a/galpy/df/isotropicHernquistdf.py b/galpy/df/isotropicHernquistdf.py index 44962f645..ed58d328e 100644 --- a/galpy/df/isotropicHernquistdf.py +++ b/galpy/df/isotropicHernquistdf.py @@ -12,29 +12,20 @@ class isotropicHernquistdf(isotropicsphericaldf): def __init__(self, pot=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an isotropic Hernquist distribution function - - INPUT: - - pot= (None) Hernquist Potential instance - - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-08-09 - Written - Lane (UofT) - + Initialize an isotropic Hernquist distribution function. + + Parameters + ---------- + pot : HernquistPotential instance + Hernquist potential instance. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-08-09 - Written - Lane (UofT) """ assert isinstance( pot, HernquistPotential @@ -53,25 +44,21 @@ def __init__(self, pot=None, ro=None, vo=None): 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: + Calculate the energy portion of an isotropic Hernquist distribution function - fE - The value of the energy portion of the DF + Parameters + ---------- + E : float or Quantity + The energy. - HISTORY: + Returns + ------- + numpy.ndarray + The value of the energy portion of the DF. - 2020-08-09 - Written - James Lane (UofT) + Notes + ----- + - 2020-08-09 - Written - James Lane (UofT) """ Etilde = -numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo) / self._psi0) # Handle E out of bounds diff --git a/galpy/df/isotropicNFWdf.py b/galpy/df/isotropicNFWdf.py index c2daffae7..e7dc23508 100644 --- a/galpy/df/isotropicNFWdf.py +++ b/galpy/df/isotropicNFWdf.py @@ -28,31 +28,24 @@ class isotropicNFWdf(isotropicsphericaldf): def __init__(self, pot=None, widrow=False, rmax=1e4, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an isotropic NFW distribution function - - INPUT: - - pot= (None) NFW Potential instance - - widrow= (False) if True, use the approximate form from Widrow (2000), otherwise use improved fit that has <~1e-5 relative density errors - - rmax= (1e4) maximum radius to consider (can be Quantity); set to numpy.inf to evaluate NFW w/o cut-off - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2021-02-01 - Written - Bovy (UofT) + Initialize an isotropic NFW distribution function + + Parameters + ---------- + pot : NFWPotential instance + NFW Potential instance + widrow : bool, optional + If True, use the approximate form from Widrow (2000), otherwise use improved fit that has <~1e-5 relative density errors + rmax : float or Quantity, optional + Maximum radius to consider; set to numpy.inf to evaluate NFW w/o cut-off + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2021-02-01 - Written - Bovy (UofT) """ assert isinstance(pot, NFWPotential), "pot= must be potential.NFWPotential" @@ -66,25 +59,21 @@ def __init__(self, pot=None, widrow=False, rmax=1e4, ro=None, vo=None): def fE(self, E): """ - NAME: - - fE - - PURPOSE - - Calculate the energy portion of an isotropic NFW distribution function - - INPUT: - - E - The energy (can be Quantity) - - OUTPUT: + Calculate the energy portion of an isotropic NFW distribution function - fE - The value of the energy portion of the DF + Parameters + ---------- + E : float or Quantity + The energy. - HISTORY: + Returns + ------- + ndarray + The value of the energy portion of the DF. - 2021-02-01 - Written - Bovy (UofT) + Notes + ----- + - 2021-02-01 - Written - Bovy (UofT) """ Etilde = -conversion.parse_energy(E, vo=self._vo) / self._Etildemax out = numpy.zeros_like(Etilde) diff --git a/galpy/df/isotropicPlummerdf.py b/galpy/df/isotropicPlummerdf.py index d40e62c48..1eea5cb95 100644 --- a/galpy/df/isotropicPlummerdf.py +++ b/galpy/df/isotropicPlummerdf.py @@ -18,28 +18,20 @@ class isotropicPlummerdf(isotropicsphericaldf): def __init__(self, pot=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an isotropic Plummer distribution function - - INPUT: - - pot= (None) Plummer Potential instance - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-10-01 - Written - Bovy (UofT) - + Initialize an isotropic Plummer distribution function + + Parameters + ---------- + pot : Potential object + Plummer Potential instance + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-10-01 - Written - Bovy (UofT). """ assert isinstance( pot, PlummerPotential @@ -58,25 +50,22 @@ def __init__(self, pot=None, ro=None, vo=None): def fE(self, E): """ - NAME: - - fE - - PURPOSE - - Calculate the energy portion of an isotropic Plummer distribution function - - INPUT: - - E - The energy (can be Quantity) + Calculate the energy portion of an isotropic Plummer distribution function. - OUTPUT: + Parameters + ---------- + E : float or Quantity + The energy. - fE - The value of the energy portion of the DF + Returns + ------- + ndarray + The value of the energy portion of the DF. - HISTORY: + Notes + ----- + - 2020-10-01 - Written - Bovy (UofT) - 2020-10-01 - Written - Bovy (UofT) """ Etilde = -conversion.parse_energy(E, vo=self._vo) out = numpy.zeros_like(Etilde) diff --git a/galpy/df/jeans.py b/galpy/df/jeans.py index 163da5866..ec7154b02 100644 --- a/galpy/df/jeans.py +++ b/galpy/df/jeans.py @@ -17,32 +17,27 @@ @physical_conversion("velocity", pop=True) def sigmar(Pot, r, dens=None, beta=0.0): """ - NAME: - - sigmar - - PURPOSE: - - Compute the radial velocity dispersion using the spherical Jeans equation - - INPUT: - - Pot - potential or list of potentials (evaluated at R=r/sqrt(2),z=r/sqrt(2), sphericity not checked) - - r - Galactocentric radius (can be Quantity) - - dens= (None) tracer density profile (function of r); if None, the density is assumed to be that corresponding to the potential - - beta= (0.) anisotropy; can be a constant or a function of r - - OUTPUT: - - sigma_r(r) - - HISTORY: - - 2018-07-05 - Written - Bovy (UofT) - + Compute the radial velocity dispersion using the spherical Jeans equation + + Parameters + ---------- + Pot : potential or list of potentials + Gravitational potential; evaluated at R=r/sqrt(2),z=r/sqrt(2), sphericity not checked. + r : float or Quantity + Galactocentric radius + dens : function, optional + tracer density profile (function of r); if None, the density is assumed to be that corresponding to the potential + beta : float or function, optional + anisotropy; can be a constant or a function of r + + Returns + ------- + float + sigma_r(r) + + Notes + ----- + - 2018-07-05 - Written - Bovy (UofT) """ Pot = flatten_pot(Pot) if dens is None: @@ -82,36 +77,31 @@ def sigmar(Pot, r, dens=None, beta=0.0): @physical_conversion("velocity", pop=True) def sigmalos(Pot, R, dens=None, surfdens=None, beta=0.0, sigma_r=None): """ - NAME: - - sigmalos - - PURPOSE: - - Compute the line-of-sight velocity dispersion using the spherical Jeans equation - - INPUT: - - Pot - potential or list of potentials (evaluated at R=r/sqrt(2),z=r/sqrt(2), sphericity not checked) - - R - Galactocentric projected radius (can be Quantity) - - dens= (None) tracer density profile (function of r); if None, the density is assumed to be that corresponding to the potential - - surfdens= (None) tracer surface density profile (value at R or function of R); if None, the surface density is assumed to be that corresponding to the density - - beta= (0.) anisotropy; can be a constant or a function of r - - sigma_r= (None) if given, the solution of the spherical Jeans equation sigma_r(r) (used instead of solving the Jeans equation as part of this routine) - - OUTPUT: - - sigma_los(R) - - HISTORY: - - 2018-08-27 - Written - Bovy (UofT) - + Compute the line-of-sight velocity dispersion using the spherical Jeans equation + + Parameters + ---------- + Pot : potential or list of potentials + Gravitational potential; evaluated at R=r/sqrt(2),z=r/sqrt(2), sphericity not checked. + R : float or Quantity + Galactocentric projected radius + dens : function, optional + tracer density profile (function of r); if None, the density is assumed to be that corresponding to the potential + surfdens : float or function, optional + tracer surface density profile (value at R or function of R); if None, the surface density is assumed to be that corresponding to the density + beta : float or function, optional + anisotropy; can be a constant or a function of r + sigma_r : float or function, optional + if given, the solution of the spherical Jeans equation sigma_r(r) (used instead of solving the Jeans equation as part of this routine) + + Returns + ------- + float + sigma_los(R) + + Notes + ----- + - 2018-08-27 - Written - Bovy (UofT) """ Pot = flatten_pot(Pot) if dens is None: diff --git a/galpy/df/kingdf.py b/galpy/df/kingdf.py index fc13de160..72b503ebc 100644 --- a/galpy/df/kingdf.py +++ b/galpy/df/kingdf.py @@ -23,34 +23,26 @@ class kingdf(isotropicsphericaldf): def __init__(self, W0, M=1.0, rt=1.0, 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) - + Initialize a King DF + + Parameters + ---------- + W0 : float + Dimensionless central potential :math:`W_0 = \\Psi(0)/\\sigma^2` (in practice, needs to be :math:`\\lesssim 200`, where the DF is essentially isothermal). + M : float or Quantity, optional + Total mass. + rt : float or Quantity, optional + Tidal radius. + npt : int, optional + Number of points to use to solve for :math:`\\Psi(r)`. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-07-09 - Written - Bovy (UofT) """ # Just run df init to set up unit-conversion parameters df.__init__(self, ro=ro, vo=vo) diff --git a/galpy/df/osipkovmerrittHernquistdf.py b/galpy/df/osipkovmerrittHernquistdf.py index 6fd763334..8a8c34b7e 100644 --- a/galpy/df/osipkovmerrittHernquistdf.py +++ b/galpy/df/osipkovmerrittHernquistdf.py @@ -20,29 +20,22 @@ class osipkovmerrittHernquistdf(_osipkovmerrittdf): def __init__(self, pot=None, ra=1.4, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a Hernquist DF with Osipkov-Merritt anisotropy - - INPUT: - - pot - Hernquist potential which determines the DF - - ra - anisotropy radius (can be a Quantity) - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-11-12 - Written - Bovy (UofT) + Initialize a Hernquist DF with Osipkov-Merritt anisotropy + + Parameters + ---------- + pot : potential.HernquistPotential + Hernquist potential which determines the DF + ra : float or Quantity, optional + Anisotropy radius + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-11-12 - Written - Bovy (UofT) """ assert isinstance( pot, HernquistPotential @@ -62,26 +55,21 @@ def __init__(self, pot=None, ra=1.4, ro=None, vo=None): def fQ(self, Q): """ - NAME: - - fQ - - PURPOSE - - Calculate the f(Q) portion of an Osipkov-Merritt Hernquist distribution function - - INPUT: - - Q - The Osipkov-Merritt 'energy' E-L^2/[2ra^2] (can be Quantity) - - OUTPUT: - - fQ - The value of the f(Q) portion of the DF + Calculate the f(Q) portion of an Osipkov-Merritt Hernquist distribution function - HISTORY: + Parameters + ---------- + Q : float or numpy.ndarray + The Osipkov-Merritt 'energy' E-L^2/[2ra^2] (can be Quantity) - 2020-11-12 - Written - Bovy (UofT) + Returns + ------- + float or numpy.ndarray + The value of the f(Q) portion of the DF + Notes + ----- + - 2020-11-12 - Written - Bovy (UofT) """ Qtilde = numpy.atleast_1d(conversion.parse_energy(Q, vo=self._vo) / self._psi0) # Handle potential Q outside of bounds diff --git a/galpy/df/osipkovmerrittNFWdf.py b/galpy/df/osipkovmerrittNFWdf.py index 592a82ded..a149fd428 100644 --- a/galpy/df/osipkovmerrittNFWdf.py +++ b/galpy/df/osipkovmerrittNFWdf.py @@ -35,31 +35,24 @@ class osipkovmerrittNFWdf(_osipkovmerrittdf): def __init__(self, pot=None, ra=1.4, rmax=1e4, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a NFW DF with Osipkov-Merritt anisotropy - - INPUT: - - pot - NFW potential which determines the DF - - ra - anisotropy radius (can be a Quantity) - - rmax= (1e4) maximum radius to consider (can be Quantity); set to numpy.inf to evaluate NFW w/o cut-off - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-11-12 - Written - Bovy (UofT) + Initialize a NFW DF with Osipkov-Merritt anisotropy + + Parameters + ---------- + pot : potential.NFWPotential + NFW potential which determines the DF + ra : float or Quantity, optional + Anisotropy radius + rmax : float or Quantity, optional + Maximum radius to consider (set to numpy.inf to evaluate NFW w/o cut-off) + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-11-12 - Written - Bovy (UofT) """ assert isinstance(pot, NFWPotential), "pot= must be potential.NFWPotential" _osipkovmerrittdf.__init__(self, pot=pot, ra=ra, rmax=rmax, ro=ro, vo=vo) @@ -74,25 +67,21 @@ def __init__(self, pot=None, ra=1.4, rmax=1e4, ro=None, vo=None): def fQ(self, Q): """ - NAME: - - fQ - - PURPOSE - - Calculate the f(Q) portion of an Osipkov-Merritt NFW distribution function - - INPUT: - - Q - The Osipkov-Merritt 'energy' E-L^2/[2ra^2] (can be Quantity) - - OUTPUT: + Calculate the f(Q) portion of an Osipkov-Merritt NFW distribution function - fQ - The value of the f(Q) portion of the DF + Parameters + ---------- + Q : float or Quantity + The Osipkov-Merritt 'energy' E-L^2/[2ra^2] - HISTORY: + Returns + ------- + ndarray + The value of the f(Q) portion of the DF - 2021-02-09 - Written - Bovy (UofT) + Notes + ----- + - 2021-02-09 - Written - Bovy (UofT) """ Qtilde = conversion.parse_energy(Q, vo=self._vo) / self._Qtildemax diff --git a/galpy/df/osipkovmerrittdf.py b/galpy/df/osipkovmerrittdf.py index 3065151a2..b8e4905bc 100644 --- a/galpy/df/osipkovmerrittdf.py +++ b/galpy/df/osipkovmerrittdf.py @@ -18,37 +18,28 @@ def __init__( self, pot=None, denspot=None, ra=1.4, rmax=None, scale=None, ro=None, vo=None ): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a DF with Osipkov-Merritt anisotropy - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - ra - anisotropy radius (can be a Quantity) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale - Characteristic scale radius to aid sampling calculations. - Not necessary, and will also be overridden by value from pot - if available. - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-11-12 - Written - Bovy (UofT) + Initialize a DF with Osipkov-Merritt anisotropy. + + Parameters + ---------- + pot : Potential instance or list thereof, optional + Default: None + denspot : Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot), optional + Default: None + ra : float or Quantity, optional + Anisotropy radius. Default: 1.4 + rmax : float or Quantity, optional + Maximum radius to consider; DF is cut off at E = Phi(rmax). Default: None + scale : float or Quantity, optional + Characteristic scale radius to aid sampling calculations. Not necessary, and will also be overridden by value from pot if available. Default: None + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-11-12 - Written - Bovy (UofT) """ anisotropicsphericaldf.__init__( @@ -59,27 +50,23 @@ def __init__( def _call_internal(self, *args): """ - NAME: - - _call_internal - - PURPOSE: - - Evaluate the DF for an Osipkov-Merritt-anisotropy DF - - INPUT: - - E - The energy - - L - The angular momentum - - OUTPUT: + Evaluate the DF for an Osipkov-Merritt-anisotropy DF - fH - The value of the DF + Parameters + ---------- + E : float + The energy + L : float + The angular momentum - HISTORY: + Returns + ------- + float + The value of the DF - 2020-11-12 - Written - Bovy (UofT) + Notes + ----- + - 2020-11-12 - Written - Bovy (UofT) """ E, L, _ = args @@ -217,36 +204,28 @@ def __init__( self, pot=None, denspot=None, ra=1.4, rmax=1e4, scale=None, ro=None, vo=None ): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a DF with Osipkov-Merritt anisotropy - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - ra - anisotropy radius (can be a Quantity) - - rmax= (1e4) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale - Characteristic scale radius to aid sampling calculations. Optionaland will also be overridden by value from pot if available. - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2021-02-07 - Written - Bovy (UofT) - + Initialize a DF with Osipkov-Merritt anisotropy. + + Parameters + ---------- + pot : Potential instance or list thereof, optional + Default: None + denspot : Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot), optional + Default: None + ra : float or Quantity, optional + Anisotropy radius. Default: 1.4 + rmax : float or Quantity, optional + Maximum radius to consider; DF is cut off at E = Phi(rmax). Default: None + scale : float or Quantity, optional + Characteristic scale radius to aid sampling calculations. Not necessary, and will also be overridden by value from pot if available. Default: None + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2021-02-07 - Written - Bovy (UofT) """ _osipkovmerrittdf.__init__( self, pot=pot, denspot=denspot, ra=ra, rmax=rmax, scale=scale, ro=ro, vo=vo @@ -316,25 +295,21 @@ def sample(self, R=None, z=None, phi=None, n=1, return_orbit=True, rmin=0.0): def fQ(self, Q): """ - NAME: - - fQ - - PURPOSE - - Calculate the f(Q) portion of an Osipkov-Merritt Hernquist distribution function - - INPUT: - - Q - The Osipkov-Merritt 'energy' E-L^2/[2ra^2] (can be Quantity) - - OUTPUT: + Calculate the f(Q) portion of an Osipkov-Merritt Hernquist distribution function - fQ - The value of the f(Q) portion of the DF + Parameters + ---------- + Q : float + The Osipkov-Merritt 'energy' E-L^2/[2ra^2] - HISTORY: + Returns + ------- + float + The value of the f(Q) portion of the DF - 2021-02-07 - Written - Bovy (UofT) + Notes + ----- + - 2021-02-07 - Written - Bovy (UofT) """ return self._edf.fE(-Q) diff --git a/galpy/df/quasiisothermaldf.py b/galpy/df/quasiisothermaldf.py index 3d3bde05c..7b8451917 100644 --- a/galpy/df/quasiisothermaldf.py +++ b/galpy/df/quasiisothermaldf.py @@ -54,56 +54,44 @@ def __init__( 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) - + Initialize a quasi-isothermal DF + + Parameters + ---------- + hr : float or Quantity + Radial scale length. + sr : float or Quantity + Radial velocity dispersion at the solar radius. + sz : float or Quantity + Vertical velocity dispersion at the solar radius. + hsr : float or Quantity + Radial-velocity-dispersion scale length. + hsz : float or Quantity + Vertial-velocity-dispersion scale length. + pot : Potential or list thereof + Potential or list of potentials that represents the underlying potential. + aA : actionAngle instance + 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 : bool, optional + If True, set counter-rotating stars' DF to zero. + refr : float or Quantity, optional + Reference radius for dispersions (can be different from ro). + lo : float or Quantity, optional + Reference angular momentum below where there are significant numbers of retrograde stars. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + _precomputerg : bool, optional + If True (default), pre-compute the rL(L). + _precomputergrmax : float or Quantity, optional + If set, this is the maximum R for which to pre-compute rg (default: 5*hr). + _precomputergnLz : int, optional + If set, number of Lz to pre-compute rg for (default: 51). + + Notes + ----- + - 2012-07-25 - Started - Bovy (IAS@MPIA) """ df.__init__(self, ro=ro, vo=vo) self._hr = parse_length(hr, ro=self._ro) @@ -176,56 +164,40 @@ def __init__( @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 - + Evaluate the DF + + Parameters + ---------- + args: tuple or Orbit + 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: bool, optional + If True, return the natural log. + func: function of (jr,lz,jz), optional + Function of the actions to multiply the DF with (useful for moments). + _return_actions: bool, optional + If True, return the actions as well. + _return_freqs: bool, optional + If True, return the frequencies as well. + _return_rgr: bool, optional + If True, return the rg as well. + kwargs: dict, optional + scipy.integrate.quadrature kwargs. + + Returns + ------- + float + Value of DF. + + Notes + ----- + - 2012-07-25 - Written - Bovy (IAS@MPIA) """ # First parse log log = kwargs.pop("log", False) @@ -346,34 +318,28 @@ def __call__(self, *args, **kwargs): @physical_conversion("position", pop=True) def estimate_hr(self, R, z=0.0, dR=10.0**-8.0, **kwargs): """ - NAME: + Estimate the exponential scale length at R. - 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) + Parameters + ---------- + R : float or Quantity + Galactocentric radius. + z : float or Quantity, optional + Height (default: 0 pc). + dR : float or Quantity, optional + Range in R to use. + **kwargs + Density kwargs. - 2013-01-28 - Re-written - Bovy + Returns + ------- + float or Quantity + Estimated hR. + Notes + ----- + - 2012-09-11 - Written - Bovy (IAS) + - 2013-01-28 - Re-written - Bovy """ Rs = [R - dR / 2.0, R + dR / 2.0] if z is None: @@ -391,32 +357,28 @@ def estimate_hr(self, R, z=0.0, dR=10.0**-8.0, **kwargs): @physical_conversion("position", pop=True) def estimate_hz(self, R, z, dz=10.0**-8.0, **kwargs): """ - NAME: - - estimate_hz - - PURPOSE: + Estimate the exponential scale height at R. - 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) + Parameters + ---------- + R : float or Quantity + Galactocentric radius. + z : float or Quantity + Height above the Galactic plane. + dz : float or Quantity, optional + z range to use. + **kwargs + density kwargs. - 2013-01-28 - Re-written - Bovy + Returns + ------- + float or Quantity + Estimated hz. + Notes + ----- + - 2012-08-30 - Written - Bovy (IAS) + - 2013-01-28 - Re-written - Bovy """ if z == 0.0: zs = [z, z + dz] @@ -432,31 +394,27 @@ def estimate_hz(self, R, z, dz=10.0**-8.0, **kwargs): @physical_conversion("position", pop=True) def estimate_hsr(self, R, z=0.0, dR=10.0**-8.0, **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) + Estimate the exponential scale length of the radial dispersion at R. - dR- range in R to use (can be Quantity) - - density kwargs - - OUTPUT: - - estimated hsR + Parameters + ---------- + R : float or Quantity + Galactocentric radius. + z : float or Quantity, optional + Height (default: 0 pc). + dR : float or Quantity, optional + Range in R to use. + **kwargs + Density kwargs. - HISTORY: + Returns + ------- + float or Quantity + Estimated hsR. - 2013-03-08 - Written - Bovy (IAS) + Notes + ----- + - 2013-03-08 - Written - Bovy (IAS) """ Rs = [R - dR / 2.0, R + dR / 2.0] @@ -468,31 +426,27 @@ def estimate_hsr(self, R, z=0.0, dR=10.0**-8.0, **kwargs): @physical_conversion("position", pop=True) def estimate_hsz(self, R, z=0.0, dR=10.0**-8.0, **kwargs): """ - NAME: - - estimate_hsz - - PURPOSE: - - estimate the exponential scale length of the vertical dispersion at R - - INPUT: + Estimate the exponential scale length of the vertical dispersion at R. - 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 + Parameters + ---------- + R : float or Quantity + Galactocentric radius. + z : float or Quantity, optional + Height (default: 0 pc). + dR : float or Quantity, optional + Range in R to use. + **kwargs + Density kwargs. - HISTORY: + Returns + ------- + float or Quantity + Estimated hsz. - 2013-03-08 - Written - Bovy (IAS) + Notes + ----- + - 2013-03-08 - Written - Bovy (IAS) """ Rs = [R - dR / 2.0, R + dR / 2.0] @@ -505,37 +459,32 @@ def estimate_hsz(self, R, z=0.0, dR=10.0**-8.0, **kwargs): def surfacemass_z( self, R, nz=7, zmax=1.0, fixed_quad=True, fixed_order=8, **kwargs ): - r""" - 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) - + """ + Calculate the vertically-integrated surface density. + + Parameters + ---------- + R : float or Quantity + Galactocentric radius. + nz : int, optional + Number of zs to use to estimate. Default is 7. + zmax : float or Quantity, optional + Maximum z to use. Default is 1.0. + fixed_quad : bool, optional + If True (default), use Gauss-Legendre integration. + fixed_order : int, optional + Order of GL integration to use. Default is 8. + **kwargs : dict + Density kwargs. + + Returns + ------- + float or Quantity + Surface density at R. + + Notes + ----- + - 2012-08-30 - Written - Bovy (IAS) """ if fixed_quad: return ( @@ -559,50 +508,45 @@ def surfacemass_z( 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: - - at R,z (no support for units) - - HISTORY: - - 2012-08-06 - Written - Bovy (IAS@MPIA) + Calculate the an arbitrary moment of the velocity distribution at R times the density + + Parameters + ---------- + R : float + radius at which to calculate the moment(/ro) + z : float + height at which to calculate the moment(/ro) + n : int + vR^n + m : int + vT^m + o : int + vz^o + nsigma : int, optional + number of sigma to integrate the vR and vz velocities over (when doing explicit numerical integral; default: 4) + vTmax : float, optional + upper limit for integration over vT (default: 1.5) + mc : bool, optional + if True, calculate using Monte Carlo integration + nmc : int, optional + if mc, use nmc samples + gl : bool, optional + use Gauss-Legendre + _returngl : bool, optional + if True, return the evaluated DF + _return_actions : bool, optional + if True, return the evaluated actions (does not work with _returngl currently) + _return_freqs : bool, optional + if True, return the evaluated frequencies and rg (does not work with _returngl currently) + + Returns + ------- + float + at R,z (no support for units) + + Notes + ----- + - 2012-08-06 - Written - Bovy (IAS@MPIA) """ use_physical = kwargs.pop("use_physical", True) @@ -969,39 +913,35 @@ def _vmomentdensity( 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: - - at R (no support for units) - - HISTORY: - - 2012-08-09 - Written - Bovy (IAS@MPIA) + Calculate the an arbitrary moment of an action of the velocity distribution at R times the surfacmass. + + Parameters + ---------- + R : float + radius at which to calculate the moment(/ro) + z : float + height at which to calculate the moment(/ro) + n : int + jr^n + m : int + lz^m + o : int + jz^o + nsigma : int, optional + Number of sigma to integrate the velocities over (when doing explicit numerical integral). Default is None. + mc : bool, optional + If True, calculate using Monte Carlo integration. Default is False. + nmc : int, optional + If mc is True, use nmc samples. Default is None. + + Returns + ------- + float or Quantity + at R (no support for units) + + Notes + ----- + - 2012-08-09 - Written - Bovy (IAS@MPIA) """ use_physical = kwargs.pop("use_physical", True) @@ -1131,41 +1071,35 @@ 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) + Calculate the density at R,z by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate the density. + z : float or Quantity + Height at which to calculate the density. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + Density at (R,z). + + Notes + ----- + - 2012-07-26 - Written - Bovy (IAS@MPIA) """ return self._vmomentdensity( @@ -1178,41 +1112,35 @@ 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) + Calculate sigma_R^2 by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : int, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + sigma_R^2. + + Notes + ----- + - 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: @@ -1269,41 +1197,35 @@ 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 + Calculate sigma_RZ^2 by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : int, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs + scipy.integrate.tplquad kwargs epsabs and epsrel. - HISTORY: + Returns + ------- + float + sigma_Rz^2. - 2012-07-30 - Written - Bovy (IAS@MPIA) + Notes + ----- + - 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: @@ -1360,49 +1282,35 @@ 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) - + Calculate the tilt of the velocity ellipsoid by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : int, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + + Returns + ------- + float + Tilt in radians. + + Notes + ----- + - 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, @@ -1503,41 +1411,35 @@ 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) + Calculate sigma_z^2 by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : int, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + sigma_z^2. + + Notes + ----- + - 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: @@ -1594,41 +1496,35 @@ 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) + Calculate the mean rotational velocity by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + Mean rotational velocity. + + Notes + ----- + - 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: @@ -1685,41 +1581,35 @@ 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) + Calculate the mean radial velocity by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + Mean radial velocity. + + Notes + ----- + - 2012-12-23 - Written - Bovy (IAS) """ if mc: @@ -1776,42 +1666,35 @@ 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) - + Calculate the mean vertical velocity by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs : dict, optional + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + Mean vertical velocity + + Notes + ----- + - 2012-12-23 - Written - Bovy (IAS) """ if mc: surfmass, vrs, vts, vzs = self._vmomentdensity( @@ -1867,41 +1750,35 @@ 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 + Calculate sigma_T^2 by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : int, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc is True, use nmc samples. + gl : bool, optional + If True, calculate using Gauss-Legendre integration. + ngl : int, optional + If gl is True, use ngl-th order Gauss-Legendre integration for each dimension. + **kwargs + scipy.integrate.tplquad kwargs epsabs and epsrel. - HISTORY: + Returns + ------- + float + sigma_T^2. - 2012-07-30 - Written - Bovy (IAS@MPIA) + Notes + ----- + - 2012-07-30 - Written - Bovy (IAS@MPIA) """ if mc: @@ -1994,37 +1871,31 @@ def sigmaT2( @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) + Calculate the mean radial action by marginalizing over velocity + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this + z : float or Quantity + Height at which to calculate this + nsigma : float, optional + Number of sigma to integrate the velocities over + mc : bool, optional + If True, calculate using Monte Carlo integration + nmc : int, optional + If mc, use nmc samples + **kwargs : dict + scipy.integrate.tplquad kwargs epsabs and epsrel + + Returns + ------- + float + Mean jr + + Notes + ----- + - 2012-08-09 - Written - Bovy (IAS@MPIA) """ if mc: @@ -2069,39 +1940,34 @@ def meanjr(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): @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 + Calculate the mean angular momentum by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + **kwargs + scipy.integrate.tplquad kwargs epsabs and epsrel. - HISTORY: + Returns + ------- + float + Mean angular momentum. - 2012-08-09 - Written - Bovy (IAS@MPIA) + Notes + ----- + - 2012-08-09 - Written - Bovy (IAS@MPIA) """ + if mc: surfmass, vrs, vts, vzs = self._vmomentdensity( R, @@ -2144,37 +2010,31 @@ def meanlz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): @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) + Calculate the mean vertical action by marginalizing over velocity. + + Parameters + ---------- + R : float or Quantity + Radius at which to calculate this. + z : float or Quantity + Height at which to calculate this. + nsigma : float, optional + Number of sigma to integrate the velocities over. + mc : bool, optional + If True, calculate using Monte Carlo integration. + nmc : int, optional + If mc, use nmc samples. + **kwargs : dict + scipy.integrate.tplquad kwargs epsabs and epsrel. + + Returns + ------- + float + Mean jz. + + Notes + ----- + - 2012-08-09 - Written - Bovy (IAS@MPIA) """ if mc: @@ -2218,30 +2078,25 @@ def meanjz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): @potential_physical_input def sampleV(self, R, z, n=1, **kwargs): """ - NAME: - - sampleV + Sample a radial, azimuthal, and vertical velocity at R,z - PURPOSE: + Parameters + ---------- + R : float or Quantity + Galactocentric distance. + z : float or Quantity + Height. + n : int, optional + Number of distances to sample. - 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) + Returns + ------- + list + List of samples. + Notes + ----- + - 2012-12-17 - Written - Bovy (IAS@MPIA) """ use_physical = kwargs.pop("use_physical", True) vo = kwargs.pop("vo", None) @@ -2317,38 +2172,35 @@ def sampleV_interpolate( **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) - + Sample radial, azimuthal, and vertical velocity at R,z using interpolation. + + Parameters + ---------- + R : numpy.ndarray or Quantity + Galactocentric distance. + z : numpy.ndarray or Quantity + Height. + R_pixel : float + The pixel size for creating the grid for interpolation (in natural units). + z_pixel : float + The pixel size for creating the grid for interpolation (in natural units). + num_std : float, optional + Number of standard deviation to be considered outliers sampled separately from interpolation. + R_min : float, optional + Minimum R value for the grid. + R_max : float, optional + Maximum R value for the grid. + z_max : float, optional + Maximum z value for the grid. + + Returns + ------- + numpy.ndarray + A numpy array containing the sampled velocity, (vR, vT, vz), where each row corresponds to the row of (R,z). + + Notes + ----- + - 2018-08-10 - Written - Samuel Wong (University of Toronto) """ use_physical = kwargs.pop("use_physical", True) vo = kwargs.pop("vo", None) @@ -2441,31 +2293,25 @@ def sampleV_interpolate( 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 + Sample a radial, azimuthal, and vertical velocity at R,z. - OUTPUT: + Parameters + ---------- + R : float or numpy.ndarray + Galactocentric distance. + z : float or numpy.ndarray + Height. + maxVT : numpy.ndarray + An array of pre-optimized maximum vT at corresponding R,z. - a numpy array containing the sampled velocity, (vR, vT, vz), - where each row correspond to the row of (R,z) + Returns + ------- + numpy.ndarray + 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) + Notes + ----- + - 2018-08-10 - Written - Samuel Wong (University of Toronto) """ length = numpy.size(R) @@ -2518,37 +2364,33 @@ def _sampleV_preoptimized(self, R, z, maxVT): @physical_conversion("phasespacedensityvelocity2", pop=True) def pvR(self, vR, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0, 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) + Calculate the marginalized vR probability at this location (NOT normalized by the density). + + Parameters + ---------- + vR : float or Quantity + Radial velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + nsigma : float, optional + Number of sigma to integrate the velocities over. + vTmax : float, optional + Sets integration limits to [0,vTmax] for integration over vT. + + Returns + ------- + float + p(vR,R,z). + + Notes + ----- + - 2012-12-22 - Written - Bovy (IAS@MPIA) """ sigmaz1 = self._sz * numpy.exp((self._refr - R) / self._hsz) @@ -2614,36 +2456,32 @@ def pvR(self, vR, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0, vTmax=1.5): @physical_conversion("phasespacedensityvelocity2", pop=True) def pvT(self, vT, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0): """ - NAME: - - pvT - - PURPOSE: - - calculate the marginalized vT probability at this location (NOT normalized by the density) - - INPUT: - - vT - rotational 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) + Calculate the marginalized vT probability at this location (NOT normalized by the density). + + Parameters + ---------- + vT : float or Quantity + Azimuthal velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + nsigma : float, optional + Number of sigma to integrate the velocities over. + + Returns + ------- + float + p(vT,R,z). + + Notes + ----- + - 2012-12-22 - Written - Bovy (IAS@MPIA) + - 2018-01-12 - Added Gauss-Legendre integration prefactor nsigma^2/4 - Trick (MPA) """ sigmaR1 = self._sr * numpy.exp((self._refr - R) / self._hsr) @@ -2741,38 +2579,33 @@ def pvz( _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) - + Calculate the marginalized vz probability at this location (NOT normalized by the density). + + Parameters + ---------- + vz : float or Quantity + Vertical velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + nsigma : float, optional + Number of sigma to integrate the velocities over. + vTmax : float, optional + Sets integration limits to [0,vTmax] for integration over vT. + + Returns + ------- + float + p(vz,R,z). + + Notes + ----- + - 2012-12-22 - Written - Bovy (IAS) """ if _sigmaR1 is None: sigmaR1 = self._sr * numpy.exp((self._refr - R) / self._hsr) @@ -2920,39 +2753,34 @@ def pvz( @physical_conversion("phasespacedensityvelocity", pop=True) def pvRvT(self, vR, vT, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0): """ - 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 - rotational 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) - + Calculate the marginalized (vR,vT) probability at this location (NOT normalized by the density). + + Parameters + ---------- + vR : float or Quantity + Radial velocity. + vT : float or Quantity + Azimuthal velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + nsigma : float, optional + Number of sigma to integrate the velocities over. + + Returns + ------- + float + p(vR,vT,R,z). + + Notes + ----- + - 2012-12-22 - 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: @@ -3007,38 +2835,34 @@ def pvRvT(self, vR, vT, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0): @physical_conversion("phasespacedensityvelocity", pop=True) def pvTvz(self, vT, vz, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0): """ - NAME: - - pvTvz - - PURPOSE: - - calculate the marginalized (vT,vz) probability at this location (NOT normalized by the density) - - INPUT: - - vT - rotational 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) + Calculate the marginalized (vT,vz) probability at this location (NOT normalized by the density). + + Parameters + ---------- + vT : float or Quantity + Azimuthal velocity. + vz : float or Quantity + Vertical velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + nsigma : float, optional + Number of sigma to integrate the velocities over. + + Returns + ------- + float or Quantity + p(vT,vz,R,z). + + Notes + ----- + - 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) @@ -3094,39 +2918,34 @@ def pvTvz(self, vT, vz, R, z, gl=True, ngl=_DEFAULTNGL2, nsigma=4.0): @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) - + Calculate the marginalized (vR,vz) probability at this location (NOT normalized by the density). + + Parameters + ---------- + vR : float or Quantity + Radial velocity. + vz : float or Quantity + Vertical velocity. + R : float or Quantity + Radius. + z : float or Quantity + Height. + gl : bool, optional + If True, use Gauss-Legendre integration. + ngl : int, optional + If gl, use ngl-th order Gauss-Legendre integration for each dimension. + vTmax : float, optional + Sets integration limits to [0,vTmax] for integration over vT. + + Returns + ------- + float or Quantity + p(vR,vz,R,z). + + Notes + ----- + - 2013-01-02 - Written - Bovy (IAS) + - 2018-01-12 - Added Gauss-Legendre integration prefactor vTmax/2 - Trick (MPA) """ if gl: if ngl % 2 == 1: @@ -3174,56 +2993,61 @@ def pvRvz(self, vR, vz, R, z, gl=True, ngl=_DEFAULTNGL2, vTmax=1.5): 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 + Calculate the epicycle frequency at r. + + Parameters + ---------- + r : float + Radius. + + Returns + ------- + float + Epicycle frequency. + + Notes + ----- + - 2012-07-25 - Written - Bovy (IAS@MPIA) """ 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 + Calculate the vertical frequency at r. + + Parameters + ---------- + r : float + Radius. + + Returns + ------- + float + Vertical frequency. + + Notes + ----- + - 2012-07-25 - Written - Bovy (IAS@MPIA) """ 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... + Calculate the radius of a circular orbit of Lz. + + Parameters + ---------- + lz : float + Angular momentum. + + Returns + ------- + float + Radius. + + Notes + ----- + - 2012-07-25 - Written - Bovy (IAS@MPIA) """ if isinstance(lz, numpy.ndarray): indx = (lz > self._precomputergLzmax) * (lz < self._precomputergLzmin) diff --git a/galpy/df/sphericaldf.py b/galpy/df/sphericaldf.py index 7ed327dd4..7179df040 100644 --- a/galpy/df/sphericaldf.py +++ b/galpy/df/sphericaldf.py @@ -41,34 +41,26 @@ class sphericaldf(df): def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initializes a spherical DF - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale= (None) length-scale parameter to be used internally - - ro= ,vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-07-22 - Written - Lane (UofT) - + Initializes a spherical DF + + Parameters + ---------- + pot : Potential instance or list thereof + The potential. Default is None. + denspot : Potential instance or list thereof, optional + The potential that represents the density of the tracers (assumed to be spherical). If None, set equal to pot. Default is None. + rmax : float or Quantity, optional + The maximum radius to consider. DF is cut off at E = Phi(rmax). Default is None. + scale : float or Quantity, optional + The length-scale parameter to be used internally. Default is None. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-07-22 - Written - Lane (UofT) """ df.__init__(self, ro=ro, vo=vo) if not conversion.physical_compatible(self, pot): @@ -113,34 +105,24 @@ def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=No @physical_conversion("phasespacedensity", pop=True) def __call__(self, *args, **kwargs): """ - NAME: - - __call__ - - PURPOSE: - - return the DF - - INPUT: + Evaluate the DF + Parameters + ---------- + *args: tuple Either: + a) (E,L,Lz): tuple of E and (optionally) L and (optionally) Lz. Each may be Quantity + b) R,vR,vT,z,vz,phi: cylindrical coordinates (can be Quantity) + c) Orbit instance: orbit.Orbit instance and if specific time then orbit.Orbit(t) - a) (E,L,Lz): tuple of E and (optionally) L and (optionally) Lz. - Each may be Quantity - - b) R,vR,vT,z,vz,phi: - - c) Orbit instance: orbit.Orbit instance and if specific time - then orbit.Orbit(t) - - OUTPUT: - + Returns + ------- + ndarray or Quantity Value of DF - HISTORY: - - 2020-07-22 - Written - Lane (UofT) - + Notes + ----- + - 2020-07-22 - Written - Lane (UofT) """ # Get E,L,Lz if len(args) == 1: @@ -181,25 +163,21 @@ def __call__(self, *args, **kwargs): @physical_conversion("energydensity", pop=True) def dMdE(self, E): """ - NAME: - - dMdE + Compute the differential energy distribution dM/dE: the amount of mass per unit energy - PURPOSE: + Parameters + ---------- + E : float or array_like + Energy; can be a Quantity - Compute the differential energy distribution dM/dE: the amount of mass per unit energy + Returns + ------- + float, numpy.ndarray, or Quantity + The differential energy distribution - INPUT: - - E - energy; can be a Quantity - - OUTPUT: - - dM/dE - - HISTORY: - - 2023-05-23 - Written - Bovy (UofT) + Notes + ----- + - 2023-05-23 - Written - Bovy (UofT) """ return self._dMdE( @@ -208,30 +186,25 @@ def dMdE(self, E): def vmomentdensity(self, r, n, m, **kwargs): """ - NAME: - - vmomentdensity - - PURPOSE: - - calculate an arbitrary moment of the velocity distribution - at r times the density - - INPUT: - - r - spherical radius at which to calculate the moment - - n - vr^n, where vr = v x cos eta - - m - vt^m, where vt = v x sin eta - - OUTPUT: - - at r - - HISTORY: - - 2020-09-04 - Written - Bovy (UofT) + Calculate an arbitrary moment of the velocity distribution at r times the density. + + Parameters + ---------- + r : float + Spherical radius at which to calculate the moment. + n : float + vr^n, where vr = v x cos eta. + m : float + vt^m, where vt = v x sin eta. + + Returns + ------- + float or Quantity + at r. + + Notes + ----- + - 2020-09-04 - Written - Bovy (UofT) """ r = conversion.parse_length(r, ro=self._ro) use_physical = kwargs.pop("use_physical", True) @@ -281,25 +254,21 @@ def _vmomentdensity(self, r, n, m): @physical_conversion("velocity", pop=True) def sigmar(self, r): """ - NAME: + Calculate the radial velocity dispersion at radius r. - sigmar + Parameters + ---------- + r : float + Spherical radius at which to calculate the radial velocity dispersion. - PURPOSE: + Returns + ------- + float or Quantity + The radial velocity dispersion at radius r. - calculate the radial velocity dispersion at radius r - - INPUT: - - r - spherical radius at which to calculate the radial velocity dispersion - - OUTPUT: - - sigma_r(r) - - HISTORY: - - 2020-09-04 - Written - Bovy (UofT) + Notes + ----- + - 2020-09-04 - Written - Bovy (UofT) """ r = conversion.parse_length(r, ro=self._ro) return numpy.sqrt(self._vmomentdensity(r, 2, 0) / self._vmomentdensity(r, 0, 0)) @@ -307,50 +276,44 @@ def sigmar(self, r): @physical_conversion("velocity", pop=True) def sigmat(self, r): """ - NAME: - - sigmar - - PURPOSE: + Calculate the tangential velocity dispersion at radius r. - calculate the tangential velocity dispersion at radius r + Parameters + ---------- + r : float + Spherical radius at which to calculate the tangential velocity dispersion. - INPUT: + Returns + ------- + float or Quantity + The tangential velocity dispersion at radius r. - r - spherical radius at which to calculate the tangential velocity dispersion + Notes + ----- + - 2020-09-04 - Written - Bovy (UofT) - OUTPUT: - - sigma_t(r) - - HISTORY: - - 2020-09-04 - Written - Bovy (UofT) """ r = conversion.parse_length(r, ro=self._ro) return numpy.sqrt(self._vmomentdensity(r, 0, 2) / self._vmomentdensity(r, 0, 0)) def beta(self, r): """ - NAME: - - sigmar - - PURPOSE: - - calculate the anisotropy at radius r - - INPUT: - - r - spherical radius at which to calculate the anisotropy + Calculate the anisotropy at radius r. - OUTPUT: + Parameters + ---------- + r : float + Spherical radius at which to calculate the anisotropy. - beta(r) + Returns + ------- + float + Anisotropy at radius r. - HISTORY: + Notes + ----- + - 2020-09-04 - Written - Bovy (UofT) - 2020-09-04 - Written - Bovy (UofT) """ r = conversion.parse_length(r, ro=self._ro) return 1.0 - self._vmomentdensity(r, 0, 2) / 2.0 / self._vmomentdensity(r, 2, 0) @@ -358,45 +321,32 @@ def beta(self, r): ############################### SAMPLING THE DF################################ def sample(self, R=None, z=None, phi=None, n=1, return_orbit=True, rmin=0.0): """ - NAME: - - sample - - PURPOSE: - - Return full 6D samples of the DF - - INPUT: - - R= cylindrical radius at which to generate samples (can be Quantity) - - z= height at which to generate samples (can be Quantity) - - phi= azimuth at which to generate samples (can be Quantity) - - n= number of samples to generate - - rmin= (0.) only sample r > rmin (can be Quantity) - - OPTIONAL INPUT: - - return_orbit= (True) If True output is an orbit.Orbit object, if False output is (R,vR,vT,z,vz,phi) - - OUTPUT: - - List of samples. Either vector (R,vR,vT,z,vz,phi) or orbit.Orbit; the (R,vR,vT,z,vz,phi) is either in internal units or is a set of Quantities - - NOTES: - - If R,z,phi are None then sample positions with CMF. If R,z,phi are - floats then sample n velocities at location. If array then sample - velocities at radii, ignoring n. phi can be None if R,z are set - by any above mechanism, will then sample phi for output. - - HISTORY: - - 2020-07-22 - Written - Lane (UofT) - + Sample the DF + + Parameters + ---------- + R : float, array_like, Quantity, or None, optional + If set, sample velocities at this radius. If array, sample velocities at these radii, ignoring n. + z : float, array_like, Quantity, or None, optional + If set, sample velocities at this height. If array, sample velocities at these heights, ignoring n. + phi : float, array_like, Quantity, or None, optional + If set, sample velocities at this azimuth. If array, sample velocities at these azimuths, ignoring n. + n : int, optional + Number of samples to generate. Default is 1. + return_orbit : bool, optional + If True, return an orbit.Orbit instance. If False, return a tuple of (R,vR,vT,z,vz,phi). Default is True. + rmin : float, Quantity, optional + Minimum radius at which to sample. Default is 0. + + Returns + ------- + orbit.Orbit instance or tuple + If return_orbit is True, an orbit.Orbit instance. Otherwise, a tuple of (R,vR,vT,z,vz,phi). + + Notes + ----- + - When specifying position, it is necessary to specify both R and z; if phi is not set in this case, it is sampled + - 2020-07-22 - Written - Lane (UofT) """ rmin = conversion.parse_length(rmin, ro=self._ro) if hasattr(self, "_rmin_sampling") and rmin != self._rmin_sampling: @@ -541,12 +491,6 @@ def _vmax_at_r(self, pot, r, **kwargs): def _make_pvr_interpolator(self, r_a_start=-3, r_a_end=3, n_r_a=120, n_v_vesc=100): """ - NAME: - - _make_pvr_interpolator - - PURPOSE: - Calculate a grid of the velocity sampling function v^2*f(E) over many radii. The radii are fractional with respect to some scale radius which characteristically describes the size of the potential, @@ -555,23 +499,25 @@ def _make_pvr_interpolator(self, r_a_start=-3, r_a_end=3, n_r_a=120, n_v_vesc=10 represents the inverse cumulative distribution at many radii. This allows for sampling of v/vesc given an input r/a - INPUT: - - r_a_start= radius grid start location in units of log10(r/a) - - r_a_end= radius grid end location in units of log10(r/a) - - n_r_a= number of radius grid points to use - - n_v_vesc= number of velocity grid points to use - - OUTPUT: - - None (But sets self._v_vesc_pvr_interpolator) - - HISTORY: - - Written 2020-07-24 - James Lane (UofT) + Parameters + ---------- + r_a_start : float, optional + Radius grid start location in units of log10(r/a). Default is -3. + r_a_end : float, optional + Radius grid end location in units of log10(r/a). Default is 3. + n_r_a : int, optional + Number of radius grid points to use. Default is 120. + n_v_vesc : int, optional + Number of velocity grid points to use. Default is 100. + + Returns + ------- + scipy.interpolate.RectBivariateSpline + Interpolator for v/vesc given an input r/a. + + Notes + ----- + - 2020-07-24 - Written - Lane (UofT) """ # Check that interpolated potential has appropriate grid range if ( @@ -640,30 +586,27 @@ def _make_pvr_interpolator(self, r_a_start=-3, r_a_end=3, n_r_a=120, n_v_vesc=10 def _setup_rphi_interpolator(self, r_a_min=1e-6, r_a_max=1e6, nra=10001): """ - NAME: - - _setup_rphi_interpolator - - PURPOSE: - Set up the interpolator for r(phi) - INPUT: - - r_a_min= minimum r/a - - r_a_max= maximum r/a - - nra= number of points to use in the r/a grid - - OUTPUT: - - _rphi_interpolator (scipy.interpolate.InterpolatedUnivariateSpline) - - HISTORY: - - Written 2023-02-23 - James Lane (UofT) + Parameters + ---------- + r_a_min : float, optional + Minimum r/a. Default is 1e-6. + r_a_max : float, optional + Maximum r/a. Default is 1e6. + nra : int, optional + Number of points to use in the r/a grid. Default is 10001. + + Returns + ------- + scipy.interpolate.InterpolatedUnivariateSpline + Interpolator for r(phi). + + Notes + ----- + - 2023-02-23 - Written - Lane (UofT) """ + r_a_values = numpy.concatenate( (numpy.array([0.0]), numpy.geomspace(r_a_min, r_a_max, nra)) ) @@ -686,33 +629,27 @@ class isotropicsphericaldf(sphericaldf): def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an isotropic distribution function - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale= scale parameter to be used internally - - ro=, vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-09-02 - Written - Bovy (UofT) + Initialize an isotropic distribution function + + Parameters + ---------- + pot : Potential instance or list thereof + Default: None + denspot : Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot), optional + Default: None + rmax : float or Quantity, optional + Maximum radius to consider; DF is cut off at E = Phi(rmax) + Default: None + scale : float, optional + Scale parameter to be used internally + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-09-02 - Written - Bovy (UofT) """ sphericaldf.__init__( @@ -721,25 +658,20 @@ def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=No def _call_internal(self, *args): """ - NAME: - - _call_internal - - PURPOSE - - Calculate the distribution function for an isotropic DF + Calculate the distribution function for an isotropic DF. - INPUT: + Parameters + ---------- + *args : tuple of (E,L,Lz) with L and Lz optionalA - E,L,Lz - The energy, angular momemtum magnitude, and its z component (only E is used) + Returns + ------- + float + The distribution function evaluated at E. - OUTPUT: - - f(x,v) = f(E[x,v]) - - HISTORY: - - 2020-07 - Written - Lane (UofT) + Notes + ----- + - 2020-07 - Written - Lane (UofT) """ return self.fE(args[0]) @@ -811,33 +743,26 @@ class anisotropicsphericaldf(sphericaldf): def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=None): """ - NAME: - - __init__ - - PURPOSE: - - Initialize an anisotropic distribution function - - INPUT: - - pot= (None) Potential instance or list thereof - - denspot= (None) Potential instance or list thereof that represent the density of the tracers (assumed to be spherical; if None, set equal to pot) - - rmax= (None) maximum radius to consider (can be Quantity); DF is cut off at E = Phi(rmax) - - scale= (None) length-scale parameter to be used internally - - ro= ,vo= galpy unit parameters - - OUTPUT: - - None - - HISTORY: - - 2020-07-22 - Written - Lane (UofT) + Initialize an anisotropic distribution function + + Parameters + ---------- + pot : Potential instance or list thereof + The potential. Default: None. + denspot : Potential instance or list thereof, optional + The potential representing the density of the tracers (assumed to be spherical). If None, set equal to pot. Default: None. + rmax : float or Quantity, optional + Maximum radius to consider. DF is cut off at E = Phi(rmax). Default: None. + scale : float, optional + Length-scale parameter to be used internally. Default: None. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2020-07-22 - Written - Lane (UofT) """ sphericaldf.__init__( diff --git a/galpy/df/streamdf.py b/galpy/df/streamdf.py index 5e3cba9aa..72103e940 100644 --- a/galpy/df/streamdf.py +++ b/galpy/df/streamdf.py @@ -95,91 +95,71 @@ def __init__( custom_transform=None, ): """ - NAME: - - __init__ - - PURPOSE: - - Initialize a quasi-isothermal DF - - INPUT: - - sigv - radial velocity dispersion of the progenitor (can be Quantity) - - tdisrupt= (5 Gyr) time since start of disruption (can be Quantity) - - leading= (True) if True, model the leading part of the stream - if False, model the trailing part - - progenitor= progenitor orbit as Orbit instance (will be re-integrated, so don't bother integrating the orbit before) - - progIsTrack= (False) if True, then the progenitor (x,v) is actually the (x,v) of the stream track at zero angle separation; useful when initializing with an orbit fit; the progenitor's position will be calculated - - pot= Potential instance or list thereof - - aA= actionAngle instance used to convert (x,v) to actions - - useTM= (False) if set to an actionAngleTorus instance, use this to speed up calculations - - sigMeanOffset= (6.) offset between the mean of the frequencies - and the progenitor, in units of the largest - eigenvalue of the frequency covariance matrix - (along the largest eigenvector), should be positive; - to model the trailing part, set leading=False - - sigangle= (sigv/122/[1km/s]=1.8sigv in natural coordinates) - estimate of the angle spread of the debris initially (can be Quantity) - - deltaAngleTrack= (None) angle to estimate the stream track over (rad; or can be Quantity) - - nTrackChunks= (floor(deltaAngleTrack/0.15)+1) number of chunks to divide the progenitor track in - - nTrackIterations= Number of iterations to perform when establishing the track; each iteration starts from a previous approximation to the track in (x,v) and calculates a new track based on the deviation between the previous track and the desired track in action-angle coordinates; if not set, an appropriate value is determined based on the magnitude of the misalignment between stream and orbit, with larger numbers of iterations for larger misalignments - - interpTrack= (might change), interpolate the stream track while - setting up the instance (can be done by hand by - calling self._interpolate_stream_track() and - self._interpolate_stream_track_aA()) - - useInterp= (might change), use interpolation by default when - calculating approximated frequencies and angles - - nosetup= (False) if True, don't setup the stream track and anything - else that is expensive - - nospreadsetup= (False) if True, don't setup the spread around the stream track (only for nosetup is False) - - multi= (None) if set, use multi-processing - - Coordinate transformation inputs: - - vo= (220) circular velocity to normalize velocities with [used to be Vnorm; can be Quantity] - - ro= (8) Galactocentric radius to normalize positions with [used to be Rnorm; can be Quantity] - - R0= (8) Galactocentric radius of the Sun (kpc) [can be different from ro; can be Quantity] - - Zsun= (0.0208) Sun's height above the plane (kpc; can be Quantity) - - vsun= ([-11.1,241.92,7.25]) Sun's motion in cylindrical coordinates (vR positive away from center) (can be Quantity array, but not a list of Quantities) - - custom_transform= (None) matrix implementing the rotation from (ra,dec) to a custom set of sky coordinates - - approxConstTrackFreq= (False) if True, approximate the stream assuming that the frequency is constant along the stream (only works with useTM, for which this leads to a significant speed-up) - - useTMHessian= (False) if True, compute the basic Hessian dO/dJ_prog using TM; otherwise use aA - - OUTPUT: - - object - - HISTORY: - - 2013-09-16 - Started - Bovy (IAS) - - 2013-11-25 - Started over - Bovy (IAS) - + Initialize the DF of a tidal stream + + Parameters + ---------- + sigv : float or Quantity + Radial velocity dispersion of the progenitor. + progenitor : galpy.orbit.Orbit + Progenitor orbit as Orbit instance (will be re-integrated, so don't bother integrating the orbit before). + pot : galpy.potential.Potential or list thereof, optional + Potential instance or list thereof. + aA : actionAngle instance + ActionAngle instance used to convert (x,v) to actions. Generally a actionAngleIsochroneApprox instance. + useTM : bool, optional + If set to an actionAngleTorus instance, use this to speed up calculations. + tdisrupt : float or Quantity, optional + Time since start of disruption (default: 5 Gyr). + sigMeanOffset : float, optional + Offset between the mean of the frequencies and the progenitor, in units of the largest eigenvalue of the frequency covariance matrix (along the largest eigenvector), should be positive; to model the trailing part, set leading=False (default: 6.0). + leading : bool, optional + If True, model the leading part of the stream; if False, model the trailing part (default: True). + sigangle : float or Quantity, optional + Estimate of the angle spread of the debris initially (default: sigv/122/[1km/s]=1.8sigv in natural coordinates). + deltaAngleTrack : float or Quantity, optional + Angle to estimate the stream track over (rad; or can be Quantity) (default: None). + nTrackChunks : int, optional + Number of chunks to divide the progenitor track in (default: floor(deltaAngleTrack/0.15)+1). + nTrackIterations : int, optional + Number of iterations to perform when establishing the track; each iteration starts from a previous approximation to the track in (x,v) and calculates a new track based on the deviation between the previous track and the desired track in action-angle coordinates; if not set, an appropriate value is determined based on the magnitude of the misalignment between stream and orbit, with larger numbers of iterations for larger misalignments (default: None). + progIsTrack : bool, optional + If True, then the progenitor (x,v) is actually the (x,v) of the stream track at zero angle separation; useful when initializing with an orbit fit; the progenitor's position will be calculated (default: False). + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + Vnorm : float or Quantity, optional + Deprecated. Use vo instead (default: None). + Rnorm : float or Quantity, optional + Deprecated. Use ro instead (default: None). + R0 : float or Quantity, optional + Galactocentric radius of the Sun (kpc) (can be different from ro) (default: 8.0). + Zsun : float or Quantity, optional + Sun's height above the plane (kpc) (default: 0.0208). + vsun : array_like or Quantity, optional + Sun's motion in cylindrical coordinates (vR positive away from center) (can be Quantity array, but not a list of Quantities) (default: [-11.1, 8.0 * 30.24, 7.25]). + multi : int, optional + If set, use multi-processing (default: None). + interpTrack : bool, optional + Interpolate the stream track while setting up the instance (can be done by hand by calling self._interpolate_stream_track() and self._interpolate_stream_track_aA()) (default: _INTERPDURINGSETUP). + useInterp : bool, optional + Use interpolation by default when calculating approximated frequencies and angles (default: _USEINTERP). + nosetup : bool, optional + If True, don't setup the stream track and anything else that is expensive (default: False). + nospreadsetup : bool, optional + If True, don't setup the spread around the stream track (only for nosetup is False) (default: False). + approxConstTrackFreq : bool, optional + If True, approximate the stream assuming that the frequency is constant along the stream (only works with useTM, for which this leads to a significant speed-up) (default: False). + useTMHessian : bool, optional + If True, compute the basic Hessian dO/dJ_prog using TM; otherwise use aA (default: False). + custom_transform : array_like, optional + Matrix implementing the rotation from (ra,dec) to a custom set of sky coordinates (default: None). + + Notes + ----- + - 2013-09-16 - Started - Bovy (IAS) + - 2013-11-25 - Started over - Bovy (IAS) """ if ro is None and not Rnorm is None: warnings.warn( @@ -442,28 +422,27 @@ def _setup_progIsTrack(self): @physical_conversion("angle", pop=True) def misalignment(self, isotropic=False, **kwargs): """ - NAME: - - misalignment - - PURPOSE: - - calculate the misalignment between the progenitor's frequency - and the direction along which the stream disrupts - - INPUT: + Calculate the misalignment between the progenitor's frequency + and the direction along which the stream disrupts. - isotropic= (False), if True, return the misalignment assuming an isotropic action distribution + Parameters + ---------- + isotropic : bool, optional + If True, return the misalignment assuming an isotropic action distribution. - OUTPUT: + Returns + ------- + float + Misalignment in radians. - misalignment in rad + Warnings + -------- + In versions >1.3, the output unit of streamdf.misalignment has been changed to radian (from degree before). - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) - - 2017-10-28 - Changed output unit to rad - Bovy (UofT) + Notes + ----- + - 2013-12-05 - Written - Bovy (IAS) + - 2017-10-28 - Changed output unit to rad - Bovy (UofT) """ warnings.warn( @@ -485,27 +464,21 @@ def misalignment(self, isotropic=False, **kwargs): def freqEigvalRatio(self, isotropic=False): """ - NAME: - - freqEigvalRatio - - PURPOSE: - - calculate the ratio between the largest and 2nd-to-largest (in abs) - eigenvalue of sqrt(dO/dJ^T V_J dO/dJ) - (if this is big, a 1D stream will form) + Calculate the ratio between the largest and 2nd-to-largest (in abs) eigenvalue of sqrt(dO/dJ^T V_J dO/dJ) (if this is big, a 1D stream will form) - INPUT: + Parameters + ---------- + isotropic : bool, optional + If True, return the ratio assuming an isotropic action distribution (i.e., just of dO/dJ) - isotropic= (False), if True, return the ratio assuming an isotropic action distribution (i.e., just of dO/dJ) + Returns + ------- + float + Ratio between eigenvalues of fabs(dO / dJ) - OUTPUT: - - ratio between eigenvalues of fabs(dO / dJ) - - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-05 - Written - Bovy (IAS) """ if isotropic: @@ -519,25 +492,21 @@ def freqEigvalRatio(self, isotropic=False): @physical_conversion("time", pop=True) def estimateTdisrupt(self, deltaAngle): """ - NAME: - - estimateTdisrupt - - PURPOSE: + Estimate the time of disruption. - estimate the time of disruption + Parameters + ---------- + deltaAngle : float + Spread in angle since disruption. - INPUT: + Returns + ------- + float or Quantity + Time of disruption. - deltaAngle- spread in angle since disruption - - OUTPUT: - - time in natural units - - HISTORY: - - 2013-11-27 - Written - Bovy (IAS) + Notes + ----- + - 2013-11-27 - Written - Bovy (IAS) """ return deltaAngle / numpy.sqrt(numpy.sum(self._dsigomeanProg**2.0)) @@ -546,33 +515,29 @@ def subhalo_encounters( self, venc=numpy.inf, sigma=150.0 / 220.0, nsubhalo=0.3, bmax=0.025, yoon=False ): """ - NAME: - - subhalo_encounters - - PURPOSE: - - estimate the number of encounters with subhalos over the lifetime of this stream, using a formalism similar to that of Yoon et al. (2011) - - INPUT: - - venc= (numpy.inf) count encounters with (relative) speeds less than this (relative radial velocity in cylindrical stream frame, unless yoon is True) (can be Quantity) - - sigma= (150/220) velocity dispersion of the DM subhalo population (can be Quantity) - - nsubhalo= (0.3) spatial number density of subhalos (can be Quantity) - - bmax= (0.025) maximum impact parameter (if larger than width of stream) (can be Quantity) - - yoon= (False) if True, use erroneous Yoon et al. formula - - OUTPUT: - - number of encounters - - HISTORY: - - 2016-01-19 - Written - Bovy (UofT) + Estimate the number of encounters with subhalos over the lifetime of this stream, using a formalism similar to that of Yoon et al. (2011) + + Parameters + ---------- + venc : float, optional + Count encounters with (relative) speeds less than this (relative radial velocity in cylindrical stream frame, unless yoon is True) (can be Quantity) + sigma : float, optional + Velocity dispersion of the DM subhalo population (can be Quantity) + nsubhalo : float, optional + Spatial number density of subhalos (can be Quantity) + bmax : float, optional + Maximum impact parameter (if larger than width of stream) (can be Quantity) + yoon : bool, optional + If True, use erroneous Yoon et al. formula + + Returns + ------- + float + Number of encounters + + Notes + ----- + - 2016-01-19 - Written - Bovy (UofT) """ venc = conversion.parse_velocity(venc, vo=self._vo) @@ -621,37 +586,32 @@ def plotTrack( self, d1="x", d2="z", interp=True, spread=0, simple=_USESIMPLE, *args, **kwargs ): """ - NAME: - - plotTrack - - PURPOSE: - - plot the stream track - - INPUT: - - d1= plot this on the X axis ('x','y','z','R','phi','vx','vy','vz','vR','vt','ll','bb','dist','pmll','pmbb','vlos') - - d2= plot this on the Y axis (same list as for d1) - - interp= (True) if True, use the interpolated stream track - - spread= (0) if int > 0, also plot the spread around the track as spread x sigma - - scaleToPhysical= (False), if True, plot positions in kpc and velocities in km/s - - simple= (False), if True, use a simple estimate for the spread in perpendicular angle - - galpy.util.plot.plotplot args and kwargs - - OUTPUT: - - plot to output device - - HISTORY: - - 2013-12-09 - Written - Bovy (IAS) + Plot the stream track. + + Parameters + ---------- + d1 : str, optional + Plot this on the X axis ('x','y','z','R','phi','vx','vy','vz','vR','vt','ll','bb','dist','pmll','pmbb','vlos'). Default is 'x'. + d2 : str, optional + Plot this on the Y axis (same list as for d1). Default is 'z'. + interp : bool, optional + If True, use the interpolated stream track. Default is True. + spread : int, optional + If int > 0, also plot the spread around the track as spread x sigma. Default is 0. + simple : bool, optional + If True, use a simple estimate for the spread in perpendicular angle. Default is _USESIMPLE. + *args : tuple + Arguments passed to galpy.util.plot.plot. + **kwargs : dict + Keyword arguments passed to galpy.util.plot.plot. + + Returns + ------- + None + + Notes + ----- + - 2013-12-09 - Written - Bovy (IAS) """ if not hasattr(self, "_ObsTrackLB") and ( @@ -717,32 +677,28 @@ def plotTrack( def plotProgenitor(self, d1="x", d2="z", *args, **kwargs): """ - NAME: - - plotProgenitor - - PURPOSE: - - plot the progenitor orbit - - INPUT: - - d1= plot this on the X axis ('x','y','z','R','phi','vx','vy','vz','vR','vt','ll','bb','dist','pmll','pmbb','vlos') - - d2= plot this on the Y axis (same list as for d1) - - scaleToPhysical= (False), if True, plot positions in kpc and velocities in km/s - - galpy.util.plot.plot args and kwargs - - OUTPUT: - - plot to output device - - HISTORY: - - 2013-12-09 - Written - Bovy (IAS) - + Plot the progenitor orbit. + + Parameters + ---------- + d1 : str, optional + Plot this on the X axis ('x','y','z','R','phi','vx','vy','vz','vR','vt','ll','bb','dist','pmll','pmbb','vlos'). Default is 'x'. + d2 : str, optional + Plot this on the Y axis (same list as for d1). Default is 'z'. + scaleToPhysical : bool, optional + If True, plot positions in kpc and velocities in km/s. Default is False. + *args : tuple + Arguments passed to `galpy.util.plot.plot`. + **kwargs : dict + Keyword arguments passed to `galpy.util.plot.plot`. + + Returns + ------- + None + + Notes + ----- + - 2013-12-09 - Written - Bovy (IAS) """ tts = self._progenitor.t[ self._progenitor.t < self._trackts[self._nTrackChunks - 1] @@ -1013,26 +969,20 @@ def _parse_track_spread(self, d1, d2, interp=True, phys=False, simple=_USESIMPLE def plotCompareTrackAAModel(self, **kwargs): """ - NAME: - - plotCompareTrackAAModel - - PURPOSE: - - plot the comparison between the underlying model's dOmega_perp vs. dangle_r (line) and the track in (x,v)'s dOmega_perp vs. dangle_r (dots; explicitly calculating the track's action-angle coordinates) - - INPUT: - - galpy.util.plot.plot kwargs - - OUTPUT: + Plot the comparison between the underlying model's dOmega_perp vs. dangle_r (line) and the track in (x,v)'s dOmega_perp vs. dangle_r (dots; explicitly calculating the track's action-angle coordinates) - plot + Parameters + ---------- + **kwargs: dict + galpy.util.plot.plot kwargs - HISTORY: - - 2014-08-27 - Written - Bovy (IAS) + Returns + ------- + None + Notes + ----- + - 2014-08-27 - Written - Bovy (IAS) """ # First calculate the model model_adiff = (self._ObsTrackAA[:, 3:] - self._progenitor_angle)[ @@ -1766,37 +1716,28 @@ def _interpolate_stream_track_aA(self): def calc_stream_lb(self, vo=None, ro=None, R0=None, Zsun=None, vsun=None): """ - NAME: - - calc_stream_lb - - PURPOSE: - - convert the stream track to observational coordinates and store - - INPUT: - - Coordinate transformation inputs (all default to the instance-wide - values): - - vo= circular velocity to normalize velocities with - - ro= Galactocentric radius to normalize positions with - - R0= Galactocentric radius of the Sun (kpc) - - Zsun= Sun's height above the plane (kpc) - - vsun= Sun's motion in cylindrical coordinates (vR positive away from center) - - OUTPUT: - - (none) - - HISTORY: - - 2013-12-02 - Written - Bovy (IAS) - + Convert the stream track to observational coordinates and store + + Parameters + ---------- + ro : float or Quantity, optional + Distance scale for translation into internal units (object-wide default. + vo : float or Quantity, optional + Velocity scale for translation into internal units (object-wide default. + R0 : float or Quantity, optional + Distance to the Galactic center (object-wide default. + Zsun : float or Quantity, optional + Distance to the Galactic plane (object-wide default. + vsun : float or Quantity, optional + Velocity of the Sun (object-wide default. + + Returns + ------- + None + + Notes + ----- + - 2013-12-02 - Written - Bovy (IAS) """ if vo is None: vo = self._vo @@ -1890,32 +1831,27 @@ def find_closest_trackpoint( self, R, vR, vT, z, vz, phi, interp=True, xy=False, usev=False ): """ - NAME: - - find_closest_trackpoint - - PURPOSE: - - find the closest point on the stream track to a given point - - INPUT: - - R,vR,vT,z,vz,phi - phase-space coordinates of the given point - - interp= (True), if True, return the index of the interpolated track - - xy= (False) if True, input is X,Y,Z,vX,vY,vZ in Galactocentric rectangular coordinates; if xy, some coordinates may be missing (given as None) and they will not be used - - usev= (False) if True, also use velocities to find the closest point - - OUTPUT: - - index into the track of the closest track point - - HISTORY: - - 2013-12-04 - Written - Bovy (IAS) - + Find the closest point on the stream track to a given point + + Parameters + ---------- + R,vR,vT,z,vz,phi : float + Phase-space coordinates of the given point + interp : bool, optional + If True, return the index of the interpolated track + xy : bool, optional + If True, input is X,Y,Z,vX,vY,vZ in Galactocentric rectangular coordinates; if xy, some coordinates may be missing (given as None) and they will not be used + usev : bool, optional + If True, also use velocities to find the closest point + + Returns + ------- + int + Index into the track of the closest track point + + Notes + ----- + - 2013-12-04 - Written - Bovy (IAS) """ if xy: X = R @@ -1986,28 +1922,35 @@ def find_closest_trackpointLB( self, l, b, D, vlos, pmll, pmbb, interp=True, usev=False ): """ - NAME: - - find_closest_trackpointLB - - PURPOSE: - find the closest point on the stream track to a given point in (l,b,...) coordinates - - INPUT: - - l,b,D,vlos,pmll,pmbb- coordinates in (deg,deg,kpc,km/s,mas/yr,mas/yr) - - interp= (True) if True, return the closest index on the interpolated track - - usev= (False) if True, also use the velocity components (default is to only use the positions) - - OUTPUT: - - index of closest track point on the interpolated or not-interpolated track - - HISTORY: - - 2013-12-17- Written - Bovy (IAS) + Find the closest point on the stream track to a given point in (l,b,...) coordinates + + Parameters + ---------- + l : float + Galactic longitude in degrees + b : float + Galactic latitude in degrees + D : float + Distance in kpc + vlos : float + Line-of-sight velocity in km/s + pmll : float + Proper motion in Galactic longitude in mas/yr + pmbb : float + Proper motion in Galactic latitude in mas/yr + interp : bool, optional + If True, return the closest index on the interpolated track (default is True) + usev : bool, optional + If True, also use the velocity components (default is False) + + Returns + ------- + int + Index of closest track point on the interpolated or not-interpolated track + + Notes + ----- + - 2013-12-17 - Written - Bovy (IAS) """ if interp: @@ -2089,18 +2032,34 @@ def find_closest_trackpointLB( def _find_closest_trackpointaA(self, Or, Op, Oz, ar, ap, az, interp=True): """ - NAME: - _find_closest_trackpointaA - PURPOSE: - find the closest point on the stream track to a given point in - frequency-angle coordinates - INPUT: - Or,Op,Oz,ar,ap,az - phase-space coordinates of the given point - interp= (True), if True, return the index of the interpolated track - OUTPUT: - index into the track of the closest track point - HISTORY: - 2013-12-22 - Written - Bovy (IAS) + Find the closest point on the stream track to a given point in + frequency-angle coordinates + + Parameters + ---------- + Or : float + Radial frequency + Op : float + Azimuthal frequency + Oz : float + Vertical frequency + ar : float + Radial angle + ap : float + Azimuthal angle + az : float + Vertical angle + interp : bool, optional + If True, return the index of the interpolated track (default is True) + + Returns + ------- + int + Index into the track of the closest track point + + Notes + ----- + - 2013-12-22 - Written - Bovy (IAS) """ # Calculate angle offset along the stream parallel to the stream track, # finding first the angle among a few wraps where the point is @@ -2133,29 +2092,27 @@ def _find_closest_trackpointaA(self, Or, Op, Oz, ar, ap, az, interp=True): #########DISTRIBUTION AS A FUNCTION OF ANGLE ALONG THE STREAM################## def pOparapar(self, Opar, apar, tdisrupt=None): """ - NAME: - - pOparapar - - PURPOSE: - - return the probability of a given parallel (frequency,angle) offset pair - - INPUT: - - Opar - parallel frequency offset (array) (can be Quantity) - - apar - parallel angle offset along the stream (scalar) (can be Quantity) - - OUTPUT: - - p(Opar,apar) - - HISTORY: - - 2015-12-07 - Written - Bovy (UofT) - + Return the probability of a given parallel (frequency,angle) offset pair. + + Parameters + ---------- + Opar : numpy.ndarray or Quantity + Parallel frequency offset. + apar : float or Quantity + Parallel angle offset along the stream. + tdisrupt : float, optional + Time since the start of the simulation at which the stream is disrupted (in galpy time units). If None, the value set at the initialization of the object is used. + + Returns + ------- + ndarray + Probability of a given parallel (frequency,angle) offset pair. + + Notes + ----- + - 2015-11-17 - Written - Bovy (UofT). """ + Opar = conversion.parse_frequency(Opar, ro=self._ro, vo=self._vo) apar = conversion.parse_angle(apar) if tdisrupt is None: @@ -2174,28 +2131,27 @@ def pOparapar(self, Opar, apar, tdisrupt=None): def density_par(self, dangle, coord="apar", tdisrupt=None, **kwargs): """ - NAME: - - density_par - - PURPOSE: - - calculate the density as a function of a parallel coordinate - - INPUT: - - dangle - parallel angle offset for this coordinate value - - coord - coordinate to return the density in ('apar' [default], - 'll','ra','customra','phi') - - OUTPUT: - - density(angle) - - HISTORY: - - 2015-11-17 - Written - Bovy (UofT) + Calculate the density as a function of a parallel coordinate + + Parameters + ---------- + dangle : float + Parallel angle offset for this coordinate value + coord : str, optional + Coordinate to return the density in ('apar' [default], + 'll','ra','customra','phi') + tdisrupt : float, optional + Time since the disruption to calculate the density at (in galpy + natural units). If None, use the current time (default: None). + + Returns + ------- + float + Density(angle) + + Notes + ----- + - 2015-11-17 - Written - Bovy (UofT) """ if coord.lower() != "apar": @@ -2276,32 +2232,27 @@ def _density_par(self, dangle, tdisrupt=None): def length(self, threshold=0.2, phys=False, ang=False, tdisrupt=None, **kwargs): """ - NAME: - - length - - PURPOSE: - - calculate the length of the stream - - INPUT: - - threshold - threshold down from the density near the progenitor at which to define the 'end' of the stream - - phys= (False) if True, return the length in physical kpc - - ang= (False) if True, return the length in sky angular arc length in degree - - coord - coordinate to return the density in ('apar' [default], - 'll','ra','customra','phi') - - OUTPUT: - - length (rad for parallel angle; kpc for physical length; deg for sky arc length) - - HISTORY: - - 2015-12-22 - Written - Bovy (UofT) + Calculate the length of the stream. + + Parameters + ---------- + threshold : float, optional + Threshold down from the density near the progenitor at which to define the 'end' of the stream. Default is 0.2. + phys : bool, optional + If True, return the length in physical kpc. Default is False. + ang : bool, optional + If True, return the length in sky angular arc length in degree. Default is False. + tdisrupt : float, optional + Time since disruption in natural units. Default is None. + + Returns + ------- + float + Length (rad for parallel angle; kpc for physical length; deg for sky arc length). + + Notes + ----- + - 2015-12-22 - Written - Bovy (UofT) """ peak_dens = self.density_par( @@ -2405,29 +2356,27 @@ def length(self, threshold=0.2, phys=False, ang=False, tdisrupt=None, **kwargs): @physical_conversion("frequency", pop=True) def meanOmega(self, dangle, oned=False, offset_sign=None, tdisrupt=None): """ - NAME: - - meanOmega - - PURPOSE: - - calculate the mean frequency as a function of angle, assuming a uniform time distribution up to a maximum time - - INPUT: - - dangle - angle offset - - oned= (False) if True, return the 1D offset from the progenitor (along the direction of disruption) - - offset_sign= sign of the frequency offset (shouldn't be set) - - OUTPUT: - - mean Omega - - HISTORY: - - 2013-12-01 - Written - Bovy (IAS) + Calculate the mean frequency as a function of angle, assuming a uniform time distribution up to a maximum time. + + Parameters + ---------- + dangle : float + Angle offset. + oned : bool, optional + If True, return the 1D offset from the progenitor (along the direction of disruption). Default is False. + offset_sign : None, optional + Sign of the frequency offset. Default is None. + tdisrupt : float, optional + Maximum time. Default is None. + + Returns + ------- + float or Quantity + Mean Omega. + + Notes + ----- + - 2013-12-01 - Written - Bovy (IAS) """ if offset_sign is None: @@ -2458,25 +2407,21 @@ def meanOmega(self, dangle, oned=False, offset_sign=None, tdisrupt=None): @physical_conversion("frequency", pop=True) def sigOmega(self, dangle): """ - NAME: - - sigmaOmega - - PURPOSE: + Calculate the 1D sigma in frequency as a function of angle, assuming a uniform time distribution up to a maximum time - calculate the 1D sigma in frequency as a function of angle, assuming a uniform time distribution up to a maximum time + Parameters + ---------- + dangle : float + Angle offset - INPUT: + Returns + ------- + float or Quantity + Sigma Omega - dangle - angle offset - - OUTPUT: - - sigma Omega - - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-01 - Written - Bovy (IAS) """ dOmin = dangle / self._tdisrupt @@ -2502,27 +2447,23 @@ def sigOmega(self, dangle): def ptdAngle(self, t, dangle): """ - NAME: - - ptdangle + Return the probability of a given stripping time at a given angle along the stream. - PURPOSE: + Parameters + ---------- + t : numpy.ndarray + Stripping time. + dangle : float + Angle offset along the stream. - return the probability of a given stripping time at a given angle along the stream + Returns + ------- + numpy.ndarray + p(td|dangle) - INPUT: - - t - stripping time - - dangle - angle offset along the stream - - OUTPUT: - - p(td|dangle) - - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-05 - Written - Bovy (IAS). """ t = numpy.array(t) @@ -2540,25 +2481,21 @@ def ptdAngle(self, t, dangle): @physical_conversion("time", pop=True) def meantdAngle(self, dangle): """ - NAME: - - meantdAngle - - PURPOSE: - - calculate the mean stripping time at a given angle - - INPUT: + Calculate the mean stripping time at a given angle. - dangle - angle offset along the stream + Parameters + ---------- + dangle : float + Angle offset along the stream. - OUTPUT: + Returns + ------- + float or Quantity + Mean stripping time at this dangle. - mean stripping time at this dangle - - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-05 - Written - Bovy (IAS) """ Tlow = dangle / (self._meandO + 3.0 * numpy.sqrt(self._sortedSigOEig[2])) @@ -2577,25 +2514,21 @@ def meantdAngle(self, dangle): @physical_conversion("time", pop=True) def sigtdAngle(self, dangle): """ - NAME: - - sigtdAngle - - PURPOSE: - - calculate the dispersion in the stripping times at a given angle + Calculate the dispersion in the stripping times at a given angle. - INPUT: + Parameters + ---------- + dangle : float + Angle offset along the stream. - dangle - angle offset along the stream + Returns + ------- + float or Quantity + Dispersion in the stripping times at this angle. - OUTPUT: - - dispersion in the stripping times at this angle - - HISTORY: - - 2013-12-05 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-05 - Written - Bovy (IAS) """ Tlow = dangle / (self._meandO + 3.0 * numpy.sqrt(self._sortedSigOEig[2])) @@ -2612,28 +2545,25 @@ def sigtdAngle(self, dangle): def pangledAngle(self, angleperp, dangle, smallest=False): """ - NAME: - - pangledAngle - - PURPOSE: - return the probability of a given perpendicular angle at a given angle along the stream - - INPUT: - - angleperp - perpendicular angle - - dangle - angle offset along the stream - - smallest= (False) calculate for smallest eigenvalue direction rather than for middle - - OUTPUT: - - p(angle_perp|dangle) - - HISTORY: - - 2013-12-06 - Written - Bovy (IAS) + Return the probability of a given perpendicular angle at a given angle along the stream. + + Parameters + ---------- + angleperp : numpy.ndarray + Perpendicular angle. + dangle : float + Angle offset along the stream. + smallest : bool, optional + Calculate for smallest eigenvalue direction rather than for middle. Default is False. + + Returns + ------- + numpy.ndarray + p(angle_perp|dangle) + + Notes + ----- + - 2013-12-06 - Written - Bovy (IAS) """ angleperp = numpy.array(angleperp) @@ -2651,27 +2581,23 @@ def pangledAngle(self, angleperp, dangle, smallest=False): @physical_conversion("angle", pop=True) def meanangledAngle(self, dangle, smallest=False): """ - NAME: - - meanangledAngle - - PURPOSE: + Calculate the mean perpendicular angle at a given angle. - calculate the mean perpendicular angle at a given angle + Parameters + ---------- + dangle : float + Angle offset along the stream. + smallest : bool, optional + Calculate for smallest eigenvalue direction rather than for middle. Default is False. - INPUT: + Returns + ------- + float or Quantity + Mean perpendicular angle. - dangle - angle offset along the stream - - smallest= (False) calculate for smallest eigenvalue direction rather than for middle - - OUTPUT: - - mean perpendicular angle - - HISTORY: - - 2013-12-06 - Written - Bovy (IAS) + Notes + ----- + - 2013-12-06 - Written - Bovy (IAS) """ if smallest: @@ -2696,32 +2622,27 @@ def meanangledAngle(self, dangle, smallest=False): @physical_conversion("angle", pop=True) def sigangledAngle(self, dangle, assumeZeroMean=True, smallest=False, simple=False): """ - NAME: - - sigangledAngle - - PURPOSE: - - calculate the dispersion in the perpendicular angle at a given angle - - INPUT: - - dangle - angle offset along the stream - - assumeZeroMean= (True) if True, assume that the mean is zero (should be) - - smallest= (False) calculate for smallest eigenvalue direction rather than for middle - - simple= (False), if True, return an even simpler estimate - - OUTPUT: - - dispersion in the perpendicular angle at this angle - - HISTORY: - - 2013-12-06 - Written - Bovy (IAS) - + Calculate the dispersion in the perpendicular angle at a given angle. + + Parameters + ---------- + dangle : float + Angle offset along the stream. + assumeZeroMean : bool, optional + If True, assume that the mean is zero (should be). Default is True. + smallest : bool, optional + Calculate for smallest eigenvalue direction rather than for middle. Default is False. + simple : bool, optional + If True, return an even simpler estimate. Default is False. + + Returns + ------- + float or Quantity + Dispersion in the perpendicular angle at this angle. + + Notes + ----- + - 2013-12-06 - Written - Bovy (IAS) """ if smallest: eigIndx = 0 @@ -2779,20 +2700,36 @@ def _pangledAnglet(self, t, angleperp, dangle, smallest): ################APPROXIMATE FREQUENCY-ANGLE TRANSFORMATION##################### def _approxaA(self, R, vR, vT, z, vz, phi, interp=True, cindx=None): """ - NAME: - _approxaA - PURPOSE: - return action-angle coordinates for a point based on the linear - approximation around the stream track - INPUT: - R,vR,vT,z,vz,phi - phase-space coordinates of the given point - interp= (True), if True, use the interpolated track - cindx= index of the closest point on the (interpolated) stream track if not given, determined from the dimensions given - OUTPUT: - (Or,Op,Oz,ar,ap,az) - HISTORY: - 2013-12-03 - Written - Bovy (IAS) - 2015-11-12 - Added weighted sum of two nearest Jacobians to help with smoothness - Bovy (UofT) + Return action-angle coordinates for a point based on the linear approximation around the stream track + + Parameters + ---------- + R : float + Radius + vR : float + Radial velocity + vT : float + Tangential velocity + z : float + Vertical height + vz : float + Vertical velocity + phi : float + Azimuth + interp : bool, optional + If True, use the interpolated track. Default is True. + cindx : int, optional + Index of the closest point on the (interpolated) stream track. If not given, determined from the dimensions given. + + Returns + ------- + tuple + (Or,Op,Oz,ar,ap,az) + + Notes + ----- + - 2013-12-03 - Written - Bovy (IAS) + - 2015-11-12 - Added weighted sum of two nearest Jacobians to help with smoothness - Bovy (UofT) """ if isinstance(R, (int, float, numpy.float32, numpy.float64)): # Scalar input R = numpy.array([R]) @@ -2905,19 +2842,33 @@ def _approxaA(self, R, vR, vT, z, vz, phi, interp=True, cindx=None): def _approxaAInv(self, Or, Op, Oz, ar, ap, az, interp=True): """ - NAME: - _approxaAInv - PURPOSE: - return R,vR,... coordinates for a point based on the linear - approximation around the stream track - INPUT: - Or,Op,Oz,ar,ap,az - phase space coordinates in frequency-angle - space - interp= (True), if True, use the interpolated track - OUTPUT: - (R,vR,vT,z,vz,phi) - HISTORY: - 2013-12-22 - Written - Bovy (IAS) + Return R,vR,... coordinates for a point based on the linear approximation around the stream track + + Parameters + ---------- + Or : float + Radial action + Op : float + Azimuthal action + Oz : float + Vertical action + ar : float + Radial angle + ap : float + Azimuthal angle + az : float + Vertical angle + interp : bool, optional + If True, use the interpolated track. Default is True. + + Returns + ------- + tuple + (R,vR,vT,z,vz,phi) + + Notes + ----- + - 2013-12-22 - Written - Bovy (IAS) """ if isinstance(Or, (int, float, numpy.float32, numpy.float64)): # Scalar input Or = numpy.array([Or]) @@ -3009,50 +2960,34 @@ def _approxaAInv(self, Or, Op, Oz, ar, ap, az, interp=True): ################################EVALUATE THE DF################################ def __call__(self, *args, **kwargs): """ - NAME: - - __call__ - - PURPOSE: - - evaluate the DF - - INPUT: - - Either: - - a) R,vR,vT,z,vz,phi ndarray [nobjects] - - b) (Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) tuple if aAInput - - where: - - Omegar - radial frequency - - Omegaphi - azimuthal frequency - - Omegaz - vertical frequency - - angler - radial angle - - anglephi - azimuthal angle - - anglez - vertical angle - - c) Orbit instance or list thereof - - log= if True, return the natural log - - aaInput= (False) if True, option b above - - OUTPUT: - - value of DF - - HISTORY: - - 2013-12-03 - Written - Bovy (IAS) - + Evaluate the distribution function (DF). + + Parameters + ---------- + *args : tuple + Either: + a) R,vR,vT,z,vz,phi ndarray [nobjects] + b) (Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) tuple if aAInput where: + * Omegar - radial frequency + * Omegaphi - azimuthal frequency + * Omegaz - vertical frequency + * angler - radial angle + * anglephi - azimuthal angle + * anglez - vertical angle + c) Orbit instance or list thereof + log : bool, optional + If True, return the natural log. Default is True. + aaInput : bool, optional + If True, option b above. Default is False. + + Returns + ------- + ndarray + Value of DF. + + Notes + ----- + - 2013-12-03 - Written - Bovy (IAS) """ # First parse log log = kwargs.pop("log", True) @@ -3097,16 +3032,36 @@ def __call__(self, *args, **kwargs): def prepData4Call(self, *args, **kwargs): """ - NAME: - prepData4Call - PURPOSE: - prepare stream data for the __call__ method - INPUT: - __call__ inputs - OUTPUT: - (dOmega,dangle); wrt the progenitor; each [3,nobj] - HISTORY: - 2013-12-04 - Written - Bovy (IAS) + Prepare stream data for the __call__ method. + + Parameters + ---------- + *args : tuple + Either: + a) R,vR,vT,z,vz,phi ndarray [nobjects] + b) (Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) tuple if aAInput + where: + * Omegar - radial frequency + * Omegaphi - azimuthal frequency + * Omegaz - vertical frequency + * angler - radial angle + * anglephi - azimuthal angle + * anglez - vertical angle + c) Orbit instance or list thereof + log : bool, optional + If True, return the natural log. Default is True. + aaInput : bool, optional + If True, option b above. Default is False. + + Returns + ------- + tuple + Tuple containing two arrays (dOmega, dangle), each of shape [3, nobj], wrt the progenitor. + + Notes + ----- + - 2013-12-04 - Written - Bovy (IAS). + """ # First calculate the actionAngle coordinates if they're not given # as such @@ -3175,45 +3130,41 @@ def _parse_call_args(self, *args, **kwargs): def callMarg(self, xy, **kwargs): """ - NAME: - - callMarg - - PURPOSE: - evaluate the DF, marginalizing over some directions, in Galactocentric rectangular coordinates (or in observed l,b,D,vlos,pmll,pmbb) coordinates) - - INPUT: - - xy - phase-space point [X,Y,Z,vX,vY,vZ]; the distribution of the dimensions set to None is returned - - interp= (object-wide interp default) if True, use the interpolated stream track - - cindx= index of the closest point on the (interpolated) stream track if not given, determined from the dimensions given - - nsigma= (3) number of sigma to marginalize the DF over (approximate sigma) - - ngl= (5) order of Gauss-Legendre integration - - lb= (False) if True, xy contains [l,b,D,vlos,pmll,pmbb] in [deg,deg,kpc,km/s,mas/yr,mas/yr] and the marginalized PDF in these coordinates is returned - - vo= (220) circular velocity to normalize with when lb=True - - ro= (8) Galactocentric radius to normalize with when lb=True - - R0= (8) Galactocentric radius of the Sun (kpc) - - Zsun= (0.0208) Sun's height above the plane (kpc) - - vsun= ([-11.1,241.92,7.25]) Sun's motion in cylindrical coordinates (vR positive away from center) - - OUTPUT: - - p(xy) marginalized over missing directions in xy - - HISTORY: - - 2013-12-16 - Written - Bovy (IAS) - + Evaluate the distribution function (DF), marginalizing over some directions, in Galactocentric rectangular coordinates (or in observed l,b,D,vlos,pmll,pmbb) coordinates) + + Parameters + ---------- + xy : array_like + Phase-space point [X,Y,Z,vX,vY,vZ]; the distribution of the dimensions set to None is returned. + interp : bool, optional + If True, use the interpolated stream track. Default is True. + cindx : int, optional + Index of the closest point on the (interpolated) stream track if not given, determined from the dimensions given. + nsigma : int, optional + Number of sigma to marginalize the DF over (approximate sigma). Default is 3. + ngl : int, optional + Order of Gauss-Legendre integration. Default is 5. + lb : bool, optional + If True, xy contains [l,b,D,vlos,pmll,pmbb] in [deg,deg,kpc,km/s,mas/yr,mas/yr] and the marginalized PDF in these coordinates is returned. Default is False. + ro : float or Quantity, optional + Distance scale for translation into internal units (object-wide default). + vo : float or Quantity, optional + Velocity scale for translation into internal units (object-wide default). + R0 : float or Quantity, optional + Distance to the Galactic center (object-wide default). + Zsun : float or Quantity, optional + Sun's height above the plane (object-wide default). + vsun : float or Quantity, optional + Sun's motion in cylindrical coordinates (object-wide default). + + Returns + ------- + ndarray + Value of DF. + + Notes + ----- + - 2013-12-16 - Written - Bovy (IAS) """ coordGiven = numpy.array([not x is None for x in xy], dtype="bool") if numpy.sum(coordGiven) == 6: @@ -3361,32 +3312,27 @@ def callMarg(self, xy, **kwargs): def gaussApprox(self, xy, **kwargs): """ - NAME: - - gaussApprox - - PURPOSE: - - return the mean and variance of a Gaussian approximation to the stream DF at a given phase-space point in Galactocentric rectangular coordinates (distribution is over missing directions) - - INPUT: - - xy - phase-space point [X,Y,Z,vX,vY,vZ]; the distribution of the dimensions set to None is returned - - interp= (object-wide interp default) if True, use the interpolated stream track - - cindx= index of the closest point on the (interpolated) stream track if not given, determined from the dimensions given - - lb= (False) if True, xy contains [l,b,D,vlos,pmll,pmbb] in [deg,deg,kpc,km/s,mas/yr,mas/yr] and the Gaussian approximation in these coordinates is returned - - OUTPUT: - - (mean,variance) of the approximate Gaussian DF for the missing directions in xy - - HISTORY: - - 2013-12-12 - Written - Bovy (IAS) - + Return the mean and variance of a Gaussian approximation to the stream DF at a given phase-space point in Galactocentric rectangular coordinates (distribution is over missing directions) + + Parameters + ---------- + xy : array_like + Phase-space point [X,Y,Z,vX,vY,vZ]; the distribution of the dimensions set to None is returned. + interp : bool, optional + If True, use the interpolated stream track. Default is True. + cindx : int, optional + Index of the closest point on the (interpolated) stream track if not given, determined from the dimensions given. + lb : bool, optional + If True, xy contains [l,b,D,vlos,pmll,pmbb] in [deg,deg,kpc,km/s,mas/yr,mas/yr] and the Gaussian approximation in these coordinates is returned. Default is False. + + Returns + ------- + tuple + (mean,variance) of the approximate Gaussian DF for the missing directions in xy. + + Notes + ----- + - 2013-12-12 - Written - Bovy (IAS). """ interp = kwargs.get("interp", self._useInterp) lb = kwargs.get("lb", False) @@ -3456,36 +3402,31 @@ def sample( self, n, returnaAdt=False, returndt=False, interp=None, xy=False, lb=False ): """ - NAME: - - sample - - PURPOSE: - - sample from the DF - - INPUT: - - n - number of points to return - - returnaAdt= (False) if True, return (Omega,angle,dt) - - returndT= (False) if True, also return the time since the star was stripped - - interp= (object-wide default) use interpolation of the stream track - - xy= (False) if True, return Galactocentric rectangular coordinates - - lb= (False) if True, return Galactic l,b,d,vlos,pmll,pmbb coordinates - - OUTPUT: - - (R,vR,vT,z,vz,phi) of points on the stream in 6,N array - - HISTORY: - - 2013-12-22 - Written - Bovy (IAS) - + Sample from the DF. + + Parameters + ---------- + n : int + Number of points to return. + returnaAdt : bool, optional + If True, return (Omega,angle,dt). Default is False. + returndt : bool, optional + If True, also return the time since the star was stripped. Default is False. + interp : object-wide default, optional + Use interpolation of the stream track. Default is object-wide default. + xy : bool, optional + If True, return Galactocentric rectangular coordinates. Default is False. + lb : bool, optional + If True, return Galactic l,b,d,vlos,pmll,pmbb coordinates. Default is False. + + Returns + ------- + ndarray or tuple + (R,vR,vT,z,vz,phi) of points on the stream in 6,N array. Or (X,Y,Z,vX,vY,vZ) if xy is True. Or (l,b,d,vlos,pmll,pmbb) if lb is True. Or (Omega,angle,dt) if returnaAdt is True. If returndt is set, also return the time since the star was stripped. + + Notes + ----- + - 2013-12-22 - Written - Bovy (IAS) """ # First sample frequencies Om, angle, dt = self._sample_aAt(n) @@ -3670,16 +3611,21 @@ def _sample_aAt(self, n): def sample_t(self, n): """ - NAME: - sample_t - PURPOSE: - generate a stripping time (time since stripping); simple implementation could be replaced by more complicated distributions in sub-classes of streamdf - INPUT: - n - number of points to return - OUTPUT: - array of n stripping times - HISTORY: - 2015-09-16 - Written - Bovy (UofT) + Sample the time since the progenitor was stripped + + Parameters + ---------- + n : int + Number of points to return + + Returns + ------- + ndarray + Array of n stripping times + + Notes + ----- + - 2015-09-16 - Written - Bovy (UofT) """ return numpy.random.uniform(size=n) * self._tdisrupt @@ -3921,39 +3867,49 @@ def calcaAJac( _initacfs=None, ): """ - NAME: - calcaAJac - PURPOSE: - calculate the Jacobian d(J,theta)/d(x,v) - INPUT: - xv - phase-space point: Either - 1) [R,vR,vT,z,vz,phi] - 2) [l,b,D,vlos,pmll,pmbb] (if lb=True, see below) - 3) list/array of 6 numbers that can be transformed into (normalized) R,vR,vT,z,vz,phi using coordFunc - - aA - actionAngle instance - - dxv - infinitesimal to use (rescaled for lb, so think fractionally)) - - freqs= (False) if True, go to frequencies rather than actions - - dOdJ= (False), actually calculate d Frequency / d action - - actionsFreqsAngles= (False) if True, calculate d(action,freq.,angle)/d (xv) - - lb= (False) if True, start with (l,b,D,vlos,pmll,pmbb) in (deg,deg,kpc,km/s,mas/yr,mas/yr) - vo= (220) circular velocity to normalize with when lb=True - ro= (8) Galactocentric radius to normalize with when lb=True - R0= (8) Galactocentric radius of the Sun (kpc) - Zsun= (0.0208) Sun's height above the plane (kpc) - vsun= ([-11.1,241.92,7.25]) Sun's motion in cylindrical coordinates (vR positive away from center) - - coordFunc= (None) if set, this is a function that takes xv and returns R,vR,vT,z,vz,phi in normalized units (units where vc=1 at r=1 if the potential is normalized that way, for example) - - OUTPUT: - Jacobian matrix - HISTORY: - 2013-11-25 - Written - Bovy (IAS) + Calculate the Jacobian d(J,theta)/d(x,v) + + Parameters + ---------- + xv: numpy.ndarray + Either: + 1) [R,vR,vT,z,vz,phi] + 2) [l,b,D,vlos,pmll,pmbb] (if lb=True, see below) + 3) list/array of 6 numbers that can be transformed into (normalized) R,vR,vT,z,vz,phi using coordFunc + aA: actionAngle instance + dxv: float, optional + Infinitesimal to use (rescaled for lb, so think fractionally)) + freqs: bool, optional + If True, go to frequencies rather than actions + dOdJ: bool, optional + Actually calculate d Frequency / d action + actionsFreqsAngles: bool, optional + If True, calculate d(action,freq.,angle)/d (xv) + lb: bool, optional + If True, start with (l,b,D,vlos,pmll,pmbb) in (deg,deg,kpc,km/s,mas/yr,mas/yr) + vo: float, optional + Circular velocity to normalize with when lb=True + ro: float, optional + Galactocentric radius to normalize with when lb=True + R0: float, optional + Galactocentric radius of the Sun (kpc) + Zsun: float, optional + Sun's height above the plane (kpc) + vsun: list, optional + Sun's motion in cylindrical coordinates (vR positive away from center) + coordFunc: function, optional + If set, this is a function that takes xv and returns R,vR,vT,z,vz,phi in normalized units (units where vc=1 at r=1 if the potential is normalized that way, for example) + _initacfs: list, optional + If set, this is a list of the actions, frequencies, and angles at xv + + Returns + ------- + jac: numpy.ndarray + Jacobian matrix + + Notes + ----- + - 2013-11-25 - Written - Bovy (IAS) """ if lb: coordFunc = lambda x: lbCoordFunc(xv, vo, ro, R0, Zsun, vsun) diff --git a/galpy/df/streamgapdf.py b/galpy/df/streamgapdf.py index 21133a18b..13ab58f97 100644 --- a/galpy/df/streamgapdf.py +++ b/galpy/df/streamgapdf.py @@ -40,58 +40,105 @@ class streamgapdf(streamdf.streamdf): def __init__(self, *args, **kwargs): """ - NAME: - - __init__ - - PURPOSE: - - Initialize the DF of a gap in a stellar stream - - INPUT: - - streamdf args and kwargs - - Subhalo and impact parameters: - - impactb= impact parameter (can be Quantity) - - subhalovel= velocity of the subhalo shape=(3) (can be Quantity) - - timpact time since impact (can be Quantity) - - impact_angle= angle offset from progenitor at which the impact occurred (rad) (can be Quantity) - - Subhalo: specify either 1( mass and size of Plummer sphere or 2( general spherical-potential object (kick is numerically computed) - - 1( GM= mass of the subhalo (can be Quantity) - - rs= size parameter of the subhalo (can be Quantity) - - 2( subhalopot= galpy potential object or list thereof (should be spherical) - - 3( hernquist= (False) if True, use Hernquist kicks for GM/rs - - deltaAngleTrackImpact= (None) angle to estimate the stream track over to determine the effect of the impact [similar to deltaAngleTrack] (rad) - - nTrackChunksImpact= (floor(deltaAngleTrack/0.15)+1) number of chunks to divide the progenitor track in near the impact [similar to nTrackChunks] - - nKickPoints= (30xnTrackChunksImpact) number of points along the stream to compute the kicks at (kicks are then interpolated); '30' chosen such that higherorderTrack can be set to False and get calculations accurate to > 99% - - nokicksetup= (False) if True, only run as far as setting up the coordinate transformation at the time of impact (useful when using this in streampepperdf) - - spline_order= (3) order of the spline to interpolate the kicks with - - higherorderTrack= (False) if True, calculate the track using higher-order terms - - OUTPUT: - - object - - HISTORY: - - 2015-06-02 - Started - Bovy (IAS) - + Initialize the DF of a gap in a stellar stream + + Parameters + ---------- + sigv : float or Quantity + Radial velocity dispersion of the progenitor. + progenitor : galpy.orbit.Orbit + Progenitor orbit as Orbit instance (will be re-integrated, so don't bother integrating the orbit before). + pot : galpy.potential.Potential or list thereof, optional + Potential instance or list thereof. + aA : actionAngle instance + ActionAngle instance used to convert (x,v) to actions. Generally a actionAngleIsochroneApprox instance. + useTM : bool, optional + If set to an actionAngleTorus instance, use this to speed up calculations. + tdisrupt : float or Quantity, optional + Time since start of disruption (default: 5 Gyr). + sigMeanOffset : float, optional + Offset between the mean of the frequencies and the progenitor, in units of the largest eigenvalue of the frequency covariance matrix (along the largest eigenvector), should be positive; to model the trailing part, set leading=False (default: 6.0). + leading : bool, optional + If True, model the leading part of the stream; if False, model the trailing part (default: True). + sigangle : float or Quantity, optional + Estimate of the angle spread of the debris initially (default: sigv/122/[1km/s]=1.8sigv in natural coordinates). + deltaAngleTrack : float or Quantity, optional + Angle to estimate the stream track over (rad; or can be Quantity) (default: None). + nTrackChunks : int, optional + Number of chunks to divide the progenitor track in (default: floor(deltaAngleTrack/0.15)+1). + nTrackIterations : int, optional + Number of iterations to perform when establishing the track; each iteration starts from a previous approximation to the track in (x,v) and calculates a new track based on the deviation between the previous track and the desired track in action-angle coordinates; if not set, an appropriate value is determined based on the magnitude of the misalignment between stream and orbit, with larger numbers of iterations for larger misalignments (default: None). + progIsTrack : bool, optional + If True, then the progenitor (x,v) is actually the (x,v) of the stream track at zero angle separation; useful when initializing with an orbit fit; the progenitor's position will be calculated (default: False). + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + Vnorm : float or Quantity, optional + Deprecated. Use vo instead (default: None). + Rnorm : float or Quantity, optional + Deprecated. Use ro instead (default: None). + R0 : float or Quantity, optional + Galactocentric radius of the Sun (kpc) (can be different from ro) (default: 8.0). + Zsun : float or Quantity, optional + Sun's height above the plane (kpc) (default: 0.0208). + vsun : array_like or Quantity, optional + Sun's motion in cylindrical coordinates (vR positive away from center) (can be Quantity array, but not a list of Quantities) (default: [-11.1, 8.0 * 30.24, 7.25]). + multi : int, optional + If set, use multi-processing (default: None). + interpTrack : bool, optional + Interpolate the stream track while setting up the instance (can be done by hand by calling self._interpolate_stream_track() and self._interpolate_stream_track_aA()) (default: _INTERPDURINGSETUP). + useInterp : bool, optional + Use interpolation by default when calculating approximated frequencies and angles (default: _USEINTERP). + nosetup : bool, optional + If True, don't setup the stream track and anything else that is expensive (default: False). + nospreadsetup : bool, optional + If True, don't setup the spread around the stream track (only for nosetup is False) (default: False). + approxConstTrackFreq : bool, optional + If True, approximate the stream assuming that the frequency is constant along the stream (only works with useTM, for which this leads to a significant speed-up) (default: False). + useTMHessian : bool, optional + If True, compute the basic Hessian dO/dJ_prog using TM; otherwise use aA (default: False). + custom_transform : array_like, optional + Matrix implementing the rotation from (ra,dec) to a custom set of sky coordinates (default: None). + impactb : float or Quantity, optional + Impact parameter (can be Quantity) (default: 1.0). + subhalovel : numpy.ndarray or Quantity, optional + Velocity of the subhalo shape=(3) (default: [0.0, 1.0, 0.0]). + timpact : float or Quantity, optional + Time since impact (can be Quantity) (default: 1.0). + impact_angle : float or Quantity, optional + Angle offset from progenitor at which the impact occurred (rad) (can be Quantity) (default: 1.0). + GM : float or Quantity, optional + Mass of the subhalo when using a Plummer or Hernquist model. + rs : float or Quantity, optional + Scale parameter of the subhalo when using a Plummer or Hernquist model. + hernquist : bool, optional + If True, use Hernquist kicks for GM/rs (default: False --> Plummer). + subhalopot : Potential or list thereof, optional + Gravitational potential of the subhalo (alternative to specifying GM and rs) + deltaAngleTrackImpact : float or Quantity, optional + Angle to estimate the stream track over to determine the effect of the impact [similar to deltaAngleTrack] (rad) (default: None). + nTrackChunksImpact : int, optional + Number of chunks to divide the progenitor track in near the impact [similar to nTrackChunks] (default: floor(deltaAngleTrack/0.15)+1). + nKickPoints : int, optional + Number of points along the stream to compute the kicks at (kicks are then interpolated); '30' chosen such that higherorderTrack can be set to False and get calculations accurate to > 99% (default: 30xnTrackChunksImpact). + nokicksetup : bool, optional + If True, only run as far as setting up the coordinate transformation at the time of impact (useful when using this in streampepperdf) (default: False). + spline_order : int, optional + Order of the spline to interpolate the kicks with (default: 3). + higherorderTrack : bool, optional + If True, calculate the track using higher-order terms (default: False). + nTrackChunks : int, optional + Number of chunks to divide the progenitor track into (default: 8). + interpTrack : bool, optional + If True, interpolate the track (default: True). + useInterp : bool, optional + If True, use the interpolated track to calculate actions and angles (default: True). + + Notes + ----- + - Parameters above up to impactb are streamdf parameters used to setup the underlying smooth stream. + - 2015-06-02 - Started - Bovy (IAS) """ df.__init__(self, ro=kwargs.get("ro", None), vo=kwargs.get("vo", None)) # Parse kwargs @@ -171,27 +218,23 @@ def __init__(self, *args, **kwargs): def pOparapar(self, Opar, apar): """ - NAME: - - pOparapar - - PURPOSE: - - return the probability of a given parallel (frequency,angle) offset pair + Return the probability of a given parallel (frequency,angle) offset pair. - INPUT: + Parameters + ---------- + Opar : numpy.ndarray or Quantity + Parallel frequency offset. + apar : float or Quantity + Parallel angle offset along the stream. - Opar - parallel frequency offset (array) (can be Quantity) + Returns + ------- + numpy.ndarray + Probability of a given parallel (frequency,angle) offset pair. - apar - parallel angle offset along the stream (scalar) (can be Quantity) - - OUTPUT: - - p(Opar,apar) - - HISTORY: - - 2015-11-17 - Written - Bovy (UofT) + Notes + ----- + - 2015-11-17 - Written - Bovy (UofT). """ Opar = conversion.parse_frequency(Opar, ro=self._ro, vo=self._vo) @@ -360,25 +403,25 @@ def _densMoments_approx_higherorder_gaussxpolyInts(self, ll, ul, maxj): def minOpar(self, dangle, tdisrupt=None, _return_raw=False): """ - NAME: - - minOpar - - PURPOSE: - - return the approximate minimum parallel frequency at a given angle - - INPUT: - - dangle - parallel angle - - OUTPUT: - - minimum frequency that gets to this parallel angle - - HISTORY: - - 2015-12-28 - Written - Bovy (UofT) + Return the approximate minimum parallel frequency at a given angle + + Parameters + ---------- + dangle : float + Parallel angle + tdisrupt : float, optional + Disruption time (default is the value passed at initialization) + _return_raw : bool, optional + If True, return the index of the minimum frequency and the value of the minimum frequency (default is False) + + Returns + ------- + float or tuple + Minimum frequency that gets to this parallel angle or a tuple with the index of the minimum frequency and the value of the minimum frequency + + Notes + ----- + - 2015-12-28 - Written - Bovy (UofT) """ if tdisrupt is None: @@ -407,31 +450,29 @@ def meanOmega( self, dangle, oned=False, tdisrupt=None, approx=True, higherorder=None ): """ - NAME: - - meanOmega - - PURPOSE: - - calculate the mean frequency as a function of angle, assuming a uniform time distribution up to a maximum time - - INPUT: - - dangle - angle offset - - oned= (False) if True, return the 1D offset from the progenitor (along the direction of disruption) - - approx= (True) if True, compute the mean Omega by direct integration of the spline representation - - higherorder= (object-wide default higherorderTrack) if True, include higher-order spline terms in the approximate computation - - OUTPUT: - - mean Omega - - HISTORY: - - 2015-11-17 - Written - Bovy (UofT) + Calculate the mean frequency as a function of angle, assuming a uniform time distribution up to a maximum time. + + Parameters + ---------- + dangle : float + Angle offset. + oned : bool, optional + If True, return the 1D offset from the progenitor (along the direction of disruption). Default is False. + tdisrupt : float, optional + Maximum time. Default is None. + approx : bool, optional + If True, compute the mean Omega by direct integration of the spline representation. Default is True. + higherorder : object, optional + Higher-order spline terms in the approximate computation. Default is object-wide default higherorderTrack. + + Returns + ------- + float + Mean Omega. + + Notes + ----- + - 2015-11-17 - Written - Bovy (UofT) """ if higherorder is None: @@ -1233,36 +1274,31 @@ def _sample_aAt(self, n): def impulse_deltav_plummer(v, y, b, w, GM, rs): """ - NAME: - - impulse_deltav_plummer - - PURPOSE: - - calculate the delta velocity to due an encounter with a Plummer sphere in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - y - position along the stream (nstar) - - b - impact parameter - - w - velocity of the Plummer sphere (3) - - GM - mass of the Plummer sphere (in natural units) - - rs - size of the Plummer sphere - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-04-30 - Written based on Erkal's expressions - Bovy (IAS) - + Calculate the delta velocity to due an encounter with a Plummer sphere in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + y : numpy.ndarray + position along the stream (nstar) + b : float + impact parameter + w : numpy.ndarray + velocity of the Plummer sphere (3) + GM : float + mass of the Plummer sphere (in natural units) + rs : float + size of the Plummer sphere + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-04-30 - Written based on Erkal's expressions - Bovy (IAS) """ if len(v.shape) == 1: v = numpy.reshape(v, (1, 3)) @@ -1307,40 +1343,35 @@ def impulse_deltav_plummer(v, y, b, w, GM, rs): def impulse_deltav_plummer_curvedstream(v, x, b, w, x0, v0, GM, rs): """ - NAME: - - impulse_deltav_plummer_curvedstream - - PURPOSE: - - calculate the delta velocity to due an encounter with a Plummer sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - b - impact parameter - - w - velocity of the Plummer sphere (3) - - x0 - point of closest approach - - v0 - velocity of point of closest approach - - GM - mass of the Plummer sphere (in natural units) - - rs - size of the Plummer sphere - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-05-04 - Written based on above - SANDERS - + Calculate the delta velocity to due an encounter with a Plummer sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + b : float + impact parameter + w : numpy.ndarray + velocity of the Plummer sphere (3) + x0 : numpy.ndarray + point of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + GM : float + mass of the Plummer sphere (in natural units) + rs : float + size of the Plummer sphere + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-05-04 - Written based on above - Sanders (Cambridge) """ if len(v.shape) == 1: v = numpy.reshape(v, (1, 3)) @@ -1373,35 +1404,31 @@ def HernquistX(s): def impulse_deltav_hernquist(v, y, b, w, GM, rs): """ - NAME: - - impulse_deltav_hernquist - - PURPOSE: - - calculate the delta velocity to due an encounter with a Hernquist sphere in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - y - position along the stream (nstar) - - b - impact parameter - - w - velocity of the Hernquist sphere (3) - - GM - mass of the Hernquist sphere (in natural units) - - rs - size of the Hernquist sphere - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-08-13 SANDERS, using Wyn Evans calculation + Calculate the delta velocity to due an encounter with a Hernquist sphere in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + y : numpy.ndarray + position along the stream (nstar) + b : float + impact parameter + w : numpy.ndarray + velocity of the Hernquist sphere (3) + GM : float + mass of the Hernquist sphere (in natural units) + rs : float + size of the Hernquist sphere + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-08-13 - Written using Wyn Evans calculation - Sanders (Cambridge) """ if len(v.shape) == 1: @@ -1466,39 +1493,35 @@ def impulse_deltav_hernquist(v, y, b, w, GM, rs): def impulse_deltav_hernquist_curvedstream(v, x, b, w, x0, v0, GM, rs): """ - NAME: - - impulse_deltav_plummer_hernquist - - PURPOSE: - - calculate the delta velocity to due an encounter with a Hernquist sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - b - impact parameter - - w - velocity of the Hernquist sphere (3) - - x0 - point of closest approach - - v0 - velocity of point of closest approach - - GM - mass of the Hernquist sphere (in natural units) - - rs - size of the Hernquist sphere - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-08-13 - SANDERS, using Wyn Evans calculation + Calculate the delta velocity to due an encounter with a Hernquist sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + b : float + impact parameter + w : numpy.ndarray + velocity of the Hernquist sphere (3) + x0 : numpy.ndarray + point of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + GM : float + mass of the Hernquist sphere (in natural units) + rs : float + size of the Hernquist sphere + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-08-13 - Written using Wyn Evans calculation - Sanders (Cambridge) """ if len(v.shape) == 1: @@ -1538,37 +1561,30 @@ def _deltav_integrate(y, b, w, pot): def impulse_deltav_general(v, y, b, w, pot): """ - NAME: - - impulse_deltav_general - - - PURPOSE: - - calculate the delta velocity to due an encounter with a general spherical potential in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - y - position along the stream (nstar) - - b - impact parameter - - w - velocity of the subhalo (3) - - pot - Potential object or list thereof (should be spherical) - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-05-04 - SANDERS - - 2015-06-15 - Tweak to use galpy' potential objects - Bovy (IAS) - + Calculate the delta velocity to due an encounter with a general spherical potential in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + y : numpy.ndarray + position along the stream (nstar) + b : float + impact parameter + w : numpy.ndarray + velocity of the subhalo (3) + pot : Potential object or list thereof + Potential object or list thereof (should be spherical) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-05-04 - Written - Sanders (Cambridge) + - 2015-06-15 - Tweak to use galpy' potential objects - Bovy (IAS) """ pot = flatten_potential(pot) if len(v.shape) == 1: @@ -1596,40 +1612,34 @@ def impulse_deltav_general(v, y, b, w, pot): def impulse_deltav_general_curvedstream(v, x, b, w, x0, v0, pot): """ - NAME: - - impulse_deltav_general_curvedstream - - PURPOSE: - - calculate the delta velocity to due an encounter with a general spherical potential in the impulse approximation; allows for arbitrary velocity vectors and arbitrary shaped streams - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - b - impact parameter - - w - velocity of the subhalo (3) - - x0 - position of closest approach (3) - - v0 - velocity of stream at closest approach (3) - - pot - Potential object or list thereof (should be spherical) - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-05-04 - SANDERS - - 2015-06-15 - Tweak to use galpy' potential objects - Bovy (IAS) - + Calculate the delta velocity to due an encounter with a general spherical potential in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + b : float + impact parameter + w : numpy.ndarray + velocity of the subhalo (3) + x0 : numpy.ndarray + point of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + pot : Potential object or list thereof + Potential object or list thereof (should be spherical) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-05-04 - Written - Sanders (Cambridge) + - 2015-06-15 - Tweak to use galpy' potential objects - Bovy (IAS) """ pot = flatten_potential(pot) if len(v.shape) == 1: @@ -1659,46 +1669,43 @@ def impulse_deltav_general_orbitintegration( integrate_method="symplec4_c", ): """ - NAME: - - impulse_deltav_general_orbitintegration - - PURPOSE: - - calculate the delta velocity to due an encounter with a general spherical potential NOT in the impulse approximation by integrating each particle in the underlying galactic potential; allows for arbitrary velocity vectors and arbitrary shaped streams. - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - b - impact parameter - - w - velocity of the subhalo (3) - - x0 - position of closest approach (3) - - v0 - velocity of stream at closest approach (3) - - pot - Potential object or list thereof (should be spherical) - - tmax - maximum integration time - - galpot - galpy Potential object or list thereof - - nsamp(1000) - number of forward integration points - - integrate_method= ('symplec4_c') orbit integrator to use (see Orbit.integrate) - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-08-17 - SANDERS - + Calculate the delta velocity to due an encounter with a general spherical potential NOT in the impulse approximation by integrating each particle in the underlying galactic potential; allows for arbitrary velocity vectors and arbitrary shaped streams. + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + b : float + impact parameter + w : numpy.ndarray + velocity of the subhalo (3) + x0 : numpy.ndarray + position of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + pot : Potential object or list thereof + Potential object or list thereof (should be spherical) + tmax : float + maximum integration time + galpot : Potential object or list thereof + galpy Potential object or list thereof + tmaxfac : float + multiple of rs/fabs(w - v0) to use for time integration interval + nsamp : int + number of forward integration points + integrate_method : str + orbit integrator to use (see Orbit.integrate) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-08-17 - Written - Sanders (Cambridge) """ galpot = flatten_potential(galpot) if len(v.shape) == 1: @@ -1748,48 +1755,43 @@ def impulse_deltav_general_fullplummerintegration( integrate_method="symplec4_c", ): """ - NAME: - - impulse_deltav_general_fullplummerintegration - - PURPOSE: - - calculate the delta velocity to due an encounter with a moving Plummer sphere and galactic potential relative to just in galactic potential - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - b - impact parameter - - w - velocity of the subhalo (3) - - x0 - position of closest approach (3) - - v0 - velocity of stream at closest approach (3) - - galpot - Galaxy Potential object - - GM - mass of Plummer - - rs - scale of Plummer - - tmaxfac(10) - multiple of rs/fabs(w - v0) to use for time integration interval - - N(1000) - number of forward integration points - - integrate_method('symplec4_c') - orbit integrator to use (see Orbit.integrate) - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-08-18 - SANDERS - + Calculate the delta velocity to due an encounter with a moving Plummer sphere and galactic potential relative to just in galactic potential by integrating each particle in the underlying galactic potential; allows for arbitrary velocity vectors and arbitrary shaped streams. + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + b : float + impact parameter + w : numpy.ndarray + velocity of the subhalo (3) + x0 : numpy.ndarray + position of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + galpot : Potential object or list thereof + galpy Potential object or list thereof + GM : float + mass of Plummer + rs : float + scale of Plummer + tmaxfac : float + multiple of rs/fabs(w - v0) to use for time integration interval + N : int + number of forward integration points + integrate_method : str + orbit integrator to use (see Orbit.integrate) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-08-18 - Written - Sanders (Cambridge) """ galpot = flatten_potential(galpot) if len(v.shape) == 1: @@ -1854,38 +1856,35 @@ def _astream_integrand_z(t, y, v, b, w, b2, w2, wperp, wperp2, wpar, GSigma, rs2 def impulse_deltav_plummerstream(v, y, b, w, GSigma, rs, tmin=None, tmax=None): """ - NAME: - - impulse_deltav_plummerstream - - PURPOSE: - - calculate the delta velocity to due an encounter with a Plummer-softened stream in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream - - INPUT: - - v - velocity of the stream (nstar,3) - - y - position along the stream (nstar) - - b - impact parameter - - w - velocity of the Plummer sphere (3) - - GSigma - surface density of the Plummer-softened stream (in natural units); should be a function of time - - rs - size of the Plummer sphere - - tmin, tmax= (None) minimum and maximum time to consider for GSigma (need to be set) - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-11-14 - Written - Bovy (UofT) - + Calculate the delta velocity to due an encounter with a Plummer-softened stream in the impulse approximation; allows for arbitrary velocity vectors, but y is input as the position along the stream + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + y : numpy.ndarray + position along the stream (nstar) + b : float + impact parameter + w : numpy.ndarray + velocity of the Plummer sphere (3) + GSigma : function + surface density of the Plummer-softened stream (in natural units); should be a function of time + rs : float + size of the Plummer sphere + tmin : float + minimum time to consider for GSigma (need to be set) + tmax : float + maximum time to consider for GSigma (need to be set) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-11-14 - Written - Bovy (UofT) """ if len(v.shape) == 1: v = numpy.reshape(v, (1, 3)) @@ -2000,46 +1999,43 @@ def impulse_deltav_plummerstream_curvedstream( v, x, t, b, w, x0, v0, GSigma, rs, galpot, tmin=None, tmax=None ): """ - NAME: - - impulse_deltav_plummerstream_curvedstream - - PURPOSE: - - calculate the delta velocity to due an encounter with a Plummer sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream; velocities and positions are assumed to lie along an orbit - - INPUT: - - v - velocity of the stream (nstar,3) - - x - position along the stream (nstar,3) - - t - times at which (v,x) are reached, wrt the closest impact t=0 (nstar) - - b - impact parameter - - w - velocity of the Plummer sphere (3) - - x0 - point of closest approach - - v0 - velocity of point of closest approach - - GSigma - surface density of the Plummer-softened stream (in natural units); should be a function of time - - rs - size of the Plummer sphere - - galpot - galpy Potential object or list thereof - - tmin, tmax= (None) minimum and maximum time to consider for GSigma - - OUTPUT: - - deltav (nstar,3) - - HISTORY: - - 2015-11-20 - Written based on Plummer sphere above - Bovy (UofT) - + Calculate the delta velocity to due an encounter with a Plummer-softened stream in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream; velocities and positions are assumed to lie along an orbit + + Parameters + ---------- + v : numpy.ndarray + velocity of the stream (nstar,3) + x : numpy.ndarray + position along the stream (nstar,3) + t : numpy.ndarray + times at which (v,x) are reached, wrt the closest impact t=0 (nstar) + b : float + impact parameter + w : numpy.ndarray + velocity of the Plummer sphere (3) + x0 : numpy.ndarray + point of closest approach + v0 : numpy.ndarray + velocity of point of closest approach + GSigma : function + surface density of the Plummer-softened stream (in natural units); should be a function of time + rs : float + size of the Plummer sphere + galpot : Potential object or list thereof + galpy Potential object or list thereof + tmin : float + minimum time to consider for GSigma (need to be set) + tmax : float + maximum time to consider for GSigma (need to be set) + + Returns + ------- + numpy.ndarray + velocity kick deltav (nstar,3) + + Notes + ----- + - 2015-11-14 - Written - Bovy (UofT) """ galpot = flatten_potential(galpot) if len(v.shape) == 1: diff --git a/galpy/df/streamspraydf.py b/galpy/df/streamspraydf.py index 8626dac2b..fed4e88bd 100644 --- a/galpy/df/streamspraydf.py +++ b/galpy/df/streamspraydf.py @@ -29,46 +29,39 @@ def __init__( vo=None, ): """ - NAME: - - __init__ - PURPOSE: - - Initialize a stream spray DF model of a tidal stream - - INPUT: - - progenitor_mass - mass of the progenitor (can be Quantity) - - tdisrupt= (5 Gyr) time since start of disruption (can be Quantity) - - leading= (True) if True, model the leading part of the stream - if False, model the trailing part - - progenitor= progenitor orbit as Orbit instance (will be re-integrated, so don't bother integrating the orbit before) - - meankvec= (Fardal+2015-ish defaults) - - sigkvec= (Fardal+2015-ish defaults) - - pot = (None) potential for integrating orbits - - rtpot = (pot) potential for calculating tidal radius and circular velocity (should generally be the same as pot, but sometimes you need to drop parts of the potential that don't allow the tidal radius / circular velocity to be computed, such as velocity-dependent forces; when using center, rtpot should be the relevant potential in the frame of the center, thus, also being different from pot) - - center = (None) Orbit instance that represents the center around which the progenitor is orbiting for the purpose of stream formation; allows for a stream to be generated from a progenitor orbiting a moving object, like a satellite galaxy. Integrated internally using centerpot. - - centerpot = (pot) potential for calculating the orbit of the center; this might be different from the potential that the progenitor is integrated in if, for example, dynamical friction is important for the orbit of the center (if it's a satellite). - - OUTPUT: - - Instance - - HISTORY: - - 2018-07-31 - Written - Bovy (UofT) - - 2021-05-05 - Added center keyword - Yansong Qian (UofT) - + Initialize a stream spray DF model of a tidal stream + + Parameters + ---------- + progenitor_mass : float or Quantity + Mass of the progenitor. + progenitor : galpy.orbit.Orbit, optional + Progenitor orbit as Orbit instance (will be re-integrated, so don't bother integrating the orbit before). + pot : galpy.potential.Potential or list of such instances, optional + Potential for integrating orbits. + rtpot : galpy.potential.Potential or list of such instances, optional + Potential for calculating tidal radius and circular velocity (should generally be the same as pot, but sometimes you need to drop parts of the potential that don't allow the tidal radius / circular velocity to be computed, such as velocity-dependent forces; when using center, rtpot should be the relevant potential in the frame of the center, thus, also being different from pot). + tdisrupt : float or Quantity, optional + Time since start of disruption. Default is 5 Gyr. + leading : bool, optional + If True, model the leading part of the stream. If False, model the trailing part. Default is True. + center : galpy.orbit.Orbit, optional + Orbit instance that represents the center around which the progenitor is orbiting for the purpose of stream formation; allows for a stream to be generated from a progenitor orbiting a moving object, like a satellite galaxy. Integrated internally using centerpot. + centerpot : galpy.potential.Potential or list of such instances, optional + Potential for calculating the orbit of the center; this might be different from the potential that the progenitor is integrated in if, for example, dynamical friction is important for the orbit of the center (if it's a satellite). + meankvec : list or array, optional + Mean of the action-angle distribution. Default is [2.0, 0.0, 0.3, 0.0, 0.0, 0.0]. + sigkvec : list or array, optional + Dispersion of the action-angle distribution. Default is [0.4, 0.0, 0.4, 0.5, 0.5, 0.0]. + ro : float or Quantity, optional + Distance scale for translation into internal units (default from configuration file). + vo : float or Quantity, optional + Velocity scale for translation into internal units (default from configuration file). + + Notes + ----- + - 2018-07-31 - Written - Bovy (UofT) + - 2021-05-05 - Added center keyword - Yansong Qian (UofT) """ df.__init__(self, ro=ro, vo=vo) self._progenitor_mass = conversion.parse_mass( @@ -119,38 +112,28 @@ def __init__( def sample(self, n, return_orbit=True, returndt=False, integrate=True): """ - NAME: - - sample - - PURPOSE: - - sample from the DF - - INPUT: - - n - number of points to return - - return_orbit= (True) If True, the output phase-space positions is an orbit.Orbit object, if False, the output is (R,vR,vT,z,vz,phi) - - returndt= (False) if True, also return the time since the star was stripped - - integrate= (True) if True, integrate the orbits to the present time, if False, return positions at stripping (probably want to combine with returndt=True then to make sense of them!) - - xy= (False) if True, return Galactocentric rectangular coordinates - - lb= (False) if True, return Galactic l,b,d,vlos,pmll,pmbb coordinates - - OUTPUT: - - Orbit instance or (R,vR,vT,z,vz,phi) of points on the stream in 6,N array (set of 6 Quantities when physical output is on); optionally the time is included as well. The ro/vo unit-conversion parameters and the zo/solarmotion parameters as well as whether physical outputs are on match the settings of the progenitor Orbit given to the class initialization - - HISTORY: - - 2018-07-31 - Written - Bovy (UofT) - - 2022-05-18 - Made output Orbit ro/vo/zo/solarmotion/roSet/voSet match that of the progenitor orbit - Bovy (UofT) - + Sample from the DF + + Parameters + ---------- + n : int + Number of points to return. + return_orbit : bool, optional + If True, the output phase-space positions is an orbit.Orbit object. If False, the output is (R,vR,vT,z,vz,phi). Default is True. + returndt : bool, optional + If True, also return the time since the star was stripped. Default is False. + integrate : bool, optional + If True, integrate the orbits to the present time. If False, return positions at stripping (probably want to combine with returndt=True then to make sense of them!). Default is True. + + Returns + ------- + Orbit, numpy.ndarray, or tuple + Orbit instance or (R,vR,vT,z,vz,phi) of points on the stream in 6,N array (set of 6 Quantities when physical output is on); optionally the time is included as well. The ro/vo unit-conversion parameters and the zo/solarmotion parameters as well as whether physical outputs are on, match the settings of the progenitor Orbit given to the class initialization + + Notes + ----- + - 2018-07-31 - Written - Bovy (UofT) + - 2022-05-18 - Made output Orbit ro/vo/zo/solarmotion/roSet/voSet match that of the progenitor orbit - Bovy (UofT) """ # First sample times dt = numpy.random.uniform(size=n) * self._tdisrupt diff --git a/galpy/df/surfaceSigmaProfile.py b/galpy/df/surfaceSigmaProfile.py index 612923796..b8dc0529a 100644 --- a/galpy/df/surfaceSigmaProfile.py +++ b/galpy/df/surfaceSigmaProfile.py @@ -21,17 +21,18 @@ def __init__(self): def formatStringParams(self): """ - NAME: - formatStringParams - PURPOSE: - when writing the parameters of this profile, what - format-strings to use? - This function defaults to '%6.4f' for each parameter in self._params - INPUT: - OUTPUT: - tuple of format-strings - HISTORY: - 2010-03-28 - Written - Bovy (NYU) + Return the format-strings to use when writing the parameters of this profile. + + Returns + ------- + tuple of str + Tuple of format-strings to use for each parameter in self._params. + + Notes + ----- + This function defaults to '%6.4f' for each parameter in self._params. + + - 2010-03-28 - Written - Bovy (NYU) """ out = [] for param in self._params: @@ -40,31 +41,39 @@ def formatStringParams(self): def outputParams(self): """ - NAME: - outputParams - PURPOSE:' - return a list of parameters of this profile, to create filenames - INPUT: - OUTPUT: - tuple of parameters (given to % ...) - HISTORY: - 2010-03-28 - Written - Bovy + Return a tuple of parameters of this profile, to create filenames. + + Returns + ------- + tuple + Tuple of parameters (given to % ...). + + Notes + ----- + - 2010-03-28 - Written - Bovy """ return tuple(self._params) def surfacemass(self, R, log=False): """ - NAME: - surfacemass - PURPOSE: - return the surface density profile at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the log (default: False) - OUTPUT: - Sigma(R) - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Return the surface density profile at this R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + float + Sigma(R) + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU) + """ raise NotImplementedError( "'surfacemass' function not implemented for this surfaceSigmaProfile class" @@ -72,17 +81,23 @@ def surfacemass(self, R, log=False): def sigma2(self, R, log=False): """ - NAME: - sigma2 - PURPOSE: - return the radial velocity variance at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the log (default: False) - OUTPUT: - sigma^2(R) - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Return the radial velocity variance at this R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + float + sigma^2(R). + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU) """ raise NotImplementedError( "'sigma2' function not implemented for this surfaceSigmaProfile class" @@ -94,33 +109,40 @@ class expSurfaceSigmaProfile(surfaceSigmaProfile): def __init__(self, params=(1.0 / 3.0, 1.0, 0.2)): """ - NAME: - __init__ - PURPOSE: - initialize an exponential surface density and sigma_R^2 profile - INPUT: - params - tuple/list of (surface scale-length,sigma scale-length, - sigma(ro)) (NOTE: *not* sigma2 scale-length) - OUTPUT: - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Initialize an exponential surface density and sigma_R^2 profile. + + Parameters + ---------- + params : tuple/list + Tuple/list of (surface scale-length, sigma scale-length, sigma(ro)) (note: *not* sigma2 scale-length) + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU) + """ surfaceSigmaProfile.__init__(self) self._params = params def surfacemass(self, R, log=False): """ - NAME: - surfacemass - PURPOSE: - return the surface density profile at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the log (default: False) - OUTPUT: - Sigma(R) - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Return the surface density profile at this R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + Sigma(R) : float + Surface density profile at this R. + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU) """ if log: return -R / self._params[0] @@ -129,17 +151,23 @@ def surfacemass(self, R, log=False): def surfacemassDerivative(self, R, log=False): """ - NAME: - surfacemassDerivative - PURPOSE: - return the derivative wrt R of the surface density profile at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the derivative of the log (default: False) - OUTPUT: - Sigma'(R) or (log Sigma(r) )' - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Return the derivative wrt R of the surface density profile at this R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the derivative of the log (default: False). + + Returns + ------- + float + Sigma'(R) or (log Sigma(r) )'. + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU). """ if log: return -1.0 / self._params[0] @@ -148,17 +176,23 @@ def surfacemassDerivative(self, R, log=False): def sigma2(self, R, log=False): """ - NAME: - sigma2 - PURPOSE: - return the radial velocity variance at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the log (default: False) - OUTPUT: - sigma^2(R) - HISTORY: - 2010-03-26 - Written - Bovy (NYU) + Return the radial velocity variance at this R. + + Parameters + ---------- + R : float + Galactocentric radius (/ro). + log : bool, optional + If True, return the log (default: False). + + Returns + ------- + float + sigma^2(R). + + Notes + ----- + - 2010-03-26 - Written - Bovy (NYU) """ if log: return 2.0 * numpy.log(self._params[2]) - 2.0 * (R - 1.0) / self._params[1] @@ -169,18 +203,25 @@ def sigma2(self, R, log=False): def sigma2Derivative(self, R, log=False): """ - NAME: - sigmaDerivative - PURPOSE: - return the derivative wrt R of the sigma_R^2 profile at this R - INPUT: - R - Galactocentric radius (/ro) - log - if True, return the derivative of the log (default: False) - OUTPUT: - Sigma_R^2'(R) or (log Sigma_R^2(r) )' - HISTORY: - 2011-03-24 - Written - Bovy (NYU) + Return the derivative wrt R of the sigma_R^2 profile at this R + + Parameters + ---------- + R : float + Galactocentric radius (/ro) + log : bool, optional + If True, return the derivative of the log (default: False) + + Returns + ------- + float + Sigma_R^2'(R) or (log Sigma_R^2(r) )' + + Notes + ----- + - 2011-03-24 - Written - Bovy (NYU) """ + if log: return -2.0 / self._params[1] else: