###############################################################################
#
# conversion: utilities to convert from galpy 'natural units' to physical
# units
#
###############################################################################
from functools import wraps
import warnings
import copy
import math as m
from ..util.config import __config__
_APY_UNITS= __config__.getboolean('astropy','astropy-units')
_APY_LOADED= True
try:
from astropy import units, constants
except ImportError:
_APY_UNITS= False
_APY_LOADED= False
_G= 4.302*10.**-3. #pc / Msolar (km/s)^2
_kmsInPcMyr= 1.0227121655399913
_PCIN10p18CM= 3.08567758 #10^18 cm
_CIN10p5KMS= 2.99792458 #10^5 km/s
_MSOLAR10p30KG= 1.9891 #10^30 kg
_EVIN10m19J= 1.60217657 #10^-19 J
else:
_G= constants.G.to(units.pc/units.Msun*units.km**2/units.s**2).value
_kmsInPcMyr= (units.km/units.s).to((units.pc/units.Myr))
_PCIN10p18CM= units.pc.to(units.cm)/10.**18. #10^18 cm
_CIN10p5KMS= constants.c.to((units.km/units.s)).value/10.**5. #10^5 km/s
_MSOLAR10p30KG= units.Msun.to(units.kg)/10.**30. #10^30 kg
_EVIN10m19J= units.eV.to(units.J)*10.**19. #10^-19 J
_MyrIn1013Sec= 3.65242198*0.24*3.6 #use tropical year, like for pms
_TWOPI= 2.*m.pi
[docs]def dens_in_criticaldens(vo,ro,H=70.):
"""
NAME:
dens_in_criticaldens
PURPOSE:
convert density to units of the critical density
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
H= (default: 70) Hubble constant in km/s/Mpc
OUTPUT:
conversion from units where vo=1. at ro=1. to units of the critical density
HISTORY:
2014-01-28 - Written - Bovy (IAS)
"""
return vo**2./ro**2.*10.**6./H**2.*8.*m.pi/3.
[docs]def dens_in_meanmatterdens(vo,ro,H=70.,Om=0.3):
"""
NAME:
dens_in_meanmatterdens
PURPOSE:
convert density to units of the mean matter density
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
H= (default: 70) Hubble constant in km/s/Mpc
Om= (default: 0.3) Omega matter
OUTPUT:
conversion from units where vo=1. at ro=1. to units of the mean matter density
HISTORY:
2014-01-28 - Written - Bovy (IAS)
"""
return dens_in_criticaldens(vo,ro,H=H)/Om
[docs]def dens_in_gevcc(vo,ro):
"""
NAME:
dens_in_gevcc
PURPOSE:
convert density to GeV / cm^3
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1. to GeV/cm^3
HISTORY:
2014-06-16 - Written - Bovy (IAS)
"""
return vo**2./ro**2./_G*_MSOLAR10p30KG*_CIN10p5KMS**2./_EVIN10m19J/ _PCIN10p18CM**3.*10.**-4.
[docs]def dens_in_msolpc3(vo,ro):
"""
NAME:
dens_in_msolpc3
PURPOSE:
convert density to Msolar / pc^3
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1. to Msolar/pc^3
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2./ro**2./_G*10.**-6.
[docs]def force_in_2piGmsolpc2(vo,ro):
"""
NAME:
force_in_2piGmsolpc2
PURPOSE:
convert a force or acceleration to 2piG x Msolar / pc^2
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2./ro/_G*10.**-3./_TWOPI
[docs]def force_in_pcMyr2(vo,ro):
"""
NAME:
force_in_pcMyr2
PURPOSE:
convert a force or acceleration to pc/Myr^2
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2./ro*_kmsInPcMyr**2.*10.**-3.
[docs]def force_in_kmsMyr(vo,ro):
"""
NAME:
force_in_kmsMyr
PURPOSE:
convert a force or acceleration to km/s/Myr
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2./ro*_kmsInPcMyr*10.**-3.
[docs]def force_in_10m13kms2(vo,ro):
"""
NAME:
force_in_10m13kms2
PURPOSE:
convert a force or acceleration to 10^(-13) km/s^2
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2014-01-22 - Written - Bovy (IAS)
"""
return vo**2./ro*_kmsInPcMyr*10.**-3./_MyrIn1013Sec
[docs]def freq_in_Gyr(vo,ro):
"""
NAME:
freq_in_Gyr
PURPOSE:
convert a frequency to 1/Gyr
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo/ro*_kmsInPcMyr
[docs]def freq_in_kmskpc(vo,ro):
"""
NAME:
freq_in_kmskpc
PURPOSE:
convert a frequency to km/s/kpc
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo/ro
[docs]def surfdens_in_msolpc2(vo,ro):
"""
NAME:
surfdens_in_msolpc2
PURPOSE:
convert a surface density to Msolar / pc^2
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2./ro/_G*10.**-3.
[docs]def mass_in_msol(vo,ro):
"""
NAME:
mass_in_msol
PURPOSE:
convert a mass to Msolar
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2.*ro/_G*10.**3.
[docs]def mass_in_1010msol(vo,ro):
"""
NAME:
mass_in_1010msol
PURPOSE:
convert a mass to 10^10 x Msolar
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return vo**2.*ro/_G*10.**-7.
[docs]def time_in_Gyr(vo,ro):
"""
NAME:
time_in_Gyr
PURPOSE:
convert a time to Gyr
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2013-09-01 - Written - Bovy (IAS)
"""
return ro/vo/_kmsInPcMyr
[docs]def velocity_in_kpcGyr(vo,ro):
"""
NAME:
velocity_in_kpcGyr
PURPOSE:
convert a velocity to kpc/Gyr
INPUT:
vo - velocity unit in km/s
ro - length unit in kpc
OUTPUT:
conversion from units where vo=1. at ro=1.
HISTORY:
2014-12-19 - Written - Bovy (IAS)
"""
return vo*_kmsInPcMyr
[docs]def get_physical(obj,include_set=False):
"""
NAME:
get_physical
PURPOSE:
return the velocity and length units for converting between physical and internal units as a dictionary for any galpy object, so they can easily be fed to galpy routines
INPUT:
obj - a galpy object or list of such objects (e.g., a Potential, list of Potentials, Orbit, actionAngle instance, DF instance)
include_set= (False) if True, also include roSet and voSet, flags of whether the unit is explicitly set in the object
OUTPUT:
Dictionary {'ro':length unit in kpc,'vo':velocity unit in km/s}; note that this routine will *always* return these conversion units, even if the obj you provide does not have units turned on
HISTORY:
2019-08-03 - Written - Bovy (UofT)
"""
# Try flattening the object in case it's a nested list of Potentials
from ..potential import flatten as flatten_pot
from ..potential import Force, planarPotential, linearPotential
try:
new_obj= flatten_pot(obj)
except: # pragma: no cover
pass # hope for the best!
else: # only apply flattening for potentials
if isinstance(new_obj,(Force,planarPotential,linearPotential)) \
or (isinstance(new_obj,list) and len(new_obj) > 0 \
and isinstance(new_obj[0],(Force,planarPotential,linearPotential))):
obj= new_obj
if isinstance(obj,list):
out_obj= obj[0]
else:
out_obj= obj
out= {'ro':out_obj._ro,'vo':out_obj._vo}
if include_set:
out.update({'roSet':out_obj._roSet,'voSet':out_obj._voSet})
return out
[docs]def physical_compatible(obj,other_obj):
"""
NAME:
physical_compatible
PURPOSE:
test whether the velocity and length units for converting between physical and internal units are compatible for two galpy objects
INPUT:
obj - a galpy object or list of such objects (e.g., a Potential, list of Potentials, Orbit, actionAngle instance, DF instance)
other_obj - another galpy object or list of such objects (e.g., a Potential, list of Potentials, Orbit, actionAngle instance, DF instance)
OUTPUT:
True if the units are compatible, False if not (compatible means that the units are the same when they are set for both objects)
HISTORY:
2020-04-22 - Written - Bovy (UofT)
"""
if obj is None or other_obj is None: # if one is None, just state compat
return True
phys= get_physical(obj,include_set=True)
other_phys= get_physical(other_obj,include_set=True)
out= True
if phys['roSet'] and other_phys['roSet']:
out= out and m.fabs((phys['ro']-other_phys['ro'])/phys['ro']) < 1e-8
if phys['voSet'] and other_phys['voSet']:
out= out and m.fabs((phys['vo']-other_phys['vo'])/phys['vo']) < 1e-8
return out
# Parsers of different inputs with units
def parse_length(x,ro=None,vo=None):
return x.to(units.kpc).value/ro \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_length_kpc(x):
return x.to(units.kpc).value \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_velocity(x,ro=None,vo=None):
return x.to(units.km/units.s).value/vo \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_velocity_kms(x):
return x.to(units.km/units.s).value \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_angle(x):
return x.to(units.rad).value \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_time(x,ro=None,vo=None):
return x.to(units.Gyr).value/time_in_Gyr(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_mass(x,ro=None,vo=None):
try:
return x.to(units.pc*units.km**2/units.s**2).value\
/mass_in_msol(vo,ro)/_G \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
except units.UnitConversionError: pass
return x.to(1e10*units.Msun).value/mass_in_1010msol(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_energy(x,ro=None,vo=None):
return x.to(units.km**2/units.s**2).value/vo**2. \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_angmom(x,ro=None,vo=None):
return x.to(units.kpc*units.km/units.s).value/ro/vo \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_frequency(x,ro=None,vo=None):
return x.to(units.km/units.s/units.kpc).value/\
freq_in_kmskpc(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_force(x,ro=None,vo=None):
try:
return x.to(units.pc/units.Myr**2).value\
/force_in_pcMyr2(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
except units.UnitConversionError: pass
return x.to(units.Msun/units.pc**2).value\
/force_in_2piGmsolpc2(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_dens(x,ro=None,vo=None):
try:
return x.to(units.Msun/units.pc**3).value\
/dens_in_msolpc3(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
except units.UnitConversionError: pass
# Try Gxdens
return x.to(units.km**2/units.s**2/units.pc**2).value\
/dens_in_msolpc3(vo,ro)/_G \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_surfdens(x,ro=None,vo=None):
try:
return x.to(units.Msun/units.pc**2).value\
/surfdens_in_msolpc2(vo,ro) \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
except units.UnitConversionError: pass
# Try Gxsurfdens
return x.to(units.km**2/units.s**2/units.pc).value\
/surfdens_in_msolpc2(vo,ro)/_G \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
def parse_numdens(x,ro=None,vo=None):
return x.to(1/units.kpc**3).value*ro**3 \
if _APY_LOADED and isinstance(x,units.Quantity) \
else x
#Decorator to apply these transformations
# NOTE: names with underscores in them signify return values that *always* have
# units, which is depended on in the Orbit returns (see issue #326)
_roNecessary= {'time': True,
'position': True,
'position_kpc': True,
'velocity': False,
'velocity2': False,
'velocity2surfacendensity': False,
'velocity_kms': False,
'energy': False,
'density': True,
'numberdensity': True,
'force': True,
'velocity2surfacedensity': True,
'surfacedensity': True,
'numbersurfacedensity': True,
'surfacedensitydistance': True,
'mass': True,
'action': True,
'frequency':True,
'frequency-kmskpc':True,
'forcederivative':True,
'angle':True,
'angle_deg':True,
'proper-motion_masyr':True,
'phasespacedensity':True,
'phasespacedensity2d':True,
'phasespacedensityvelocity':True,
'phasespacedensityvelocity2':True,
'dimensionless':False}
_voNecessary= copy.copy(_roNecessary)
_voNecessary['position']= False
_voNecessary['position_kpc']= False
_voNecessary['angle']= False
_voNecessary['angle_deg']= False
_voNecessary['velocity']= True
_voNecessary['velocity2']= True
_voNecessary['velocity_kms']= True
_voNecessary['energy']= True
def physical_conversion(quantity,pop=False):
"""Decorator to convert to physical coordinates:
quantity = [position,velocity,time]"""
def wrapper(method):
@wraps(method)
def wrapped(*args,**kwargs):
use_physical= kwargs.get('use_physical',True) and \
not kwargs.get('log',False)
# Parse whether ro or vo should be considered to be set, because
# the return value will have units anyway
# (like in Orbit methods that return numbers with units, like ra)
roSet= '_' in quantity # _ in quantity name means always units
voSet= '_' in quantity # _ in quantity name means always units
use_physical= use_physical or '_' in quantity # _ in quantity name means always units
ro= kwargs.get('ro',None)
if ro is None and \
(roSet or (hasattr(args[0],'_roSet') and args[0]._roSet)):
ro= args[0]._ro
if ro is None and isinstance(args[0],list) \
and hasattr(args[0][0],'_roSet') and args[0][0]._roSet:
# For lists of Potentials
ro= args[0][0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
vo= kwargs.get('vo',None)
if vo is None and \
(voSet or (hasattr(args[0],'_voSet') and args[0]._voSet)):
vo= args[0]._vo
if vo is None and isinstance(args[0],list) \
and hasattr(args[0][0],'_voSet') and args[0][0]._voSet:
# For lists of Potentials
vo= args[0][0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
# Override Quantity output?
_apy_units= kwargs.get('quantity',_APY_UNITS)
#Remove ro, vo, use_physical, and quantity kwargs if necessary
if pop and 'use_physical' in kwargs: kwargs.pop('use_physical')
if pop and 'quantity' in kwargs: kwargs.pop('quantity')
if pop and 'ro' in kwargs: kwargs.pop('ro')
if pop and 'vo' in kwargs: kwargs.pop('vo')
if use_physical and \
not (_voNecessary[quantity.lower()] and vo is None) and \
not (_roNecessary[quantity.lower()] and ro is None):
from ..orbit import Orbit
if quantity.lower() == 'time':
fac= time_in_Gyr(vo,ro)
if _apy_units:
u= units.Gyr
elif quantity.lower() == 'position':
fac= ro
if _apy_units:
u= units.kpc
elif quantity.lower() == 'position_kpc': # already in kpc
fac= 1.
if _apy_units:
u= units.kpc
elif quantity.lower() == 'velocity':
fac= vo
if _apy_units:
u= units.km/units.s
elif quantity.lower() == 'velocity2':
fac= vo**2.
if _apy_units:
u= (units.km/units.s)**2
elif quantity.lower() == 'velocity_kms': # already in km/s
fac= 1.
if _apy_units:
u= units.km/units.s
elif quantity.lower() == 'frequency':
if kwargs.get('kmskpc',False) and not _apy_units:
fac= freq_in_kmskpc(vo,ro)
else:
fac= freq_in_Gyr(vo,ro)
if _apy_units:
u= units.Gyr**-1.
elif quantity.lower() == 'frequency-kmskpc':
fac= freq_in_kmskpc(vo,ro)
if _apy_units:
u= units.km/units.s/units.kpc
elif quantity.lower() == 'action':
fac= ro*vo
if _apy_units:
u= units.kpc*units.km/units.s
elif quantity.lower() == 'energy':
fac= vo**2.
if _apy_units:
u= units.km**2./units.s**2.
elif quantity.lower() == 'angle': # in rad
fac= 1.
if _apy_units:
u= units.rad
elif quantity.lower() == 'angle_deg': # already in deg
fac= 1.
if _apy_units:
u= units.deg
elif quantity.lower() == 'proper-motion_masyr': # already in mas/yr
fac= 1.
if _apy_units:
u= units.mas/units.yr
elif quantity.lower() == 'force':
fac= force_in_kmsMyr(vo,ro)
if _apy_units:
u= units.km/units.s/units.Myr
elif quantity.lower() == 'density':
fac= dens_in_msolpc3(vo,ro)
if _apy_units:
u= units.Msun/units.pc**3
elif quantity.lower() == 'numberdensity':
fac= 1/ro**3.
if _apy_units:
u= 1/units.kpc**3
elif quantity.lower() == 'velocity2surfacedensity':
fac= surfdens_in_msolpc2(vo,ro)*vo**2
if _apy_units:
u= units.Msun/units.pc**2*(units.km/units.s)**2
elif quantity.lower() == 'surfacedensity':
fac= surfdens_in_msolpc2(vo,ro)
if _apy_units:
u= units.Msun/units.pc**2
elif quantity.lower() == 'numbersurfacedensity':
fac= 1./ro**2.
if _apy_units:
u= 1/units.kpc**2
elif quantity.lower() == 'surfacedensitydistance':
fac= surfdens_in_msolpc2(vo,ro)*ro*1000.
if _apy_units:
u= units.Msun/units.pc
elif quantity.lower() == 'mass':
fac= mass_in_msol(vo,ro)
if _apy_units:
u= units.Msun
elif quantity.lower() == 'forcederivative':
fac= freq_in_Gyr(vo,ro)**2.
if _apy_units:
u= units.Gyr**-2.
elif quantity.lower() == 'phasespacedensity':
fac= 1./vo**3./ro**3.
if _apy_units:
u= 1/(units.km/units.s)**3/units.kpc**3
elif quantity.lower() == 'phasespacedensity2d':
fac= 1./vo**2./ro**2.
if _apy_units:
u= 1/(units.km/units.s)**2/units.kpc**2
elif quantity.lower() == 'phasespacedensityvelocity':
fac= 1./vo**2./ro**3.
if _apy_units:
u= 1/(units.km/units.s)**2/units.kpc**3
elif quantity.lower() == 'phasespacedensityvelocity2':
fac= 1./vo/ro**3.
if _apy_units:
u= 1/(units.km/units.s)/units.kpc**3
elif quantity.lower() == 'dimensionless':
fac= 1.
if _apy_units:
u= units.dimensionless_unscaled
out= method(*args,**kwargs)
if out is None:
return out
if _apy_units:
return units.Quantity(out*fac,unit=u)
else:
return out*fac
else:
return method(*args,**kwargs)
return wrapped
return wrapper
def potential_physical_input(method):
"""Decorator to convert inputs to Potential functions from physical
to internal coordinates"""
@wraps(method)
def wrapper(*args,**kwargs):
from ..potential import flatten as flatten_potential
Pot= flatten_potential(args[0])
ro= kwargs.get('ro',None)
if ro is None and hasattr(Pot,'_ro'):
ro= Pot._ro
if ro is None and isinstance(Pot,list) \
and hasattr(Pot[0],'_ro'):
# For lists of Potentials
ro= Pot[0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
if 't' in kwargs or 'M' in kwargs:
vo= kwargs.get('vo',None)
if vo is None and hasattr(Pot,'_vo'):
vo= Pot._vo
if vo is None and isinstance(Pot,list) \
and hasattr(Pot[0],'_vo'):
# For lists of Potentials
vo= Pot[0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
# Loop through args
newargs= (Pot,)
for ii in range(1,len(args)):
if _APY_LOADED and isinstance(args[ii],units.Quantity):
newargs= newargs+(args[ii].to(units.kpc).value/ro,)
else:
newargs= newargs+(args[ii],)
args= newargs
# phi and t kwargs, also do R, z, and x in case these are given as kwargs
if 'phi' in kwargs and _APY_LOADED \
and isinstance(kwargs['phi'],units.Quantity):
kwargs['phi']= kwargs['phi'].to(units.rad).value
if 't' in kwargs and _APY_LOADED \
and isinstance(kwargs['t'],units.Quantity):
kwargs['t']= kwargs['t'].to(units.Gyr).value\
/time_in_Gyr(vo,ro)
if 'R' in kwargs and _APY_LOADED \
and isinstance(kwargs['R'],units.Quantity):
kwargs['R']= kwargs['R'].to(units.kpc).value/ro
if 'z' in kwargs and _APY_LOADED \
and isinstance(kwargs['z'],units.Quantity):
kwargs['z']= kwargs['z'].to(units.kpc).value/ro
if 'x' in kwargs and _APY_LOADED \
and isinstance(kwargs['x'],units.Quantity):
kwargs['x']= kwargs['x'].to(units.kpc).value/ro
# v kwarg for dissipative forces
if 'v' in kwargs and _APY_LOADED \
and isinstance(kwargs['v'],units.Quantity):
kwargs['v']= kwargs['v'].to(units.km/units.s).value/vo
# Mass kwarg for rtide
if 'M' in kwargs and _APY_LOADED \
and isinstance(kwargs['M'],units.Quantity):
try:
kwargs['M']= kwargs['M'].to(units.Msun).value\
/mass_in_msol(vo,ro)
except units.UnitConversionError:
kwargs['M']= kwargs['M'].to(units.pc*units.km**2/units.s**2)\
.value/mass_in_msol(vo,ro)/_G
# kwargs that come up in quasiisothermaldf
# z done above
if 'dz' in kwargs and _APY_LOADED \
and isinstance(kwargs['dz'],units.Quantity):
kwargs['dz']= kwargs['dz'].to(units.kpc).value/ro
if 'dR' in kwargs and _APY_LOADED \
and isinstance(kwargs['dR'],units.Quantity):
kwargs['dR']= kwargs['dR'].to(units.kpc).value/ro
if 'zmax' in kwargs and _APY_LOADED \
and isinstance(kwargs['zmax'],units.Quantity):
kwargs['zmax']= kwargs['zmax'].to(units.kpc).value/ro
return method(*args,**kwargs)
return wrapper
def physical_conversion_actionAngle(quantity,pop=False):
"""Decorator to convert to physical coordinates for the actionAngle methods:
quantity= call, actionsFreqs, or actionsFreqsAngles (or EccZmaxRperiRap for actionAngleStaeckel)"""
def wrapper(method):
@wraps(method)
def wrapped(*args,**kwargs):
use_physical= kwargs.get('use_physical',True)
ro= kwargs.get('ro',None)
if ro is None and hasattr(args[0],'_roSet') and args[0]._roSet:
ro= args[0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
vo= kwargs.get('vo',None)
if vo is None and hasattr(args[0],'_voSet') and args[0]._voSet:
vo= args[0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
#Remove ro and vo kwargs if necessary
if pop and 'use_physical' in kwargs: kwargs.pop('use_physical')
if pop and 'ro' in kwargs: kwargs.pop('ro')
if pop and 'vo' in kwargs: kwargs.pop('vo')
if use_physical and not vo is None and not ro is None:
out= method(*args,**kwargs)
if 'call' in quantity or 'actions' in quantity:
if 'actions' in quantity and len(out) < 4: # 1D system
fac= [ro*vo]
if _APY_UNITS:
u= [units.kpc*units.km/units.s]
else:
fac= [ro*vo,ro*vo,ro*vo]
if _APY_UNITS:
u= [units.kpc*units.km/units.s,
units.kpc*units.km/units.s,
units.kpc*units.km/units.s]
if 'Freqs' in quantity:
FreqsFac= freq_in_Gyr(vo,ro)
if len(out) < 4: # 1D system
fac.append(FreqsFac)
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.append(Freqsu)
else:
fac.extend([FreqsFac,FreqsFac,FreqsFac])
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.extend([Freqsu,Freqsu,Freqsu])
if 'Angles' in quantity:
if len(out) < 4: # 1D system
fac.append(1.)
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.append(units.rad)
else:
fac.extend([1.,1.,1.])
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.extend([units.rad,units.rad,units.rad])
if 'EccZmaxRperiRap' in quantity:
fac= [1.,ro,ro,ro]
if _APY_UNITS:
u= [1.,
units.kpc,
units.kpc,
units.kpc]
if _APY_UNITS:
newOut= ()
try:
for ii in range(len(out)):
newOut= newOut+(units.Quantity(out[ii]*fac[ii],
unit=u[ii]),)
except TypeError: # happens if out = scalar
newOut= units.Quantity(out*fac[0],unit=u[0])
else:
newOut= ()
try:
for ii in range(len(out)):
newOut= newOut+(out[ii]*fac[ii],)
except TypeError: # happens if out = scalar
newOut= out*fac[0]
return newOut
else:
return method(*args,**kwargs)
return wrapped
return wrapper
def actionAngle_physical_input(method):
"""Decorator to convert inputs to actionAngle functions from physical
to internal coordinates"""
@wraps(method)
def wrapper(*args,**kwargs):
if len(args) < 3: # orbit input
return method(*args,**kwargs)
ro= kwargs.get('ro',None)
if ro is None and hasattr(args[0],'_ro'):
ro= args[0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
vo= kwargs.get('vo',None)
if vo is None and hasattr(args[0],'_vo'):
vo= args[0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
# Loop through args
newargs= ()
for ii in range(len(args)):
if _APY_LOADED and isinstance(args[ii],units.Quantity):
try:
targ= args[ii].to(units.kpc).value/ro
except units.UnitConversionError:
try:
targ= args[ii].to(units.km/units.s).value/vo
except units.UnitConversionError:
try:
targ= args[ii].to(units.rad).value
except units.UnitConversionError:
raise units.UnitConversionError("Input units not understood")
newargs= newargs+(targ,)
else:
newargs= newargs+(args[ii],)
args= newargs
return method(*args,**kwargs)
return wrapper
def physical_conversion_actionAngleInverse(quantity,pop=False):
"""Decorator to convert to physical coordinates for the actionAngleInverse methods:
quantity= call, xvFreqs, or Freqs"""
def wrapper(method):
@wraps(method)
def wrapped(*args,**kwargs):
use_physical= kwargs.get('use_physical',True)
ro= kwargs.get('ro',None)
if ro is None and hasattr(args[0],'_roSet') and args[0]._roSet:
ro= args[0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
vo= kwargs.get('vo',None)
if vo is None and hasattr(args[0],'_voSet') and args[0]._voSet:
vo= args[0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
#Remove ro and vo kwargs if necessary
if pop and 'use_physical' in kwargs: kwargs.pop('use_physical')
if pop and 'ro' in kwargs: kwargs.pop('ro')
if pop and 'vo' in kwargs: kwargs.pop('vo')
if use_physical and not vo is None and not ro is None:
fac= []
u= []
out= method(*args,**kwargs)
if 'call' in quantity or 'xv' in quantity:
if 'xv' in quantity and len(out) < 4: # 1D system
fac.extend([ro,vo])
if _APY_UNITS:
u.extend([units.kpc,units.km/units.s])
else:
fac.extend([ro,vo,vo,ro,vo,1.])
if _APY_UNITS:
u.extend([units.kpc,units.km/units.s,
units.km/units.s,units.kpc,
units.km/units.s,
units.rad])
if 'Freqs' in quantity:
FreqsFac= freq_in_Gyr(vo,ro)
if isinstance(out,float): # 1D system
fac.append(FreqsFac)
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.append(Freqsu)
else:
fac.extend([FreqsFac,FreqsFac,FreqsFac])
if _APY_UNITS:
Freqsu= units.Gyr**-1.
u.extend([Freqsu,Freqsu,Freqsu])
if _APY_UNITS:
newOut= ()
try:
for ii in range(len(out)):
newOut= newOut+(units.Quantity(out[ii]*fac[ii],
unit=u[ii]),)
except TypeError: # Happens when out == scalar
newOut= units.Quantity(out*fac[0],unit=u[0])
else:
newOut= ()
try:
for ii in range(len(out)):
newOut= newOut+(out[ii]*fac[ii],)
except TypeError: # Happens when out == scalar
newOut= out*fac[0]
return newOut
else:
return method(*args,**kwargs)
return wrapped
return wrapper
def actionAngleInverse_physical_input(method):
"""Decorator to convert inputs to actionAngleInverse functions from
physical to internal coordinates"""
@wraps(method)
def wrapper(*args,**kwargs):
ro= kwargs.get('ro',None)
if ro is None and hasattr(args[0],'_ro'):
ro= args[0]._ro
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
vo= kwargs.get('vo',None)
if vo is None and hasattr(args[0],'_vo'):
vo= args[0]._vo
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
# Loop through args
newargs= ()
for ii in range(len(args)):
if _APY_LOADED and isinstance(args[ii],units.Quantity):
try:
targ= args[ii].to(units.kpc*units.km/units.s).value/ro/vo
except units.UnitConversionError:
try:
targ= args[ii].to(units.rad).value
except units.UnitConversionError:
raise units.UnitConversionError("Input units not understood")
newargs= newargs+(targ,)
else:
newargs= newargs+(args[ii],)
args= newargs
return method(*args,**kwargs)
return wrapper