From 3231eb6607e64d43626189015178eec83f13a805 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 May 2019 00:58:14 -0400 Subject: [PATCH] Retire old Orbit interface in favor of new Orbits Orbit interface --- galpy/orbit/FullOrbit.py | 820 -------- galpy/orbit/Orbit.py | 3968 ------------------------------------ galpy/orbit/OrbitTop.py | 2207 -------------------- galpy/orbit/RZOrbit.py | 535 ----- galpy/orbit/linearOrbit.py | 195 -- galpy/orbit/planarOrbit.py | 855 -------- 6 files changed, 8580 deletions(-) delete mode 100644 galpy/orbit/FullOrbit.py delete mode 100644 galpy/orbit/Orbit.py delete mode 100644 galpy/orbit/OrbitTop.py delete mode 100644 galpy/orbit/RZOrbit.py delete mode 100644 galpy/orbit/linearOrbit.py delete mode 100644 galpy/orbit/planarOrbit.py diff --git a/galpy/orbit/FullOrbit.py b/galpy/orbit/FullOrbit.py deleted file mode 100644 index 1cf9a00f5..000000000 --- a/galpy/orbit/FullOrbit.py +++ /dev/null @@ -1,820 +0,0 @@ -import warnings -import math as m -import numpy as nu -from scipy import integrate, optimize -import scipy -_SCIPY_VERSION= [int(v.split('rc')[0]) - for v in scipy.__version__.split('.')] -if _SCIPY_VERSION[0] < 1 and _SCIPY_VERSION[1] < 10: #pragma: no cover - from scipy.maxentropy import logsumexp -elif _SCIPY_VERSION[0] < 1 and _SCIPY_VERSION[1] < 19: #pragma: no cover - from scipy.misc import logsumexp -else: - from scipy.special import logsumexp -from galpy.potential.Potential import _evaluateRforces, _evaluatezforces,\ - evaluatePotentials, _evaluatephiforces, evaluateDensities, _check_c -from galpy.potential.DissipativeForce import _isDissipative -from galpy.util import galpyWarning, galpyWarningVerbose -import galpy.util.bovy_plot as plot -from galpy.util.leung_dop853 import dop853 -import galpy.util.bovy_symplecticode as symplecticode -import galpy.util.bovy_coords as coords -#try: -from .integrateFullOrbit import integrateFullOrbit_c, _ext_loaded -ext_loaded= _ext_loaded -from galpy.util.bovy_conversion import physical_conversion -from .OrbitTop import OrbitTop -_ORBFITNORMRADEC= 360. -_ORBFITNORMDIST= 10. -_ORBFITNORMPMRADEC= 4. -_ORBFITNORMVLOS= 200. -class FullOrbit(OrbitTop): - """Class that holds and integrates orbits in full 3D potentials""" - def __init__(self,vxvv=[1.,0.,0.9,0.,0.1],vo=220.,ro=8.0,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): - """ - NAME: - - __init__ - - PURPOSE: - - intialize a full orbit - - INPUT: - - vxvv - initial condition [R,vR,vT,z,vz,phi] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - (none) - - HISTORY: - - 2010-08-01 - Written - Bovy (NYU) - - 2014-06-11 - Added conversion kwargs to physical coordinates - Bovy (IAS) - - """ - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=zo,vo=vo,solarmotion=solarmotion) - return None - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'dop853' for a dop853 implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - (none) (get the actual orbit using getOrbit() - HISTORY: - 2010-08-01 - Written - Bovy (NYU) - """ - #Reset things that may have been defined by a previous integration - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - self.t= nu.array(t) - self._pot= pot - self.orbit= _integrateFullOrbit(self.vxvv,pot,t,method,dt) - - @physical_conversion('energy') - def Jacobi(self,*args,**kwargs): - """ - NAME: - Jacobi - PURPOSE: - calculate the Jacobi integral of the motion - INPUT: - Omega - pattern speed of rotating frame - t= time - pot= potential instance or list of such instances - OUTPUT: - Jacobi integral - HISTORY: - 2011-04-18 - Written - Bovy (NYU) - """ - if not 'OmegaP' in kwargs or kwargs['OmegaP'] is None: - OmegaP= 1. - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - else: - pot= kwargs['pot'] - if isinstance(pot,list): - for p in pot: - if hasattr(p,'OmegaP'): - OmegaP= p.OmegaP() - break - else: - if hasattr(pot,'OmegaP'): - OmegaP= pot.OmegaP() - kwargs.pop('OmegaP',None) - else: - OmegaP= kwargs.pop('OmegaP') - #Make sure you are not using physical coordinates - old_physical= kwargs.get('use_physical',None) - kwargs['use_physical']= False - if not isinstance(OmegaP,(int,float)) and len(OmegaP) == 3: - if isinstance(OmegaP,list): thisOmegaP= nu.array(OmegaP) - else: thisOmegaP= OmegaP - out= self.E(*args,**kwargs)-nu.dot(thisOmegaP, - self.L(*args,**kwargs).T).T - else: - out= self.E(*args,**kwargs)-OmegaP*self.L(*args,**kwargs)[:,2] - if not old_physical is None: - kwargs['use_physical']= old_physical - else: - kwargs.pop('use_physical') - return out - - @physical_conversion('energy') - def E(self,*args,**kwargs): - """ - NAME: - E - PURPOSE: - calculate the energy - INPUT: - t - (optional) time at which to get the energy - pot= potential instance or list of such instances - OUTPUT: - energy - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],thiso[3], - phi=thiso[5],t=t,use_physical=False)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2.\ - +thiso[4]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],thiso[3,ii], - phi=thiso[5,ii], - t=t[ii],use_physical=False)\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2.\ - +thiso[4,ii]**2./2. for ii in range(len(t))]) - - @physical_conversion('energy') - def ER(self,*args,**kwargs): - """ - NAME: - ER - PURPOSE: - calculate the radial energy - INPUT: - t - (optional) time at which to get the energy - pot= potential instance or list of such instances - OUTPUT: - radial energy - HISTORY: - 2013-11-30 - Written - Bovy (IAS) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],0., - phi=thiso[5],t=t,use_physical=False)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],0., - phi=thiso[5,ii], - t=t[ii],use_physical=False)\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2. for ii in range(len(t))]) - - @physical_conversion('energy') - def Ez(self,*args,**kwargs): - """ - NAME: - Ez - PURPOSE: - calculate the vertical energy - INPUT: - t - (optional) time at which to get the energy - pot= potential instance or list of such instances - OUTPUT: - vertical energy - HISTORY: - 2013-11-30 - Written - Bovy (IAS) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],thiso[3], - phi=thiso[5],t=t,use_physical=False)\ - -evaluatePotentials(pot,thiso[0],0., - phi=thiso[5],t=t, - use_physical=False)\ - +thiso[4]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],thiso[3,ii], - phi=thiso[5,ii], - t=t[ii],use_physical=False)\ - -evaluatePotentials(pot,thiso[0,ii],0., - phi=thiso[5,ii], - t=t[ii],use_physical=False)\ - +thiso[4,ii]**2./2. for ii in range(len(t))]) - - def e(self,analytic=False,pot=None,**kwargs): - """ - NAME: - e - PURPOSE: - calculate the eccentricity - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - eccentricity - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[0]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate eccentricity") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs)) - - @physical_conversion('position') - def rap(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rap - PURPOSE: - return the apocenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_ap - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[3]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate rap") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return nu.amax(self.rs) - - @physical_conversion('position') - def rperi(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rperi - PURPOSE: - return the pericenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_peri - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[2]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate rperi") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return nu.amin(self.rs) - - @physical_conversion('position') - def zmax(self,analytic=False,pot=None,**kwargs): - """ - NAME: - zmax - PURPOSE: - return the maximum vertical height - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - Z_max - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - 2012-06-01 - Added analytic calculation - Bovy (IAS) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[1]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate zmax") - return nu.amax(nu.fabs(self.orbit[:,3])) - - def fit(self,vxvv,vxvv_err=None,pot=None,radec=False,lb=False, - customsky=False,lb_to_customsky=None,pmllpmbb_to_customsky=None, - tintJ=10,ntintJ=1000,integrate_method='dopr54_c', - disp=False, - **kwargs): - """ - NAME: - fit - PURPOSE: - fit an Orbit to data using the current orbit as the initial - condition - INPUT: - vxvv - [:,6] array of positions and velocities along the orbit [cannot be Quantities] - vxvv_err= [:,6] array of errors on positions and velocities along the orbit (if None, these are set to 0.01) [cannot be Quantities] - pot= Potential to fit the orbit in - - Keywords related to the input data: - radec= if True, input vxvv and vxvv_err are [ra,dec,d,mu_ra, mu_dec,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (all ICRS; mu_ra = mu_ra * cos dec); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates; Note that for speed reasons, galpy's internal transformation between (l,b) and (ra,dec) is used, rather than astropy's - lb= if True, input vxvv and vxvv_err are [long,lat,d,mu_ll, mu_bb,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (mu_ll = mu_ll * cos lat); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates - customsky= if True, input vxvv and vxvv_err are [custom long,custom lat,d,mu_customll, mu_custombb,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (mu_ll = mu_ll * cos lat) where custom longitude and custom latitude are a custom set of sky coordinates (e.g., ecliptic) and the proper motions are also expressed in these coordinats; you need to provide the functions lb_to_customsky and pmllpmbb_to_customsky to convert to the custom sky coordinates (these should have the same inputs and outputs as lb_to_radec and pmllpmbb_to_pmrapmdec); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - Cannot be an Orbit instance with the orbit of the reference point, as w/ the ra etc. functions - ro= distance in kpc corresponding to R=1. (default: taken from object) - vo= velocity in km/s corresponding to v=1. (default: taken from object) - - lb_to_customsky= function that converts l,b,degree=False to the custom sky coordinates (like lb_to_radec); needs to be given when customsky=True - - pmllpmbb_to_customsky= function that converts pmll,pmbb,l,b,degree=False to proper motions in the custom sky coordinates (like pmllpmbb_to_pmrapmdec); needs to be given when customsky=True - - Keywords related to the orbit integrations: - tintJ= (default: 10) time to integrate orbits for fitting the orbit - ntintJ= (default: 1000) number of time-integration points - integrate_method= (default: 'dopr54_c') integration method to use - disp= (False) display the optimizer's convergence message - - OUTPUT: - max of log likelihood - HISTORY: - 2014-06-17 - Written - Bovy (IAS) - - TEST: - from galpy.potential import LogarithmicHaloPotential; lp= LogarithmicHaloPotential(normalize=1.); from galpy.orbit import Orbit; o= Orbit(vxvv=[1.,0.1,1.1,0.1,0.02,0.]); ts= numpy.linspace(0,10,1000); o.integrate(ts,lp); outts= [0.,0.1,0.2,0.3,0.4]; vxvv= numpy.array([o.R(outts),o.vR(outts),o.vT(outts),o.z(outts),o.vz(outts),o.phi(outts)]).T; of= Orbit(vxvv=[1.02,0.101,1.101,0.101,0.0201,0.001]); of._orb.fit(vxvv,pot=lp,radec=False,tintJ=10,ntintJ=1000) - - """ - if pot is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit first or specify pot=") - if radec or lb or customsky: - from .Orbits import _parse_radec_kwargs - obs, ro, vo= _parse_radec_kwargs(self,kwargs,vel=True,dontpop=True) - else: - obs, ro, vo= None, None, None - if customsky \ - and (lb_to_customsky is None or pmllpmbb_to_customsky is None): - raise IOError('if customsky=True, the functions lb_to_customsky and pmllpmbb_to_customsky need to be given') - new_vxvv, maxLogL= _fit_orbit(self,vxvv,vxvv_err,pot,radec=radec,lb=lb, - customsky=customsky, - lb_to_customsky=lb_to_customsky, - pmllpmbb_to_customsky=pmllpmbb_to_customsky, - tintJ=tintJ,ntintJ=ntintJ, - integrate_method=integrate_method, - ro=ro,vo=vo,obs=obs,disp=disp) - #Setup with these new initial conditions - self.vxvv= new_vxvv - return maxLogL - - def plotEz(self,*args,**kwargs): - """ - NAME: - plotEz - PURPOSE: - plot Ez(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'Eznorm' - else: - kwargs['d2']= 'Ez' - return self.plot(*args,**kwargs) - - def plotER(self,*args,**kwargs): - """ - NAME: - plotER - PURPOSE: - plot ER(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'ERnorm' - else: - kwargs['d2']= 'ER' - return self.plot(*args,**kwargs) - - def plotEzJz(self,*args,**kwargs): - """ - NAME: - plotEzJz - PURPOSE: - plot E_z(.)/sqrt(dens(R)) along the orbit - INPUT: - pot= Potential instance or list of instances in which the orbit was - integrated - d1= - plot Ez vs d1: e.g., 't', 'z', 'R', 'vR', 'vT', 'vz' - +bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-08-08 - Written - Bovy (NYU) - """ - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'x':r'$x$','y':r'$y$','vx':r'$v_x$','vy':r'$v_y$'} - if not 'pot' in kwargs: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit first or specify pot=") - else: - pot= kwargs.pop('pot') - d1= kwargs.pop('d1','t') - self.EzJz= [(evaluatePotentials(pot,self.orbit[ii,0],self.orbit[ii,3], - t=self.t[ii],use_physical=False)- - evaluatePotentials(pot,self.orbit[ii,0],0., - phi= self.orbit[ii,5],t=self.t[ii], - use_physical=False)+ - self.orbit[ii,4]**2./2.)/\ - nu.sqrt(evaluateDensities(pot,self.orbit[ii,0],0., - phi=self.orbit[ii,5], - t=self.t[ii], - use_physical=False))\ - for ii in range(len(self.t))] - if not 'xlabel' in kwargs: - kwargs['xlabel']= labeldict[d1] - if not 'ylabel' in kwargs: - kwargs['ylabel']= r'$E_z/\sqrt{\rho}$' - if d1 == 't': - return plot.bovy_plot(nu.array(self.t), - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'z': - return plot.bovy_plot(self.orbit[:,3], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'R': - return plot.bovy_plot(self.orbit[:,0], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vR': - return plot.bovy_plot(self.orbit[:,1], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vT': - return plot.bovy_plot(self.orbit[:,2], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vz': - return plot.bovy_plot(self.orbit[:,4], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'phi': - return plot.bovy_plot(self.orbit[:,5], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - -def _integrateFullOrbit(vxvv,pot,t,method,dt): - """ - NAME: - _integrateFullOrbit - PURPOSE: - integrate an orbit in a Phi(R,z,phi) potential - INPUT: - vxvv - array with the initial conditions stacked like - [R,vR,vT,z,vz,phi]; vR outward! - pot - Potential instance - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - dt - if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - [:,5] array of [R,vR,vT,z,vz,phi] at each t - HISTORY: - 2010-08-01 - Written - Bovy (NYU) - """ - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - # Now check that we aren't trying to integrate a dissipative force - # with a symplectic integrator - if _isDissipative(pot) and ('leapfrog' in method or 'symplec' in method): - if '_c' in method: - method= 'dopr54_c' - else: - method= 'odeint' - warnings.warn("Cannot use symplectic integration because some of the included forces are dissipative (using non-symplectic integrator %s instead)" % (method), galpyWarning) - if method.lower() == 'leapfrog': - #go to the rectangular frame - this_vxvv= nu.array([vxvv[0]*nu.cos(vxvv[5]), - vxvv[0]*nu.sin(vxvv[5]), - vxvv[3], - vxvv[1]*nu.cos(vxvv[5])-vxvv[2]*nu.sin(vxvv[5]), - vxvv[2]*nu.cos(vxvv[5])+vxvv[1]*nu.sin(vxvv[5]), - vxvv[4]]) - #integrate - out= symplecticode.leapfrog(_rectForce,this_vxvv, - t,args=(pot,),rtol=10.**-8) - #go back to the cylindrical frame - R= nu.sqrt(out[:,0]**2.+out[:,1]**2.) - phi= nu.arccos(out[:,0]/R) - phi[(out[:,1] < 0.)]= 2.*nu.pi-phi[(out[:,1] < 0.)] - vR= out[:,3]*nu.cos(phi)+out[:,4]*nu.sin(phi) - vT= out[:,4]*nu.cos(phi)-out[:,3]*nu.sin(phi) - out[:,3]= out[:,2] - out[:,4]= out[:,5] - out[:,0]= R - out[:,1]= vR - out[:,2]= vT - out[:,5]= phi - elif ext_loaded and \ - (method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \ - or method.lower() == 'rk6_c' or method.lower() == 'symplec4_c' \ - or method.lower() == 'symplec6_c' or method.lower() == 'dopr54_c' \ - or method.lower() == 'dop853_c'): - warnings.warn("Using C implementation to integrate orbits", - galpyWarningVerbose) - out, msg= integrateFullOrbit_c(pot,nu.copy(vxvv),t,method,dt=dt) - elif method.lower() == 'odeint' or method.lower() == 'dop853' or not ext_loaded: - vphi= vxvv[2]/vxvv[0] - init= [vxvv[0],vxvv[1],vxvv[5],vphi,vxvv[3],vxvv[4]] - if method == 'dop853': - intOut = dop853(_FullEOM, init, t, args=(pot,)) - else: - intOut= integrate.odeint(_FullEOM,init,t,args=(pot,), - rtol=10.**-8.)#,mxstep=100000000) - out= nu.zeros((len(t),6)) - out[:,0]= intOut[:,0] - out[:,1]= intOut[:,1] - out[:,2]= out[:,0]*intOut[:,3] - out[:,3]= intOut[:,4] - out[:,4]= intOut[:,5] - out[:,5]= intOut[:,2] - #post-process to remove negative radii - neg_radii= (out[:,0] < 0.) - out[neg_radii,0]= -out[neg_radii,0] - out[neg_radii,5]+= m.pi - return out - -def _FullEOM(y,t,pot): - """ - NAME: - _FullEOM - PURPOSE: - implements the EOM, i.e., the right-hand side of the differential - equation - INPUT: - y - current phase-space position - t - current time - pot - (list of) Potential instance(s) - OUTPUT: - dy/dt - HISTORY: - 2010-04-16 - Written - Bovy (NYU) - """ - l2= (y[0]**2.*y[3])**2. - return [y[1], - l2/y[0]**3.+_evaluateRforces(pot,y[0],y[4],phi=y[2],t=t, - v=[y[1],y[0]*y[3],y[5]]), - y[3], - 1./y[0]**2.*(_evaluatephiforces(pot,y[0],y[4],phi=y[2],t=t, - v=[y[1],y[0]*y[3],y[5]]) - -2.*y[0]*y[1]*y[3]), - y[5], - _evaluatezforces(pot,y[0],y[4],phi=y[2],t=t, - v=[y[1],y[0]*y[3],y[5]])] - -def _rectForce(x,pot,t=0.): - """ - NAME: - _rectForce - PURPOSE: - returns the force in the rectangular frame - INPUT: - x - current position - t - current time - pot - (list of) Potential instance(s) - OUTPUT: - force - HISTORY: - 2011-02-02 - Written - Bovy (NYU) - """ - #x is rectangular so calculate R and phi - R= nu.sqrt(x[0]**2.+x[1]**2.) - phi= nu.arccos(x[0]/R) - sinphi= x[1]/R - cosphi= x[0]/R - if x[1] < 0.: phi= 2.*nu.pi-phi - #calculate forces - Rforce= _evaluateRforces(pot,R,x[2],phi=phi,t=t) - phiforce= _evaluatephiforces(pot,R,x[2],phi=phi,t=t) - return nu.array([cosphi*Rforce-1./R*sinphi*phiforce, - sinphi*Rforce+1./R*cosphi*phiforce, - _evaluatezforces(pot,R,x[2],phi=phi,t=t)]) - -def _fit_orbit(orb,vxvv,vxvv_err,pot,radec=False,lb=False, - customsky=False,lb_to_customsky=None, - pmllpmbb_to_customsky=None, - tintJ=100,ntintJ=1000,integrate_method='dopr54_c', - ro=None,vo=None,obs=None,disp=False): - """Fit an orbit to data in a given potential""" - # Need to turn this off for speed - coords._APY_COORDS_ORIG= coords._APY_COORDS - coords._APY_COORDS= False - #Import here, because otherwise there is an infinite loop of imports - from galpy.actionAngle import actionAngleIsochroneApprox, actionAngle - #Mock this up, bc we want to use its orbit-integration routines - class mockActionAngleIsochroneApprox(actionAngleIsochroneApprox): - def __init__(self,tintJ,ntintJ,pot,integrate_method='dopr54_c'): - actionAngle.__init__(self) - self._tintJ= tintJ - self._ntintJ=ntintJ - self._tsJ= nu.linspace(0.,self._tintJ,self._ntintJ) - self._integrate_dt= None - self._pot= pot - self._integrate_method= integrate_method - return None - tmockAA= mockActionAngleIsochroneApprox(tintJ,ntintJ,pot, - integrate_method=integrate_method) - opt_vxvv= optimize.fmin_powell(_fit_orbit_mlogl,orb.vxvv, - args=(vxvv,vxvv_err,pot,radec,lb, - customsky,lb_to_customsky, - pmllpmbb_to_customsky, - tmockAA, - ro,vo,obs), - disp=disp) - maxLogL= -_fit_orbit_mlogl(opt_vxvv,vxvv,vxvv_err,pot,radec,lb, - customsky,lb_to_customsky,pmllpmbb_to_customsky, - tmockAA, - ro,vo,obs) - coords._APY_COORDS= coords._APY_COORDS_ORIG - return (opt_vxvv,maxLogL) - -def _fit_orbit_mlogl(new_vxvv,vxvv,vxvv_err,pot,radec,lb, - customsky,lb_to_customsky,pmllpmbb_to_customsky, - tmockAA, - ro,vo,obs): - """The log likelihood for fitting an orbit""" - #Use this _parse_args routine, which does forward and backward integration - iR,ivR,ivT,iz,ivz,iphi= tmockAA._parse_args(True,False, - new_vxvv[0], - new_vxvv[1], - new_vxvv[2], - new_vxvv[3], - new_vxvv[4], - new_vxvv[5]) - if radec or lb or customsky: - #Need to transform to (l,b), (ra,dec), or a custom set - #First transform to X,Y,Z,vX,vY,vZ (Galactic) - X,Y,Z = coords.galcencyl_to_XYZ(iR.flatten(),iphi.flatten(), - iz.flatten(), - Xsun=obs[0]/ro, - Zsun=obs[2]/ro).T - vX,vY,vZ = coords.galcencyl_to_vxvyvz(ivR.flatten(),ivT.flatten(), - ivz.flatten(),iphi.flatten(), - vsun=nu.array(\ - obs[3:6])/vo,Xsun=obs[0]/ro,Zsun=obs[2]/ro).T - bad_indx= (X == 0.)*(Y == 0.)*(Z == 0.) - if True in bad_indx: X[bad_indx]+= ro/10000. - lbdvrpmllpmbb= coords.rectgal_to_sphergal(X*ro,Y*ro,Z*ro, - vX*vo,vY*vo,vZ*vo, - degree=True) - if lb: - orb_vxvv= nu.array([lbdvrpmllpmbb[:,0], - lbdvrpmllpmbb[:,1], - lbdvrpmllpmbb[:,2], - lbdvrpmllpmbb[:,4], - lbdvrpmllpmbb[:,5], - lbdvrpmllpmbb[:,3]]).T - elif radec: - #Further transform to ra,dec,pmra,pmdec - radec= coords.lb_to_radec(lbdvrpmllpmbb[:,0], - lbdvrpmllpmbb[:,1],degree=True, - epoch=None) - pmrapmdec= coords.pmllpmbb_to_pmrapmdec(lbdvrpmllpmbb[:,4], - lbdvrpmllpmbb[:,5], - lbdvrpmllpmbb[:,0], - lbdvrpmllpmbb[:,1], - degree=True, - epoch=None) - orb_vxvv= nu.array([radec[:,0],radec[:,1], - lbdvrpmllpmbb[:,2], - pmrapmdec[:,0],pmrapmdec[:,1], - lbdvrpmllpmbb[:,3]]).T - elif customsky: - #Further transform to ra,dec,pmra,pmdec - customradec= lb_to_customsky(lbdvrpmllpmbb[:,0], - lbdvrpmllpmbb[:,1],degree=True) - custompmrapmdec= pmllpmbb_to_customsky(lbdvrpmllpmbb[:,4], - lbdvrpmllpmbb[:,5], - lbdvrpmllpmbb[:,0], - lbdvrpmllpmbb[:,1], - degree=True) - orb_vxvv= nu.array([customradec[:,0],customradec[:,1], - lbdvrpmllpmbb[:,2], - custompmrapmdec[:,0],custompmrapmdec[:,1], - lbdvrpmllpmbb[:,3]]).T - else: - #shape=(2tintJ-1,6) - orb_vxvv= nu.array([iR.flatten(),ivR.flatten(),ivT.flatten(), - iz.flatten(),ivz.flatten(),iphi.flatten()]).T - out= 0. - for ii in range(vxvv.shape[0]): - sub_vxvv= (orb_vxvv-vxvv[ii,:].flatten())**2. - #print(sub_vxvv[nu.argmin(nu.sum(sub_vxvv,axis=1))]) - if not vxvv_err is None: - sub_vxvv/= vxvv_err[ii,:]**2. - else: - sub_vxvv/= 0.01**2. - out+= logsumexp(-0.5*nu.sum(sub_vxvv,axis=1)) - return -out - diff --git a/galpy/orbit/Orbit.py b/galpy/orbit/Orbit.py deleted file mode 100644 index 32b06a937..000000000 --- a/galpy/orbit/Orbit.py +++ /dev/null @@ -1,3968 +0,0 @@ -import warnings -import numpy as nu -_APY_LOADED= True -try: - from astropy import units -except ImportError: - _APY_LOADED= False -if _APY_LOADED: - import astropy - _APY3= astropy.__version__ > '3' - _APY_GE_31= tuple(map(int,(astropy.__version__.split('.')))) > (3,0,5) - from astropy.coordinates import SkyCoord, Galactocentric, \ - CartesianDifferential -_ASTROQUERY_LOADED= True -try: - from astroquery.simbad import Simbad -except ImportError: - _ASTROQUERY_LOADED= False -import galpy.util.bovy_coords as coords -from galpy.util.bovy_conversion import physical_conversion -from galpy.util import galpyWarning -from galpy.util import bovy_conversion -from galpy.util import config -_APY_UNITS= config.__config__.getboolean('astropy','astropy-units') -from .Orbits import _check_integrate_dt, _check_potential_dim, \ - _check_consistent_units -from .FullOrbit import FullOrbit -from .RZOrbit import RZOrbit -from .planarOrbit import planarOrbit, planarROrbit, \ - planarOrbitTop -from .linearOrbit import linearOrbit -from galpy.potential import flatten as flatten_potential -from galpy.potential import rl, _isNonAxi -from ..util.bovy_coords import _K -if _APY_LOADED: - vxvv_units= [units.kpc,units.km/units.s,units.km/units.s, - units.kpc,units.km/units.s,units.rad] -else: - _APY_UNITS= False -class Orbit(object): - """General orbit class representing an orbit""" - def __init__(self,vxvv=None,uvw=False,lb=False, - radec=False,vo=None,ro=None,zo=None, - solarmotion=None): - """ - NAME: - - __init__ - - PURPOSE: - - Initialize an Orbit instance - - INPUT: - - vxvv - initial conditions - 3D can be either - - 1) in Galactocentric cylindrical coordinates [R,vR,vT(,z,vz,phi)]; can be Quantities - - 2) astropy (>v3.0) SkyCoord that includes velocities (Note that this turns *on* physical output even if ro and vo are not given) - - 3) [ra,dec,d,mu_ra, mu_dec,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (all J2000.0; mu_ra = mu_ra * cos dec); can be Quantities; ICRS frame - - 4) [ra,dec,d,U,V,W] in [deg,deg,kpc,km/s,km/s,kms]; can be Quantities; ICRS frame - - 5) (l,b,d,mu_l, mu_b, vlos) in [deg,deg,kpc,mas/yr,mas/yr,km/s) (all J2000.0; mu_l = mu_l * cos b); can be Quantities - - 6) [l,b,d,U,V,W] in [deg,deg,kpc,km/s,km/s,kms]; can be Quantities - - 7) Unspecified: assumed to be the Sun (equivalent to ``vxvv= [0,0,0,0,0,0]`` and ``radec=True``) - - 5) and 6) also work when leaving out b and mu_b/W - - OPTIONAL INPUTS: - - radec= if True, input is 2) (or 3) above (Note that this turns *on* physical output even if ro and vo are not given) - - uvw= if True, velocities are UVW - - lb= if True, input is 4) or 5) above (Note that this turns *on* physical output even if ro and vo are not given) - - ro= distance from vantage point to GC (kpc; can be Quantity) - - vo= circular velocity at ro (km/s; can be Quantity) - - zo= offset toward the NGP of the Sun wrt the plane (kpc; can be Quantity; default = 25 pc) - - solarmotion= 'hogg' or 'dehnen', or 'schoenrich', or value in - [-U,V,W]; can be Quantity - - If ro and/or vo are specified, outputs involving distances or velocities (whether as instance methods or in plots) will by default be displayed in the physical coordinates implied by these scales. This can be overwritten for each individual method by using use_physical=False as a keyword for the method. - - OUTPUT: - - instance - - HISTORY: - - 2010-07-20 - Written - Bovy (NYU) - - """ - # If you change the way an Orbit object is setup, also change each of - # the methods that return Orbits - if vxvv is None: # Assume one wants the Sun - vxvv= [0.,0.,0.,0.,0.,0.] - radec= True - if _APY_LOADED and isinstance(ro,units.Quantity): - ro= ro.to(units.kpc).value - if _APY_LOADED and isinstance(zo,units.Quantity): - zo= zo.to(units.kpc).value - if _APY_LOADED and isinstance(vo,units.Quantity): - vo= vo.to(units.km/units.s).value - # if vxvv is SkyCoord, preferentially use its ro and zo - if _APY_LOADED and isinstance(vxvv,SkyCoord): - if not _APY3: # pragma: no cover - raise ImportError('Orbit initialization using an astropy SkyCoord requires astropy >3.0') - if zo is None and not vxvv.z_sun is None: - zo= vxvv.z_sun.to(units.kpc).value - elif not vxvv.z_sun is None: - if nu.fabs(zo-vxvv.z_sun.to(units.kpc).value) > 1e-8: - raise ValueError("Orbit initialization's zo different from SkyCoord's z_sun; these should be the same for consistency") - elif zo is None and not vxvv.galcen_distance is None: - zo= 0. - if ro is None and not vxvv.galcen_distance is None: - ro= nu.sqrt(vxvv.galcen_distance.to(units.kpc).value**2. - -zo**2.) - elif not vxvv.galcen_distance is None and \ - nu.fabs(ro**2.+zo**2.-vxvv.galcen_distance.to(units.kpc).value**2.) > 1e-10: - warnings.warn("Orbit's initialization normalization ro and zo are incompatible with SkyCoord's galcen_distance (should have galcen_distance^2 = ro^2 + zo^2)",galpyWarning) - # If at this point ro/vo not set, use default from config - if (_APY_LOADED and isinstance(vxvv,SkyCoord)) or radec or lb: - if ro is None: - ro= config.__config__.getfloat('normalization','ro') - if vo is None: - vo= config.__config__.getfloat('normalization','vo') - # If at this point zo not set, use default - if zo is None: zo= 0.025 - # if vxvv is SkyCoord, preferentially use its solarmotion - if _APY_LOADED and isinstance(vxvv,SkyCoord) \ - and not vxvv.galcen_v_sun is None: - sc_solarmotion= vxvv.galcen_v_sun.d_xyz.to(units.km/units.s).value - sc_solarmotion[0]= -sc_solarmotion[0] # right->left - sc_solarmotion[1]-= vo - if solarmotion is None: - solarmotion= sc_solarmotion - # If at this point solarmotion not set, use default - if solarmotion is None: solarmotion= 'schoenrich' - if isinstance(solarmotion,str) and solarmotion.lower() == 'hogg': - vsolar= nu.array([-10.1,4.0,6.7]) - elif isinstance(solarmotion,str) and solarmotion.lower() == 'dehnen': - vsolar= nu.array([-10.,5.25,7.17]) - elif isinstance(solarmotion,str) \ - and solarmotion.lower() == 'schoenrich': - vsolar= nu.array([-11.1,12.24,7.25]) - elif _APY_LOADED and isinstance(solarmotion,units.Quantity): - vsolar= solarmotion.to(units.km/units.s).value - else: - vsolar= nu.array(solarmotion) - # If both vxvv SkyCoord with vsun and solarmotion set, check the same - if _APY_LOADED and isinstance(vxvv,SkyCoord) \ - and not vxvv.galcen_v_sun is None: - if nu.any(nu.fabs(sc_solarmotion-vsolar) > 1e-8): - raise ValueError("Orbit initialization's solarmotion parameter not compatible with SkyCoord's galcen_v_sun; these should be the same for consistency (this may be because you did not set vo; galcen_v_sun = solarmotion+vo for consistency)") - # Now parse vxvv - if _APY_LOADED and isinstance(vxvv,SkyCoord): - galcen_v_sun= CartesianDifferential(\ - nu.array([-vsolar[0],vsolar[1]+vo,vsolar[2]])*units.km/units.s) - gc_frame= Galactocentric(\ - galcen_distance=nu.sqrt(ro**2.+zo**2.)*units.kpc, - z_sun=zo*units.kpc,galcen_v_sun=galcen_v_sun) - vxvvg= vxvv.transform_to(gc_frame) - if _APY_GE_31: - vxvvg.representation_type= 'cylindrical' - else: #pragma: no cover - vxvvg.representation= 'cylindrical' - R= vxvvg.rho.to(units.kpc).value/ro - phi= nu.pi-vxvvg.phi.to(units.rad).value - z= vxvvg.z.to(units.kpc).value/ro - try: - vR= vxvvg.d_rho.to(units.km/units.s).value/vo - except TypeError: - raise TypeError("SkyCoord given to Orbit initialization does not have velocity data, which is required to setup an Orbit") - vT= -(vxvvg.d_phi*vxvvg.rho)\ - .to(units.km/units.s, - equivalencies=units.dimensionless_angles()).value/vo - vz= vxvvg.d_z.to(units.km/units.s).value/vo - vxvv= [R,vR,vT,z,vz,phi] - elif radec or lb: - if radec: - if _APY_LOADED and isinstance(vxvv[0],units.Quantity): - ra, dec= vxvv[0].to(units.deg).value, \ - vxvv[1].to(units.deg).value - else: - ra, dec= vxvv[0], vxvv[1] - l,b= coords.radec_to_lb(ra,dec,degree=True,epoch=None) - _extra_rot= True - elif len(vxvv) == 4: - l, b= vxvv[0], 0. - _extra_rot= False - else: - l,b= vxvv[0],vxvv[1] - _extra_rot= True - if _APY_LOADED and isinstance(l,units.Quantity): - l= l.to(units.deg).value - if _APY_LOADED and isinstance(b,units.Quantity): - b= b.to(units.deg).value - if uvw: - if _APY_LOADED and isinstance(vxvv[2],units.Quantity): - X,Y,Z= coords.lbd_to_XYZ(l,b,vxvv[2].to(units.kpc).value, - degree=True) - else: - X,Y,Z= coords.lbd_to_XYZ(l,b,vxvv[2],degree=True) - vx= vxvv[3] - vy= vxvv[4] - vz= vxvv[5] - if _APY_LOADED and isinstance(vx,units.Quantity): - vx= vx.to(units.km/units.s).value - if _APY_LOADED and isinstance(vy,units.Quantity): - vy= vy.to(units.km/units.s).value - if _APY_LOADED and isinstance(vz,units.Quantity): - vz= vz.to(units.km/units.s).value - else: - if radec: - if _APY_LOADED and isinstance(vxvv[3],units.Quantity): - pmra, pmdec= vxvv[3].to(units.mas/units.yr).value, \ - vxvv[4].to(units.mas/units.yr).value - else: - pmra, pmdec= vxvv[3], vxvv[4] - pmll, pmbb= coords.pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec, - degree=True, - epoch=None) - d, vlos= vxvv[2], vxvv[5] - elif len(vxvv) == 4: - pmll, pmbb= vxvv[2], 0. - d, vlos= vxvv[1], vxvv[3] - else: - pmll, pmbb= vxvv[3], vxvv[4] - d, vlos= vxvv[2], vxvv[5] - if _APY_LOADED and isinstance(d,units.Quantity): - d= d.to(units.kpc).value - if _APY_LOADED and isinstance(vlos,units.Quantity): - vlos= vlos.to(units.km/units.s).value - if _APY_LOADED and isinstance(pmll,units.Quantity): - pmll= pmll.to(units.mas/units.yr).value - if _APY_LOADED and isinstance(pmbb,units.Quantity): - pmbb= pmbb.to(units.mas/units.yr).value - X,Y,Z,vx,vy,vz= coords.sphergal_to_rectgal(l,b,d, - vlos,pmll, pmbb, - degree=True) - X/= ro - Y/= ro - Z/= ro - vx/= vo - vy/= vo - vz/= vo - vsun= nu.array([0.,1.,0.,])+vsolar/vo - R, phi, z= coords.XYZ_to_galcencyl(X,Y,Z,Zsun=zo/ro, - _extra_rot=_extra_rot) - vR, vT,vz= coords.vxvyvz_to_galcencyl(vx,vy,vz, - R,phi,z, - vsun=vsun, - Xsun=1.,Zsun=zo/ro, - galcen=True, - _extra_rot=_extra_rot) - if lb and len(vxvv) == 4: vxvv= [R,vR,vT,phi] - else: vxvv= [R,vR,vT,z,vz,phi] - # Parse vxvv if it consists of Quantities - if _APY_LOADED and isinstance(vxvv[0],units.Quantity): - # Need to set ro and vo, default if not specified - if ro is None: - ro= config.__config__.getfloat('normalization','ro') - if vo is None: - vo= config.__config__.getfloat('normalization','vo') - new_vxvv= [vxvv[0].to(vxvv_units[0]).value/ro, - vxvv[1].to(vxvv_units[1]).value/vo] - if len(vxvv) > 2: - new_vxvv.append(vxvv[2].to(vxvv_units[2]).value/vo) - if len(vxvv) == 4: - new_vxvv.append(vxvv[3].to(vxvv_units[5]).value) - elif len(vxvv) > 4: - new_vxvv.append(vxvv[3].to(vxvv_units[3]).value/ro) - new_vxvv.append(vxvv[4].to(vxvv_units[4]).value/vo) - if len(vxvv) == 6: - new_vxvv.append(vxvv[5].to(vxvv_units[5]).value) - vxvv= new_vxvv - if len(vxvv) == 2: - self._orb= linearOrbit(vxvv=vxvv, - ro=ro,vo=vo) - elif len(vxvv) == 3: - self._orb= planarROrbit(vxvv=vxvv, - ro=ro,vo=vo,zo=zo,solarmotion=vsolar) - elif len(vxvv) == 4: - self._orb= planarOrbit(vxvv=vxvv, - ro=ro,vo=vo,zo=zo,solarmotion=vsolar) - elif len(vxvv) == 5: - self._orb= RZOrbit(vxvv=vxvv, - ro=ro,vo=vo,zo=zo,solarmotion=vsolar) - elif len(vxvv) == 6: - self._orb= FullOrbit(vxvv=vxvv, - ro=ro,vo=vo,zo=zo,solarmotion=vsolar) - #Store for actionAngle conversions - if vo is None: - self._vo= config.__config__.getfloat('normalization','vo') - self._voSet= False - else: - self._vo= vo - self._voSet= True - if ro is None: - self._ro= config.__config__.getfloat('normalization','ro') - self._roSet= False - else: - self._ro= ro - self._roSet= True - return None - - def setphi(self,phi): - """ - - NAME: - - setphi - - PURPOSE: - - set initial azimuth - - INPUT: - - phi - desired azimuth - - OUTPUT: - - (none) - - HISTORY: - - 2010-08-01 - Written - Bovy (NYU) - - BUGS: - - Should perform check that this orbit has phi - - """ - if len(self._orb.vxvv) == 2: - raise AttributeError("One-dimensional orbit has no azimuth") - elif len(self._orb.vxvv) == 3: - #Upgrade - vxvv= [self._orb.vxvv[0],self._orb.vxvv[1],self._orb.vxvv[2],phi] - self._orb= planarOrbit(vxvv=vxvv) - elif len(self._orb.vxvv) == 4: - self._orb.vxvv[-1]= phi - elif len(self._orb.vxvv) == 5: - #Upgrade - vxvv= [self._orb.vxvv[0],self._orb.vxvv[1],self._orb.vxvv[2], - self._orb.vxvv[3],self._orb.vxvv[4],phi] - self._orb= FullOrbit(vxvv=vxvv) - elif len(self._orb.vxvv) == 6: - self._orb.vxvv[-1]= phi - - def dim(self): - """ - NAME: - - dim - - PURPOSE: - - return the dimension of the Orbit - - INPUT: - - (none) - - OUTPUT: - - dimension - - HISTORY: - - 2011-02-03 - Written - Bovy (NYU) - - """ - if len(self._orb.vxvv) == 2: - return 1 - elif len(self._orb.vxvv) == 3 or len(self._orb.vxvv) == 4: - return 2 - elif len(self._orb.vxvv) == 5 or len(self._orb.vxvv) == 6: - return 3 - - def phasedim(self): - """ - NAME: - - phasedim - - PURPOSE: - - return the phase-space dimension of the problem (2 for 1D, 3 for 2D-axi, 4 for 2D, 5 for 3D-axi, 6 for 3D) - - INPUT: - - (none) - - OUTPUT: - - phase-space dimension - - HISTORY: - - 2018-12-20 - Written - Bovy (UofT) - - """ - return len(self._orb.vxvv) - - def turn_physical_off(self): - """ - NAME: - - turn_physical_off - - PURPOSE: - - turn off automatic returning of outputs in physical units - - INPUT: - - (none) - - OUTPUT: - - (none) - - HISTORY: - - 2014-06-17 - Written - Bovy (IAS) - - """ - self._roSet= False - self._voSet= False - self._orb.turn_physical_off() - - 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-01-19 - Written - Bovy (UofT) - - """ - self._roSet= True - self._voSet= True - if not ro is None: - if _APY_LOADED and isinstance(ro,units.Quantity): - ro= ro.to(units.kpc).value - self._ro= ro - if not vo is None: - if _APY_LOADED and isinstance(vo,units.Quantity): - vo= vo.to(units.km/units.s).value - self._vo= vo - self._orb.turn_physical_on(ro=ro,vo=vo) - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - - integrate - - PURPOSE: - - integrate the orbit - - INPUT: - - t - list of times at which to output (0 has to be in this!) (can be Quantity) - - pot - potential instance or list of instances - - method= 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'dop853' for a Dormand-Prince 8(5,3) implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'symplec4_c' for a 4th order symplectic integrator in C - 'symplec6_c' for a 6th order symplectic integrator in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince 5(4) integrator in C (generally the fastest) - 'dop853_c' for a Dormand-Prince 8(5, 3) integrator in C - - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize (only works for the C integrators that use a fixed stepsize) (can be Quantity) - - OUTPUT: - - (none) (get the actual orbit using getOrbit() - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - 2015-06-28 - Added dt keyword - Bovy (IAS) - - """ - if method.lower() not in ['odeint', 'leapfrog', 'dop853', 'leapfrog_c', - 'symplec4_c', 'symplec6_c', 'rk4_c', 'rk6_c', - 'dopr54_c', 'dop853_c']: - raise ValueError('{:s} is not a valid `method`'.format(method)) - pot= flatten_potential(pot) - _check_potential_dim(self,pot) - _check_consistent_units(self,pot) - # Parse t - if _APY_LOADED and isinstance(t,units.Quantity): - self._orb._integrate_t_asQuantity= True - t= t.to(units.Gyr).value\ - /bovy_conversion.time_in_Gyr(self._vo,self._ro) - if _APY_LOADED and not dt is None and isinstance(dt,units.Quantity): - dt= dt.to(units.Gyr).value\ - /bovy_conversion.time_in_Gyr(self._vo,self._ro) - from galpy.potential import MWPotential - if pot == MWPotential: - warnings.warn("Use of MWPotential as a Milky-Way-like potential is deprecated; galpy.potential.MWPotential2014, a potential fit to a large variety of dynamical constraints (see Bovy 2015), is the preferred Milky-Way-like potential in galpy", - galpyWarning) - if not _check_integrate_dt(t,dt): - raise ValueError('dt input (integrator stepsize) for Orbit.integrate must be an integer divisor of the output stepsize') - self._orb.integrate(t,pot,method=method,dt=dt) - - def integrate_dxdv(self,dxdv,t,pot,method='dopr54_c', - rectIn=False,rectOut=False): - """ - NAME: - - integrate_dxdv - - PURPOSE: - - integrate the orbit and a small area of phase space - - INPUT: - - dxdv - [dR,dvR,dvT,dphi] - - t - list of times at which to output (0 has to be in this!) (can be Quantity) - - pot - potential instance or list of instances - - method= 'odeint' for scipy's odeint - - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - - 'dopr54_c' is recommended, odeint is *not* recommended - - - rectIn= (False) if True, input dxdv is in rectangular coordinates - - rectOut= (False) if True, output dxdv (that in orbit_dxdv) is in rectangular coordinates - - OUTPUT: - - (none) (get the actual orbit using getOrbit_dxdv(), the orbit that is integrated alongside with dxdv is stored as usual, any previous regular orbit integration will be erased!) - - HISTORY: - - 2010-10-17 - Written - Bovy (IAS) - - 2014-06-29 - Added rectIn and rectOut - Bovy (IAS) - - """ - pot= flatten_potential(pot) - _check_potential_dim(self,pot) - _check_consistent_units(self,pot) - # Parse t - if _APY_LOADED and isinstance(t,units.Quantity): - self._orb._integrate_t_asQuantity= True - t= t.to(units.Gyr).value\ - /bovy_conversion.time_in_Gyr(self._vo,self._ro) - self._orb.integrate_dxdv(dxdv,t,pot,method=method, - rectIn=rectIn,rectOut=rectOut) - - def reverse(self): - """ - NAME: - - reverse - - PURPOSE: - - reverse an already integrated orbit (that is, make it go from end to beginning in t=0 to tend) - - INPUT: - - (none) - - OUTPUT: - - (none) - - HISTORY: - 2011-04-13 - Written - Bovy (NYU) - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - sortindx = list(range(len(self._orb.t))) - sortindx.sort(key=lambda x: self._orb.t[x],reverse=True) - for ii in range(self._orb.orbit.shape[1]): - self._orb.orbit[:,ii]= self._orb.orbit[sortindx,ii] - return None - - def flip(self,inplace=False): - """ - NAME: - - flip - - PURPOSE: - - 'flip' an orbit's initial conditions such that the velocities are minus the original velocities; useful for quick backward integration; returns a new Orbit instance - - INPUT: - - inplace= (False) if True, flip the orbit in-place, that is, without returning a new instance and also flip the velocities of the integrated orbit (if it exists) - - OUTPUT: - - Orbit instance that has the velocities of the current orbit flipped (inplace=False) or just flips all velocities of current instance (inplace=True) - - HISTORY: - - 2014-06-17 - Written - Bovy (IAS) - - 2016-07-21 - Added inplace keyword - Bovy (UofT) - - """ - if inplace: - self._orb.vxvv[1]= -self._orb.vxvv[1] - if len(self._orb.vxvv) > 2: - self._orb.vxvv[2]= -self._orb.vxvv[2] - if len(self._orb.vxvv) > 4: - self._orb.vxvv[4]= -self._orb.vxvv[4] - if hasattr(self._orb,'orbit'): - self._orb.orbit[:,1]= -self._orb.orbit[:,1] - if len(self._orb.vxvv) > 2: - self._orb.orbit[:,2]= -self._orb.orbit[:,2] - if len(self._orb.vxvv) > 4: - self._orb.orbit[:,4]= -self._orb.orbit[:,4] - if hasattr(self._orb,"_orbInterp"): - delattr(self._orb,"_orbInterp") - return None - orbSetupKwargs= {'ro':None, - 'vo':None, - 'zo':self._orb._zo, - 'solarmotion':self._orb._solarmotion} - if self._orb._roSet: - orbSetupKwargs['ro']= self._orb._ro - if self._orb._voSet: - orbSetupKwargs['vo']= self._orb._vo - if len(self._orb.vxvv) == 2: - return Orbit(vxvv= [self._orb.vxvv[0],-self._orb.vxvv[1]], - **orbSetupKwargs) - elif len(self._orb.vxvv) == 3: - return Orbit(vxvv=[self._orb.vxvv[0],-self._orb.vxvv[1], - -self._orb.vxvv[2]],**orbSetupKwargs) - elif len(self._orb.vxvv) == 4: - return Orbit(vxvv=[self._orb.vxvv[0],-self._orb.vxvv[1], - -self._orb.vxvv[2],self._orb.vxvv[3]], - **orbSetupKwargs) - elif len(self._orb.vxvv) == 5: - return Orbit(vxvv=[self._orb.vxvv[0],-self._orb.vxvv[1], - -self._orb.vxvv[2],self._orb.vxvv[3], - -self._orb.vxvv[4]],**orbSetupKwargs) - elif len(self._orb.vxvv) == 6: - return Orbit(vxvv= [self._orb.vxvv[0],-self._orb.vxvv[1], - -self._orb.vxvv[2],self._orb.vxvv[3], - -self._orb.vxvv[4],self._orb.vxvv[5]], - **orbSetupKwargs) - - def getOrbit(self): - """ - - NAME: - - getOrbit - - PURPOSE: - - return a previously calculated orbit - - INPUT: - - (none) - - OUTPUT: - - array orbit[nt,nd] - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.getOrbit() - - def getOrbit_dxdv(self): - """ - - NAME: - - getOrbit_dxdv - - PURPOSE: - - return a previously calculated integration of a small phase-space volume (with integrate_dxdv) - - INPUT: - - (none) - - OUTPUT: - - array orbit[nt,nd*2] with for each t the dxdv vector - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.getOrbit_dxdv() - - def fit(self,vxvv,vxvv_err=None,pot=None,radec=False,lb=False, - customsky=False,lb_to_customsky=None,pmllpmbb_to_customsky=None, - tintJ=10,ntintJ=1000,integrate_method='dopr54_c', - **kwargs): - """ - NAME: - - fit - - PURPOSE: - - fit an Orbit to data using the current orbit as the initial condition - - INPUT: - - vxvv - [:,6] array of positions and velocities along the orbit (if not lb=True or radec=True, these need to be in natural units [/ro,/vo], cannot be Quantities) - - vxvv_err= [:,6] array of errors on positions and velocities along the orbit (if None, these are set to 0.01) (if not lb=True or radec=True, these need to be in natural units [/ro,/vo], cannot be Quantities) - - pot= Potential to fit the orbit in - - Keywords related to the input data: - - radec= if True, input vxvv and vxvv are [ra,dec,d,mu_ra, mu_dec,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (all J2000.0; mu_ra = mu_ra * cos dec); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates - - lb= if True, input vxvv and vxvv are [long,lat,d,mu_ll, mu_bb,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (mu_ll = mu_ll * cos lat); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates - - customsky= if True, input vxvv and vxvv_err are [custom long,custom lat,d,mu_customll, mu_custombb,vlos] in [deg,deg,kpc,mas/yr,mas/yr,km/s] (mu_ll = mu_ll * cos lat) where custom longitude and custom latitude are a custom set of sky coordinates (e.g., ecliptic) and the proper motions are also expressed in these coordinats; you need to provide the functions lb_to_customsky and pmllpmbb_to_customsky to convert to the custom sky coordinates (these should have the same inputs and outputs as lb_to_radec and pmllpmbb_to_pmrapmdec); the attributes of the current Orbit are used to convert between these coordinates and Galactocentric coordinates - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s; entries can be Quantity) (default=Object-wide default) - Cannot be an Orbit instance with the orbit of the reference point, as w/ the ra etc. functions - Y is ignored and always assumed to be zero - - ro= distance in kpc corresponding to R=1. (default: taken from object; can be Quantity) - - vo= velocity in km/s corresponding to v=1. (default: taken from object; can be Quantity) - - lb_to_customsky= function that converts l,b,degree=False to the custom sky coordinates (like lb_to_radec); needs to be given when customsky=True - - pmllpmbb_to_customsky= function that converts pmll,pmbb,l,b,degree=False to proper motions in the custom sky coordinates (like pmllpmbb_to_pmrapmdec); needs to be given when customsky=True - - Keywords related to the orbit integrations: - - tintJ= (default: 10) time to integrate orbits for fitting the orbit (can be Quantity) - - ntintJ= (default: 1000) number of time-integration points - - integrate_method= (default: 'dopr54_c') integration method to use - - disp= (False) display the optimizer's convergence message - - OUTPUT: - - max of log likelihood - - HISTORY: - - 2014-06-17 - Written - Bovy (IAS) - - """ - pot= flatten_potential(pot) - _check_potential_dim(self,pot) - _check_consistent_units(self,pot) - return self._orb.fit(vxvv,vxvv_err=vxvv_err,pot=pot, - radec=radec,lb=lb, - customsky=customsky, - lb_to_customsky=lb_to_customsky, - pmllpmbb_to_customsky=pmllpmbb_to_customsky, - tintJ=tintJ,ntintJ=ntintJ, - integrate_method=integrate_method, - **kwargs) - - def E(self,*args,**kwargs): - """ - NAME: - - E - - PURPOSE: - - calculate the energy - - INPUT: - - t - (optional) time at which to get the energy (can be Quantity) - - pot= Potential instance or list of such instances - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - energy - - HISTORY: - - 2010-09-15 - Written - Bovy (NYU) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - _check_consistent_units(self,kwargs.get('pot',None)) - return self._orb.E(*args,**kwargs) - - def L(self,*args,**kwargs): - """ - NAME: - - L - - PURPOSE: - - calculate the angular momentum at time t - - INPUT: - - t - (optional) time at which to get the angular momentum (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - angular momentum - - HISTORY: - - 2010-09-15 - Written - Bovy (NYU) - - """ - return self._orb.L(*args,**kwargs) - - def Lz(self,*args,**kwargs): - """ - NAME: - - Lz - - PURPOSE: - - calculate the z-component of the angular momentum at time t - - INPUT: - - t - (optional) time at which to get the angular momentum (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - z-component of the angular momentum - - HISTORY: - - 2018-08-29 - Written - Bovy (UofT) - - """ - return self._orb.Lz(*args,**kwargs) - - def ER(self,*args,**kwargs): - """ - NAME: - - ER - - PURPOSE: - - calculate the radial energy - - INPUT: - - t - (optional) time at which to get the radial energy (can be Quantity) - - pot= Potential instance or list of such instances - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output (can be Quantity) - - OUTPUT: - - radial energy - - HISTORY: - - 2013-11-30 - Written - Bovy (IAS) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - _check_consistent_units(self,kwargs.get('pot',None)) - return self._orb.ER(*args,**kwargs) - - def Ez(self,*args,**kwargs): - """ - NAME: - - Ez - - PURPOSE: - - calculate the vertical energy - - INPUT: - - t - (optional) time at which to get the vertical energy (can be Quantity) - - pot= Potential instance or list of such instances - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output (can be Quantity) - - OUTPUT: - - vertical energy - - HISTORY: - - 2013-11-30 - Written - Bovy (IAS) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - _check_consistent_units(self,kwargs.get('pot',None)) - return self._orb.Ez(*args,**kwargs) - - def Jacobi(self,*args,**kwargs): - """ - NAME: - - Jacobi - - PURPOSE: - - calculate the Jacobi integral E - Omega L - - INPUT: - - t - (optional) time at which to get the Jacobi integral (can be Quantity) - - OmegaP= pattern speed (can be Quantity) - - pot= potential instance or list of such instances - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Jacobi integral - - HISTORY: - - 2011-04-18 - Written - Bovy (NYU) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - _check_consistent_units(self,kwargs.get('pot',None)) - out= self._orb.Jacobi(*args,**kwargs) - if not isinstance(out,float) and len(out) == 1: return out[0] - else: return out - - def e(self,analytic=False,pot=None,**kwargs): - """ - NAME: - - e - - PURPOSE: - - calculate the eccentricity, either numerically from the numerical orbit integration or using analytical means - - INPUT: - - analytic(= False) compute this analytically - - pot - potential to use for analytical calculation - - For 3D orbits different approximations for analytic=True are available (see the EccZmaxRperiRap method of actionAngle modules): - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic': assuming motion splits into R and z - - 2) 'staeckel': assuming motion splits into u and v of prolate spheroidal coordinate system, exact for Staeckel potentials (incl. all spherical potentials) - - 3) 'spherical': for spherical potentials, exact - - +actionAngle module setup kwargs for the corresponding actionAngle modules (actionAngleAdiabatic, actionAngleStaeckel, and actionAngleSpherical) - - OUTPUT: - - eccentricity - - HISTORY: - - 2010-09-15 - Written - Bovy (NYU) - - 2017-12-25 - Added Staeckel approximation and made that the default - Bovy (UofT) - - """ - - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - return self._orb.e(analytic=analytic,pot=pot,**kwargs) - - def rap(self,analytic=False,pot=None,**kwargs): - """ - NAME: - - rap - - PURPOSE: - - calculate the apocenter radius, either numerically from the numerical orbit integration or using analytical means - - INPUT: - - analytic(= False) compute this analytically - - pot - potential to use for analytical calculation - - For 3D orbits different approximations for analytic=True are available (see the EccZmaxRperiRap method of actionAngle modules): - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic': assuming motion splits into R and z - - 2) 'staeckel': assuming motion splits into u and v of prolate spheroidal coordinate system, exact for Staeckel potentials (incl. all spherical potentials) - - 3) 'spherical': for spherical potentials, exact - - +actionAngle module setup kwargs for the corresponding actionAngle modules (actionAngleAdiabatic, actionAngleStaeckel, and actionAngleSpherical) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - R_ap - - HISTORY: - - 2010-09-20 - Written - Bovy (NYU) - - 2017-12-25 - Added Staeckel approximation and made that the default - Bovy (UofT) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - return self._orb.rap(analytic=analytic,pot=pot,**kwargs) - - def rperi(self,analytic=False,pot=None,**kwargs): - """ - NAME: - - rperi - - PURPOSE: - - calculate the pericenter radius, either numerically from the numerical orbit integration or using analytical means - - INPUT: - - analytic(= False) compute this analytically - - pot - potential to use for analytical calculation - - For 3D orbits different approximations for analytic=True are available (see the EccZmaxRperiRap method of actionAngle modules): - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic': assuming motion splits into R and z - - 2) 'staeckel': assuming motion splits into u and v of prolate spheroidal coordinate system, exact for Staeckel potentials (incl. all spherical potentials) - - 3) 'spherical': for spherical potentials, exact - - +actionAngle module setup kwargs for the corresponding actionAngle modules (actionAngleAdiabatic, actionAngleStaeckel, and actionAngleSpherical) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - R_peri - - HISTORY: - - 2010-09-20 - Written - Bovy (NYU) - - 2017-12-25 - Added Staeckel approximation and made that the default - Bovy (UofT) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - return self._orb.rperi(analytic=analytic,pot=pot,**kwargs) - - @physical_conversion('position') - def rguiding(self,*args,**kwargs): - """ - NAME: - - rguiding - - PURPOSE: - - calculate the guiding-center radius (the radius of a circular orbit with the same angular momentum) - - INPUT: - - pot= potential instance or list of such instances - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - R_guiding - - HISTORY: - - 2018-08-29 - Written as thin wrapper around Potential.rl - Bovy (UofT) - - """ - pot= kwargs.get('pot',self._orb.__dict__.get('_pot',None)) - if pot is None: - raise RuntimeError("You need to specify the potential as pot= to compute the guiding-center radius") - flatten_potential(pot) - if _isNonAxi(pot): - raise RuntimeError('Potential given to rguiding is non-axisymmetric, but rguiding requires an axisymmetric potential') - _check_consistent_units(self,pot) - Lz= nu.atleast_1d(self.Lz(*args,use_physical=False)) - return nu.array([rl(pot,lz,use_physical=False) for lz in Lz]) - - def zmax(self,analytic=False,pot=None,**kwargs): - """ - NAME: - - zmax - - PURPOSE: - - calculate the maximum vertical height, either numerically from the numerical orbit integration or using analytical means - - INPUT: - - analytic(= False) compute this analytically - - pot - potential to use for analytical calculation - - For 3D orbits different approximations for analytic=True are available (see the EccZmaxRperiRap method of actionAngle modules): - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic': assuming motion splits into R and z - - 2) 'staeckel': assuming motion splits into u and v of prolate spheroidal coordinate system, exact for Staeckel potentials (incl. all spherical potentials) - - 3) 'spherical': for spherical potentials, exact - - +actionAngle module setup kwargs for the corresponding actionAngle modules (actionAngleAdiabatic, actionAngleStaeckel, and actionAngleSpherical) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Z_max - - HISTORY: - - 2010-09-20 - Written - Bovy (NYU) - - 2017-12-25 - Added Staeckel approximation and made that the default - Bovy (UofT) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - return self._orb.zmax(analytic=analytic,pot=pot,**kwargs) - - def resetaA(self,pot=None,type=None): - """ - NAME: - - resetaA - - PURPOSE: - - re-set up an actionAngle module for this Orbit - - INPUT: - - (none) - - OUTPUT: - - True if reset happened, False otherwise - - HISTORY: - - 2014-01-06 - Written - Bovy (IAS) - - """ - try: - delattr(self._orb,'_aA') - except AttributeError: - return False - else: - return True - - @physical_conversion('action') - def jr(self,pot=None,**kwargs): - """ - NAME: - - jr - - PURPOSE: - - calculate the radial action - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - jr - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA(self(),use_physical=False)[0]) - else: - return float(self._orb._aA(self,use_physical=False)[0]) - - @physical_conversion('action') - def jp(self,pot=None,**kwargs): - """ - NAME: - - jp - - PURPOSE: - - calculate the azimuthal action - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - jp - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA(self(),use_physical=False)[1]) - else: - return float(self._orb._aA(self,use_physical=False)[1]) - - @physical_conversion('action') - def jz(self,pot=None,**kwargs): - """ - NAME: - - jz - - PURPOSE: - - calculate the vertical action - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - jz - - HISTORY: - - 2012-06-01 - Written - Bovy (IAS) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA(self(),use_physical=False)[2]) - else: - return float(self._orb._aA(self,use_physical=False)[2]) - - @physical_conversion('angle') - def wr(self,pot=None,**kwargs): - """ - NAME: - - wr - - PURPOSE: - - calculate the radial angle - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - OUTPUT: - - wr - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqsAngles(self(), - use_physical=False)[6][0]) - else: - return float(self._orb._aA.actionsFreqsAngles(self, - use_physical=False)[6][0]) - - @physical_conversion('angle') - def wp(self,pot=None,**kwargs): - """ - NAME: - - wp - - PURPOSE: - - calculate the azimuthal angle - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - OUTPUT: - - wp - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqsAngles(self(), - use_physical=False)[7][0]) - else: - return float(self._orb._aA.actionsFreqsAngles(self, - use_physical=False)[7][0]) - - @physical_conversion('angle') - def wz(self,pot=None,**kwargs): - """ - NAME: - - wz - - PURPOSE: - - calculate the vertical angle - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - OUTPUT: - - wz - - HISTORY: - - 2012-06-01 - Written - Bovy (IAS) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqsAngles(self(), - use_physical=False)[8][0]) - else: - return float(self._orb._aA.actionsFreqsAngles(self, - use_physical=False)[8][0]) - - @physical_conversion('time') - def Tr(self,pot=None,**kwargs): - """ - NAME: - - Tr - - PURPOSE: - - calculate the radial period - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Tr - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self(), - use_physical=False)[3][0]) - else: - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self, - use_physical=False)[3][0]) - - @physical_conversion('time') - def Tp(self,pot=None,**kwargs): - """ - NAME: - - Tp - - PURPOSE: - - calculate the azimuthal period - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Tp - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self(), - use_physical=False)[4][0]) - else: - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self, - use_physical=False)[4][0]) - - def TrTp(self,pot=None,**kwargs): - """ - NAME: - - TrTp - - PURPOSE: - - the 'ratio' between the radial and azimuthal period Tr/Tphi*pi - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - OUTPUT: - - Tr/Tp*pi - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqs(self())[4][0]/self._orb._aA.actionsFreqs(self())[3][0]*nu.pi) - else: - return float(self._orb._aA.actionsFreqs(self)[4][0]/self._orb._aA.actionsFreqs(self)[3][0]*nu.pi) - - @physical_conversion('time') - def Tz(self,pot=None,**kwargs): - """ - NAME: - - Tz - - PURPOSE: - - calculate the vertical period - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Tz - - HISTORY: - - 2012-06-01 - Written - Bovy (IAS) - - 2013-11-27 - Re-written using new actionAngle modules - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self(), - use_physical=False)[5][0]) - else: - return float(2.*nu.pi/self._orb._aA.actionsFreqs(self, - use_physical=False)[5][0]) - - @physical_conversion('frequency') - def Or(self,pot=None,**kwargs): - """ - NAME: - - Or - - PURPOSE: - - calculate the radial frequency - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Or - - HISTORY: - - 2013-11-27 - Written - Bovy (IAS) - - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqs(self(),use_physical=False)[3][0]) - else: - return float(self._orb._aA.actionsFreqs(self,use_physical=False)[3][0]) - - @physical_conversion('frequency') - def Op(self,pot=None,**kwargs): - """ - NAME: - - Op - - PURPOSE: - - calculate the azimuthal frequency - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Op - - HISTORY: - - 2013-11-27 - Written - Bovy (IAS) - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqs(self(),use_physical=False)[4][0]) - else: - return float(self._orb._aA.actionsFreqs(self,use_physical=False)[4][0]) - - @physical_conversion('frequency') - def Oz(self,pot=None,**kwargs): - """ - NAME: - - Oz - - PURPOSE: - - calculate the vertical frequency - - INPUT: - - pot - potential - - type= ('staeckel') type of actionAngle module to use - - 1) 'adiabatic' - - 2) 'staeckel' - - 3) 'isochroneApprox' - - 4) 'spherical' - - +actionAngle module setup kwargs - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - Oz - - HISTORY: - - 2013-11-27 - Written - Bovy (IAS) - """ - if not pot is None: pot= flatten_potential(pot) - _check_consistent_units(self,pot) - self._orb._setupaA(pot=pot,**kwargs) - if self._orb._aAType.lower() == 'isochroneapprox': - return float(self._orb._aA.actionsFreqs(self(),use_physical=False)[5][0]) - else: - return float(self._orb._aA.actionsFreqs(self,use_physical=False)[5][0]) - - def time(self,*args,**kwargs): - """ - NAME: - - t - - PURPOSE: - - return the times at which the orbit is sampled - - INPUT: - - t - (default: integration times) time at which to get the time (for consistency reasons); default is to return the list of times at which the orbit is sampled (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - t(t) - - HISTORY: - - 2014-06-11 - Written - Bovy (IAS) - - """ - return self._orb.time(*args,**kwargs) - - def R(self,*args,**kwargs): - """ - NAME: - - R - - PURPOSE: - - return cylindrical radius at time t - - INPUT: - - t - (optional) time at which to get the radius (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - R(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.R(*args,**kwargs) - - def r(self,*args,**kwargs): - """ - NAME: - - r - - PURPOSE: - - return spherical radius at time t - - INPUT: - - t - (optional) time at which to get the radius (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - r(t) - - HISTORY: - - 2016-04-19 - Written - Bovy (UofT) - - """ - return self._orb.r(*args,**kwargs) - - def vR(self,*args,**kwargs): - """ - NAME: - - vR - - PURPOSE: - - return radial velocity at time t - - INPUT: - - t - (optional) time at which to get the radial velocity - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vR(t) - - HISTORY: - - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.vR(*args,**kwargs) - - def vT(self,*args,**kwargs): - """ - NAME: - - vT - - PURPOSE: - - return tangential velocity at time t - - INPUT: - - t - (optional) time at which to get the tangential velocity (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vT(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.vT(*args,**kwargs) - - def z(self,*args,**kwargs): - """ - NAME: - - z - - PURPOSE: - - return vertical height - - INPUT: - - t - (optional) time at which to get the vertical height (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - z(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.z(*args,**kwargs) - - def vz(self,*args,**kwargs): - """ - NAME: - - vz - - PURPOSE: - - return vertical velocity - - INPUT: - - t - (optional) time at which to get the vertical velocity (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vz(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.vz(*args,**kwargs) - - def phi(self,*args,**kwargs): - """ - NAME: - - phi - - PURPOSE: - - return azimuth - - INPUT: - - t - (optional) time at which to get the azimuth (can be Quantity) - - OUTPUT: - - phi(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.phi(*args,**kwargs) - - def vphi(self,*args,**kwargs): - """ - NAME: - - vphi - - PURPOSE: - - return angular velocity - - INPUT: - - t - (optional) time at which to get the angular velocity (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vphi(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - return self._orb.vphi(*args,**kwargs) - - def x(self,*args,**kwargs): - """ - NAME: - - x - - PURPOSE: - - return x - - INPUT: - - t - (optional) time at which to get x (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - x(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - out= self._orb.x(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def y(self,*args,**kwargs): - """ - NAME: - - y - - PURPOSE: - - return y - - INPUT: - - t - (optional) time at which to get y (can be Quantity) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - y(t) - - HISTORY: - - 2010-09-21 - Written - Bovy (NYU) - - """ - out= self._orb.y(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vx(self,*args,**kwargs): - """ - NAME: - - vx - - PURPOSE: - - return x velocity at time t - - INPUT: - - t - (optional) time at which to get the velocity (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vx(t) - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - """ - out= self._orb.vx(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vy(self,*args,**kwargs): - """ - - NAME: - - vy - - PURPOSE: - - return y velocity at time t - - INPUT: - - t - (optional) time at which to get the velocity (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - vy(t) - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - """ - out= self._orb.vy(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def ra(self,*args,**kwargs): - """ - NAME: - - ra - - PURPOSE: - - return the right ascension - - INPUT: - - t - (optional) time at which to get ra (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=[8.0,0.,0.]) OR Orbit object that corresponds to the orbit of the observer -(default=Object-wide default; can be Quantity) - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - ra(t) in deg - - HISTORY: - - 2011-02-23 - Written - Bovy (NYU) - - """ - out= self._orb.ra(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def dec(self,*args,**kwargs): - """ - NAME: - - dec - - PURPOSE: - - return the declination - - INPUT: - - t - (optional) time at which to get dec (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=[8.0,0.,0.]) OR Orbit object that corresponds to the orbit of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - dec(t) in deg - - HISTORY: - - 2011-02-23 - Written - Bovy (NYU) - - """ - out= self._orb.dec(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def ll(self,*args,**kwargs): - """ - NAME: - - ll - - PURPOSE: - - return Galactic longitude - - INPUT: - - t - (optional) time at which to get ll (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=[8.0,0.,0.]) OR Orbit object that corresponds to the orbit of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - l(t) in deg - - HISTORY: - - 2011-02-23 - Written - Bovy (NYU) - - """ - out= self._orb.ll(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def bb(self,*args,**kwargs): - """ - NAME: - - bb - - PURPOSE: - - return Galactic latitude - - INPUT: - - t - (optional) time at which to get bb (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=[8.0,0.,0.]) OR Orbit object that corresponds to the orbit of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - b(t) in deg - - HISTORY: - - 2011-02-23 - Written - Bovy (NYU) - - """ - out= self._orb.bb(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def dist(self,*args,**kwargs): - """ - NAME: - - dist - - PURPOSE: - - return distance from the observer - - INPUT: - - t - (optional) time at which to get dist (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=[8.0,0.,0.]) OR Orbit object that corresponds to the orbit of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - dist(t) in kpc - - HISTORY: - - 2011-02-23 - Written - Bovy (NYU) - - """ - out= self._orb.dist(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def pmra(self,*args,**kwargs): - """ - NAME: - - pmra - - PURPOSE: - - return proper motion in right ascension (in mas/yr) - - INPUT: - - t - (optional) time at which to get pmra (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantities) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - pm_ra(t) in mas/yr - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.pmra(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def pmdec(self,*args,**kwargs): - """ - NAME: - - pmdec - - PURPOSE: - - return proper motion in declination (in mas/yr) - - INPUT: - - t - (optional) time at which to get pmdec (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantities) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - pm_dec(t) in mas/yr - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.pmdec(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def pmll(self,*args,**kwargs): - """ - NAME: - - pmll - - PURPOSE: - - return proper motion in Galactic longitude (in mas/yr) - - INPUT: - - t - (optional) time at which to get pmll (can be Quantity) - -v obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantities) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - pm_l(t) in mas/yr - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.pmll(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def pmbb(self,*args,**kwargs): - """ - NAME: - - pmbb - - PURPOSE: - - return proper motion in Galactic latitude (in mas/yr) - - INPUT: - - t - (optional) time at which to get pmbb (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - pm_b(t) in mas/yr - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.pmbb(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vlos(self,*args,**kwargs): - """ - NAME: - - vlos - - PURPOSE: - - return the line-of-sight velocity (in km/s) - - INPUT: - - t - (optional) time at which to get vlos (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - vlos(t) in km/s - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.vlos(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vra(self,*args,**kwargs): - """ - NAME: - - vra - - PURPOSE: - - return velocity in right ascension (km/s) - - INPUT: - - t - (optional) time at which to get vra (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - v_ra(t) in km/s - - HISTORY: - - 2011-03-27 - Written - Bovy (NYU) - - """ - from .OrbitTop import _check_roSet, _check_voSet - _check_roSet(self,kwargs,'vra') - _check_voSet(self,kwargs,'vra') - dist= self._orb.dist(*args,**kwargs) - if _APY_UNITS and isinstance(dist,units.Quantity): - out= units.Quantity(dist.to(units.kpc).value*_K* - self._orb.pmra(*args,**kwargs)\ - .to(units.mas/units.yr).value, - unit=units.km/units.s) - else: - out= dist*_K*self._orb.pmra(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vdec(self,*args,**kwargs): - """ - NAME: - - vdec - - PURPOSE: - - return velocity in declination (km/s) - - INPUT: - - t - (optional) time at which to get vdec (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - v_dec(t) in km/s - - HISTORY: - - 2011-03-27 - Written - Bovy (NYU) - - """ - from .OrbitTop import _check_roSet, _check_voSet - _check_roSet(self,kwargs,'vdec') - _check_voSet(self,kwargs,'vdec') - dist= self._orb.dist(*args,**kwargs) - if _APY_UNITS and isinstance(dist,units.Quantity): - out= units.Quantity(dist.to(units.kpc).value*_K* - self._orb.pmdec(*args,**kwargs)\ - .to(units.mas/units.yr).value, - unit=units.km/units.s) - else: - out= dist*_K*self._orb.pmdec(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vll(self,*args,**kwargs): - """ - NAME: - - vll - - PURPOSE: - - return the velocity in Galactic longitude (km/s) - - INPUT: - - t - (optional) time at which to get vll (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - v_l(t) in km/s - - HISTORY: - - 2011-03-27 - Written - Bovy (NYU) - - """ - from .OrbitTop import _check_roSet, _check_voSet - _check_roSet(self,kwargs,'vll') - _check_voSet(self,kwargs,'vll') - dist= self._orb.dist(*args,**kwargs) - if _APY_UNITS and isinstance(dist,units.Quantity): - out= units.Quantity(dist.to(units.kpc).value*_K* - self._orb.pmll(*args,**kwargs)\ - .to(units.mas/units.yr).value, - unit=units.km/units.s) - else: - out= dist*_K*self._orb.pmll(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def vbb(self,*args,**kwargs): - """ - NAME: - - vbb - - PURPOSE: - - return velocity in Galactic latitude (km/s) - - INPUT: - - t - (optional) time at which to get vbb (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - v_b(t) in km/s - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - from .OrbitTop import _check_roSet, _check_voSet - _check_roSet(self,kwargs,'vbb') - _check_voSet(self,kwargs,'vbb') - dist= self._orb.dist(*args,**kwargs) - if _APY_UNITS and isinstance(dist,units.Quantity): - out= units.Quantity(dist.to(units.kpc).value*_K* - self._orb.pmbb(*args,**kwargs)\ - .to(units.mas/units.yr).value, - unit=units.km/units.s) - else: - out= dist*_K*self._orb.pmbb(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def helioX(self,*args,**kwargs): - """ - NAME: - - helioX - - PURPOSE: - - return Heliocentric Galactic rectangular x-coordinate (aka "X") - - INPUT: - - t - (optional) time at which to get X (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - helioX(t) in kpc - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.helioX(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def helioY(self,*args,**kwargs): - """ - NAME: - - helioY - - PURPOSE: - - return Heliocentric Galactic rectangular y-coordinate (aka "Y") - - INPUT: - - t - (optional) time at which to get Y (can be Quantity) - - obs=[X,Y,Z] - (optional) position and of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.]; entries can be Quantity)) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - helioY(t) in kpc - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.helioY(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def helioZ(self,*args,**kwargs): - """ - NAME: - - helioZ - - PURPOSE: - - return Heliocentric Galactic rectangular z-coordinate (aka "Z") - - INPUT: - - t - (optional) time at which to get Z (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - OUTPUT: - - helioZ(t) in kpc - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.helioZ(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def U(self,*args,**kwargs): - """ - NAME: - - U - - PURPOSE: - - return Heliocentric Galactic rectangular x-velocity (aka "U") - - INPUT: - - t - (optional) time at which to get U (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - U(t) in km/s - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.U(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def V(self,*args,**kwargs): - """ - NAME: - - V - - PURPOSE: - - return Heliocentric Galactic rectangular y-velocity (aka "V") - - INPUT: - - t - (optional) time at which to get V (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - V(t) in km/s - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.V(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def W(self,*args,**kwargs): - """ - NAME: - - W - - PURPOSE: - - return Heliocentric Galactic rectangular z-velocity (aka "W") - - INPUT: - - t - (optional) time at which to get W (can be Quantity) - - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - in the Galactocentric frame - (in kpc and km/s) (default=[8.0,0.,0.,0.,220.,0.]; entries can be Quantity) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - W(t) in km/s - - HISTORY: - - 2011-02-24 - Written - Bovy (NYU) - - """ - out= self._orb.W(*args,**kwargs) - if len(out) == 1: return out[0] - else: return out - - def SkyCoord(self,*args,**kwargs): - """ - NAME: - - SkyCoord - - PURPOSE: - - return the position and velocity as an astropy SkyCoord - - INPUT: - - t - (optional) time at which to get the position (can be Quantity) - - obs=[X,Y,Z] - (optional) position of observer (in kpc; entries can be Quantity) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - OUTPUT: - - SkyCoord(t) - - HISTORY: - - 2015-06-02 - Written - Bovy (IAS) - - 2018-04-19 - Added velocity output for astropy >= v3 - Bovy (UofT) - - """ - return self._orb.SkyCoord(*args,**kwargs) - - def __call__(self,*args,**kwargs): - """ - NAME: - - __call__ - - PURPOSE: - - return the orbit at time t - - INPUT: - - t - desired time (can be Quantity) - - OUTPUT: - - an Orbit instance with initial condition set to the - phase-space at time t or list of Orbit instances if multiple - times are given - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - orbSetupKwargs= {'ro':None, - 'vo':None, - 'zo':self._orb._zo, - 'solarmotion':self._orb._solarmotion} - if self._orb._roSet: - orbSetupKwargs['ro']= self._orb._ro - if self._orb._voSet: - orbSetupKwargs['vo']= self._orb._vo - thiso= self._orb(*args,**kwargs) - if len(thiso.shape) == 1: return Orbit(vxvv=thiso,**orbSetupKwargs) - else: return [Orbit(vxvv=thiso[:,ii], - **orbSetupKwargs) for ii in range(thiso.shape[1])] - - def plot(self,*args,**kwargs): - """ - NAME: - - plot - - PURPOSE: - - plot a previously calculated orbit (with reasonable defaults) - - INPUT: - - d1= first dimension to plot ('x', 'y', 'R', 'vR', 'vT', 'z', 'vz', ...); can also be a user-defined function of time (e.g., lambda t: o.R(t) for R) - - d2= second dimension to plot; can also be a user-defined function of time (e.g., lambda t: o.R(t) for R) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - matplotlib.plot inputs+bovy_plot.plot inputs - - OUTPUT: - - sends plot to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plot(*args,**kwargs) - - def plot3d(self,*args,**kwargs): - """ - NAME: - - plot3d - - PURPOSE: - - plot 3D aspects of an Orbit - - INPUT: - - d1= first dimension to plot ('x', 'y', 'R', 'vR', 'vT', 'z', 'vz', ...); can also be a user-defined function of time (e.g., lambda t: o.R(t) for R) - - d2= second dimension to plot - - d3= third dimension to plot - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot3d args and kwargs - - OUTPUT: - - plot - - HISTORY: - - 2010-07-26 - Written - Bovy (NYU) - - 2010-09-22 - Adapted to more general framework - Bovy (NYU) - - 2010-01-08 - Adapted to 3D - Bovy (NYU) - """ - return self._orb.plot3d(*args,**kwargs) - - def plotE(self,*args,**kwargs): - """ - NAME: - - plotE - - PURPOSE: - - plot E(.) along the orbit - - INPUT: - - pot= Potential instance or list of instances in which the orbit was integrated - - d1= plot Ez vs d1: e.g., 't', 'z', 'R', 'vR', 'vT', 'vz' - - normed= if set, plot E(t)/E(0) rather than E(t) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - +bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - 2014-06-16 - Changed to actually plot E rather than E/E0 - Bovy (IAS) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - return self._orb.plotE(*args,**kwargs) - - def plotEz(self,*args,**kwargs): - """ - NAME: - - plotEz - - PURPOSE: - - plot E_z(.) along the orbit - - INPUT: - - pot= Potential instance or list of instances in which the orbit was integrated - - d1= plot Ez vs d1: e.g., 't', 'z', 'R' - - normed= if set, plot E(t)/E(0) rather than E(t) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - +bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - return self._orb.plotEz(*args,**kwargs) - - def plotER(self,*args,**kwargs): - """ - NAME: - - plotER - - PURPOSE: - - plot E_R(.) along the orbit - - INPUT: - - pot= Potential instance or list of instances in which the orbit was integrated - - d1= plot ER vs d1: e.g., 't', 'z', 'R' - - normed= if set, plot E(t)/E(0) rather than E(t) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - +bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - return self._orb.plotER(*args,**kwargs) - - def plotEzJz(self,*args,**kwargs): - """ - NAME: - - plotEzJzt - - PURPOSE: - - plot E_z(t)/sqrt(dens(R)) / (E_z(0)/sqrt(dens(R(0)))) along the orbit (an approximation to the vertical action) - - INPUT: - - pot - Potential instance or list of instances in which the orbit was integrated - - d1= plot Ez vs d1: e.g., 't', 'z', 'R' - - +bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-08-08 - Written - Bovy (NYU) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - return self._orb.plotEzJz(*args,**kwargs) - - def plotJacobi(self,*args,**kwargs): - """ - NAME: - - plotJacobi - - PURPOSE: - - plot the Jacobi integral along the orbit - - INPUT: - - OmegaP= pattern speed - - pot= - Potential instance or list of instances in which the orbit - was integrated - - d1= - plot Ez vs d1: e.g., 't', 'z', 'R', 'vR', 'vT', 'vz' - - normed= if set, plot E(t)/E(0) rather than E(t) - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - +bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2011-10-10 - Written - Bovy (IAS) - - """ - if not kwargs.get('pot',None) is None: kwargs['pot']= flatten_potential(kwargs.get('pot')) - return self._orb.plotJacobi(*args,**kwargs) - - def plotR(self,*args,**kwargs): - """ - NAME: - - plotR - - PURPOSE: - - plot R(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotR(*args,**kwargs) - - def plotz(self,*args,**kwargs): - """ - NAME: - - plotz - - PURPOSE: - - plot z(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotz(*args,**kwargs) - - def plotvR(self,*args,**kwargs): - """ - NAME: - - plotvR - - PURPOSE: - - plot vR(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotvR(*args,**kwargs) - - def plotvT(self,*args,**kwargs): - """ - NAME: - - plotvT - - PURPOSE: - - plot vT(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotvT(*args,**kwargs) - - def plotphi(self,*args,**kwargs): - """ - NAME: - - plotphi - - PURPOSE: - - plot \phi(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotphi(*args,**kwargs) - - def plotvz(self,*args,**kwargs): - """ - NAME: - - plotvz - - PURPOSE: - - plot vz(.) along the orbit - - INPUT: - d1= plot vs d1: e.g., 't', 'z', 'R' - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - return self._orb.plotvz(*args,**kwargs) - - def plotx(self,*args,**kwargs): - """ - NAME: - - plotx - - PURPOSE: - - plot x(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-21 - Written - Bovy (NYU) - - """ - return self._orb.plotx(*args,**kwargs) - - def plotvx(self,*args,**kwargs): - """ - NAME: - - plotvx - - PURPOSE: - - plot vx(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-21 - Written - Bovy (NYU) - - """ - return self._orb.plotvx(*args,**kwargs) - - def ploty(self,*args,**kwargs): - """ - NAME: - - ploty - - PURPOSE: - - plot y(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-21 - Written - Bovy (NYU) - - """ - return self._orb.ploty(*args,**kwargs) - - def plotvy(self,*args,**kwargs): - """ - NAME: - - plotvy - - PURPOSE: - - plot vy(.) along the orbit - - INPUT: - - d1= plot vs d1: e.g., 't', 'z', 'R' - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - bovy_plot.bovy_plot inputs - - OUTPUT: - - figure to output device - - HISTORY: - - 2010-07-21 - Written - Bovy (NYU) - - """ - return self._orb.plotvy(*args,**kwargs) - - def toPlanar(self): - """ - NAME: - - toPlanar - - PURPOSE: - - convert a 3D orbit into a 2D orbit - - INPUT: - - (none) - - OUTPUT: - - planar Orbit - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - """ - orbSetupKwargs= {'ro':None, - 'vo':None, - 'zo':self._orb._zo, - 'solarmotion':self._orb._solarmotion} - if self._orb._roSet: - orbSetupKwargs['ro']= self._orb._ro - if self._orb._voSet: - orbSetupKwargs['vo']= self._orb._vo - if len(self._orb.vxvv) == 6: - vxvv= [self._orb.vxvv[0],self._orb.vxvv[1], - self._orb.vxvv[2],self._orb.vxvv[5]] - elif len(self._orb.vxvv) == 5: - vxvv= [self._orb.vxvv[0],self._orb.vxvv[1],self._orb.vxvv[2]] - else: - raise AttributeError("planar or linear Orbits do not have the toPlanar attribute") - return Orbit(vxvv=vxvv,**orbSetupKwargs) - - def toLinear(self): - """ - NAME: - - toLinear - - PURPOSE: - - convert a 3D orbit into a 1D orbit (z) - - INPUT: - - (none) - - OUTPUT: - - linear Orbit - - HISTORY: - - 2010-11-30 - Written - Bovy (NYU) - - """ - orbSetupKwargs= {'ro':None, - 'vo':None, - 'zo':self._orb._zo, - 'solarmotion':self._orb._solarmotion} - if self._orb._roSet: - orbSetupKwargs['ro']= self._orb._ro - if self._orb._voSet: - orbSetupKwargs['vo']= self._orb._vo - if len(self._orb.vxvv) == 6 or len(self._orb.vxvv) == 5: - vxvv= [self._orb.vxvv[3],self._orb.vxvv[4]] - else: - raise AttributeError("planar or linear Orbits do not have the toPlanar attribute") - return Orbit(vxvv=vxvv,**orbSetupKwargs) - - def __add__(self,linOrb): - """ - NAME: - - __add__ - - PURPOSE: - - add a linear orbit and a planar orbit to make a 3D orbit - - INPUT: - - linear or plane orbit instance - - OUTPUT: - - 3D orbit - - HISTORY: - - 2010-07-21 - Written - Bovy (NYU) - - """ - orbSetupKwargs= {'ro':None, - 'vo':None, - 'zo':self._orb._zo, - 'solarmotion':self._orb._solarmotion} - if self._orb._roSet: - orbSetupKwargs['ro']= self._orb._ro - if self._orb._voSet: - orbSetupKwargs['vo']= self._orb._vo - if (not (isinstance(self._orb,planarOrbitTop) and - isinstance(linOrb._orb,linearOrbit)) and - not (isinstance(self._orb,linearOrbit) and - isinstance(linOrb._orb,planarOrbitTop))): - raise AttributeError("Only planarOrbit+linearOrbit is supported") - if isinstance(self._orb,planarROrbit): - return Orbit(vxvv=[self._orb.vxvv[0],self._orb.vxvv[1], - self._orb.vxvv[2], - linOrb._orb.vxvv[0],linOrb._orb.vxvv[1]], - **orbSetupKwargs) - elif isinstance(self._orb,planarOrbit): - return Orbit(vxvv=[self._orb.vxvv[0],self._orb.vxvv[1], - self._orb.vxvv[2], - linOrb._orb.vxvv[0],linOrb._orb.vxvv[1], - self._orb.vxvv[3]], - **orbSetupKwargs) - elif isinstance(linOrb._orb,planarROrbit): - return Orbit(vxvv=[linOrb._orb.vxvv[0],linOrb._orb.vxvv[1], - linOrb._orb.vxvv[2], - self._orb.vxvv[0],self._orb.vxvv[1]], - **orbSetupKwargs) - elif isinstance(linOrb._orb,planarOrbit): - return Orbit(vxvv=[linOrb._orb.vxvv[0],linOrb._orb.vxvv[1], - linOrb._orb.vxvv[2], - self._orb.vxvv[0],self._orb.vxvv[1], - linOrb._orb.vxvv[3]], - **orbSetupKwargs) - - def animate(self,*args,**kwargs): #pragma: no cover - """ - NAME: - - animate - - PURPOSE: - - animate a previously calculated orbit (with reasonable defaults) - - INPUT: - - d1= first dimension to plot ('x', 'y', 'R', 'vR', 'vT', 'z', 'vz', ...); can be list with up to three entries for three subplots; each entry can also be a user-defined function of time (e.g., lambda t: o.R(t) for R) - - d2= second dimension to plot; can be list with up to three entries for three subplots; each entry can also be a user-defined function of time (e.g., lambda t: o.R(t) for R) - - width= (600) width of output div in px - - height= (400) height of output div in px - - xlabel= (pre-defined labels) label for the first dimension (or list of labels if d1 is a list); should only have to be specified when using a function as d1 and can then specify as, e.g., [None,'YOUR LABEL',None] if d1 is a list of three xs and the first and last are standard entries) - - ylabel= (pre-defined labels) label for the second dimension (or list of labels if d2 is a list); should only have to be specified when using a function as d2 and can then specify as, e.g., [None,'YOUR LABEL',None] if d1 is a list of three xs and the first and last are standard entries) - - json_filename= (None) if set, save the data necessary for the figure in this filename (e.g., json_filename= 'orbit_data/orbit.json'); this path is also used in the output HTML, so needs to be accessible - - ro= (Object-wide default) physical scale for distances to use to convert (can be Quantity) - - vo= (Object-wide default) physical scale for velocities to use to convert (can be Quantity) - - use_physical= use to override Object-wide default for using a physical scale for output - - OUTPUT: - - IPython.display.HTML object with code to animate the orbit; can be directly shown in jupyter notebook or embedded in HTML pages; get a text version of the HTML using the _repr_html_() function - - HISTORY: - - 2017-09-17-24 - Written - Bovy (UofT) - - """ - return self._orb.animate(*args,**kwargs) - - @classmethod - def from_name(cls, name, vo=None, ro=None, zo=None, solarmotion=None): - """ - NAME: - - from_name - - PURPOSE: - - given the name of an object, retrieve coordinate information for that object from SIMBAD and return a corresponding orbit - - INPUT: - - name - the name of the object - - +standard Orbit initialization keywords: - - ro= distance from vantage point to GC (kpc; can be Quantity) - - vo= circular velocity at ro (km/s; can be Quantity) - - zo= offset toward the NGP of the Sun wrt the plane (kpc; can be Quantity; default = 25 pc) - - solarmotion= 'hogg' or 'dehnen', or 'schoenrich', or value in [-U,V,W]; can be Quantity - - OUTPUT: - - orbit containing the phase space coordinates of the named object - - HISTORY: - - 2018-07-15 - Written - Mathew Bub (UofT) - - """ - if not _APY_LOADED: # pragma: no cover - raise ImportError('astropy needs to be installed to use ' - 'Orbit.from_name') - if not _ASTROQUERY_LOADED: # pragma: no cover - raise ImportError('astroquery needs to be installed to use ' - 'Orbit.from_name') - - # setup a SIMBAD query with the appropriate fields - simbad= Simbad() - simbad.add_votable_fields('ra(d)', 'dec(d)', 'pmra', 'pmdec', - 'rv_value', 'plx', 'distance') - simbad.remove_votable_fields('main_id', 'coordinates') - - # query SIMBAD for the named object - try: - simbad_table= simbad.query_object(name) - except OSError: # pragma: no cover - raise OSError('failed to connect to SIMBAD') - if not simbad_table: - raise ValueError('failed to find {} in SIMBAD'.format(name)) - - # check that the necessary coordinates have been found - missing= simbad_table.mask - if (any(missing['RA_d', 'DEC_d', 'PMRA', 'PMDEC', 'RV_VALUE'][0]) or - all(missing['PLX_VALUE', 'Distance_distance'][0])): - raise ValueError('failed to find some coordinates for {} in ' - 'SIMBAD'.format(name)) - ra, dec, pmra, pmdec, vlos= simbad_table['RA_d', 'DEC_d', 'PMRA', - 'PMDEC', 'RV_VALUE'][0] - - # get a distance value - if not missing['PLX_VALUE'][0]: - dist= 1/simbad_table['PLX_VALUE'][0] - else: - dist_str= str(simbad_table['Distance_distance'][0]) + \ - simbad_table['Distance_unit'][0] - dist= units.Quantity(dist_str).to(units.kpc).value - - return cls(vxvv=[ra,dec,dist,pmra,pmdec,vlos], radec=True, ro=ro, vo=vo, - zo=zo, solarmotion=solarmotion) diff --git a/galpy/orbit/OrbitTop.py b/galpy/orbit/OrbitTop.py deleted file mode 100644 index 8acf0fecb..000000000 --- a/galpy/orbit/OrbitTop.py +++ /dev/null @@ -1,2207 +0,0 @@ -import warnings -import json -from random import choice -from string import ascii_lowercase -import math as m -import numpy as nu -import scipy -from scipy import interpolate -_APY_LOADED= True -try: - from astropy import units, coordinates -except ImportError: - _APY_LOADED= False -if _APY_LOADED: - import astropy - _APY3= astropy.__version__ > '3' -from .Orbits import _check_roSet, _check_voSet, _helioXYZ, _lbd, _radec, \ - _XYZvxvyvz, _lbdvrpmllpmbb, _pmrapmdec -from galpy import actionAngle -from galpy.potential import PotentialError -import galpy.util.bovy_plot as plot -import galpy.util.bovy_coords as coords -from galpy.util.bovy_conversion import physical_conversion -from galpy.util import bovy_conversion, galpyWarning -from galpy.util import config -if int(scipy.__version__.split('.')[0]) < 1 and \ - int(scipy.__version__.split('.')[1]) < 15: #pragma: no cover - _OLD_SCIPY= True - _KWINTERP= {} #for scipy version <1.15 -else: - _OLD_SCIPY= False - _KWINTERP= {'ext':2} #for scipy version >=1.15 -class OrbitTop(object): - """General class that holds orbits and integrates them""" - def __init__(self,vxvv=None,vo=None,ro=None,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): - """ - NAME: - - __init__ - - PURPOSE: - - Initialize an orbit instance - - INPUT: - - vxvv - initial condition - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - (none) - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - """ - # If you change the way an Orbit object is setup, also change each of - # the methods that return Orbits - self.vxvv= vxvv - if vo is None: - self._vo= config.__config__.getfloat('normalization','vo') - self._voSet= False - else: - self._vo= vo - self._voSet= True - if ro is None: - self._ro= config.__config__.getfloat('normalization','ro') - self._roSet= False - else: - self._ro= ro - self._roSet= True - self._zo= zo - self._solarmotion= solarmotion - return None - - def turn_physical_off(self): - """ - NAME: - turn_physical_off - PURPOSE: - turn off automatic returning of outputs in physical units - INPUT: - (none) - OUTPUT: - (none) - HISTORY: - 2014-06-17 - Written - Bovy (IAS) - """ - self._roSet= False - self._voSet= False - return None - - 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) - vo= reference velocity (km/s) - OUTPUT: - (none) - HISTORY: - 2016-01-19 - Written - Bovy (UofT) - """ - self._roSet= True - self._voSet= True - if not ro is None: - self._ro= ro - if not vo is None: - self._vo= vo - return None - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - Potential instance or list of instances - OUTPUT: - (none) (get the actual orbit using self.getOrbit() - HISTORY: - 2010-07-10 - """ - raise NotImplementedError - - def getOrbit(self): - """ - NAME: - getOrbit - PURPOSE: - return a previously calculated orbit - INPUT: - (none) - OUTPUT: - (none) - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - return self.orbit - - def getOrbit_dxdv(self): - """ - NAME: - getOrbit_dxdv - PURPOSE: - return a previously calculated orbit_dxdv - INPUT: - (none) - OUTPUT: - (none) - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - return self.orbit_dxdv[:,4:] - - @physical_conversion('time') - def time(self,*args,**kwargs): - """ - NAME: - time - PURPOSE: - return the times at which the orbit is sampled - INPUT: - t - (default: integration times) time at which to get the time (for consistency reasons); default is to return the list of times at which the orbit is sampled - ro= (Object-wide default) physical scale for distances to use to convert - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - t(t) - HISTORY: - 2014-06-11 - Written - Bovy (IAS) - """ - if len(args) == 0: - try: - return self.t - except AttributeError: - return 0. - else: return args[0] - - @physical_conversion('position') - def R(self,*args,**kwargs): - """ - NAME: - R - PURPOSE: - return cylindrical radius at time t - INPUT: - t - (optional) time at which to get the radius - ro= (Object-wide default) physical scale for distances to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - R(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[0] - else: return thiso[0,:] - - @physical_conversion('position') - def r(self,*args,**kwargs): - """ - NAME: - r - PURPOSE: - return spherical radius at time t - INPUT: - t - (optional) time at which to get the radius - ro= (Object-wide default) physical scale for distances to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - r(t) - HISTORY: - 2016-04-19 - Written - Bovy (UofT) - """ - thiso= self(*args,**kwargs) - if len(thiso) > 4: - return nu.sqrt(thiso[0]**2.+thiso[3]**2.) - else: - return nu.fabs(thiso[0]) - - @physical_conversion('velocity') - def vR(self,*args,**kwargs): - """ - NAME: - vR - PURPOSE: - return radial velocity at time t - INPUT: - t - (optional) time at which to get the radial velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vR(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[1] - else: return thiso[1,:] - - @physical_conversion('velocity') - def vT(self,*args,**kwargs): - """ - NAME: - vT - PURPOSE: - return tangential velocity at time t - INPUT: - t - (optional) time at which to get the tangential velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vT(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[2] - else: return thiso[2,:] - - @physical_conversion('position') - def z(self,*args,**kwargs): - """ - NAME: - z - PURPOSE: - return vertical height - INPUT: - t - (optional) time at which to get the vertical height - ro= (Object-wide default) physical scale for distances to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - z(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - if len(self.vxvv) < 5: - raise AttributeError("linear and planar orbits do not have z()") - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[3] - else: return thiso[3,:] - - @physical_conversion('velocity') - def vz(self,*args,**kwargs): - """ - NAME: - vz - PURPOSE: - return vertical velocity - INPUT: - t - (optional) time at which to get the vertical velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vz(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - if len(self.vxvv) < 5: - raise AttributeError("linear and planar orbits do not have vz()") - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[4] - else: return thiso[4,:] - - @physical_conversion('angle') - def phi(self,*args,**kwargs): - """ - NAME: - phi - PURPOSE: - return azimuth - INPUT: - t - (optional) time at which to get the azimuth - OUTPUT: - phi(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - if len(self.vxvv) != 4 and len(self.vxvv) != 6: - raise AttributeError("orbit must track azimuth to use phi()") - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: return thiso[-1] - else: return thiso[-1,:] - - @physical_conversion('position') - def x(self,*args,**kwargs): - """ - NAME: - x - PURPOSE: - return x - INPUT: - t - (optional) time at which to get x - ro= (Object-wide default) physical scale for distances to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - x(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - if len(thiso[:,0]) == 2: - return thiso[0,:] - if len(thiso[:,0]) != 4 and len(thiso[:,0]) != 6: - raise AttributeError("orbit must track azimuth to use x()") - elif len(thiso[:,0]) == 4: - return thiso[0,:]*nu.cos(thiso[3,:]) - else: - return thiso[0,:]*nu.cos(thiso[5,:]) - - @physical_conversion('position') - def y(self,*args,**kwargs): - """ - NAME: - y - PURPOSE: - return y - INPUT: - t - (optional) time at which to get y - ro= (Object-wide default) physical scale for distances to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - y(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - if len(thiso[:,0]) != 4 and len(thiso[:,0]) != 6: - raise AttributeError("orbit must track azimuth to use x()") - elif len(thiso[:,0]) == 4: - return thiso[0,:]*nu.sin(thiso[3,:]) - else: - return thiso[0,:]*nu.sin(thiso[5,:]) - - @physical_conversion('velocity') - def vx(self,*args,**kwargs): - """ - NAME: - vx - PURPOSE: - return x velocity at time t - INPUT: - t - (optional) time at which to get the velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vx(t) - HISTORY: - 2010-11-30 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - if len(thiso[:,0]) == 2: - return thiso[1,:] - if len(thiso[:,0]) != 4 and len(thiso[:,0]) != 6: - raise AttributeError("orbit must track azimuth to use vx()") - elif len(thiso[:,0]) == 4: - theta= thiso[3,:] - else: - theta= thiso[5,:] - return thiso[1,:]*nu.cos(theta)-thiso[2,:]*nu.sin(theta) - - @physical_conversion('velocity') - def vy(self,*args,**kwargs): - """ - NAME: - vy - PURPOSE: - return y velocity at time t - INPUT: - t - (optional) time at which to get the velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vy(t) - HISTORY: - 2010-11-30 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - if len(thiso[:,0]) != 4 and len(thiso[:,0]) != 6: - raise AttributeError("orbit must track azimuth to use vx()") - elif len(thiso[:,0]) == 4: - theta= thiso[3,:] - else: - theta= thiso[5,:] - return thiso[2,:]*nu.cos(theta)+thiso[1,:]*nu.sin(theta) - - @physical_conversion('velocity') - def vphi(self,*args,**kwargs): - """ - NAME: - vphi - PURPOSE: - return angular velocity - INPUT: - t - (optional) time at which to get the angular velocity - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - OUTPUT: - vphi(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - return thiso[2,:]/thiso[0,:] - - @physical_conversion('angle_deg') - def ra(self,*args,**kwargs): - """ - NAME: - ra - PURPOSE: - return the right ascension - INPUT: - t - (optional) time at which to get ra - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - ra(t) - HISTORY: - 2011-02-23 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'ra') - radec= self._radec(*args,**kwargs) - return radec[:,0] - - @physical_conversion('angle_deg') - def dec(self,*args,**kwargs): - """ - NAME: - dec - PURPOSE: - return the declination - INPUT: - t - (optional) time at which to get dec - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - dec(t) - HISTORY: - 2011-02-23 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'dec') - radec= self._radec(*args,**kwargs) - return radec[:,1] - - @physical_conversion('angle_deg') - def ll(self,*args,**kwargs): - """ - NAME: - ll - PURPOSE: - return Galactic longitude - INPUT: - t - (optional) time at which to get ll - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - l(t) - HISTORY: - 2011-02-23 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'ll') - lbd= self._lbd(*args,**kwargs) - return lbd[:,0] - - @physical_conversion('angle_deg') - def bb(self,*args,**kwargs): - """ - NAME: - bb - PURPOSE: - return Galactic latitude - INPUT: - t - (optional) time at which to get bb - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - b(t) - HISTORY: - 2011-02-23 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'bb') - lbd= self._lbd(*args,**kwargs) - return lbd[:,1] - - @physical_conversion('position_kpc') - def dist(self,*args,**kwargs): - """ - NAME: - dist - PURPOSE: - return distance from the observer in kpc - INPUT: - t - (optional) time at which to get dist - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - dist(t) in kpc - HISTORY: - 2011-02-23 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'dist') - lbd= self._lbd(*args,**kwargs) - return lbd[:,2].astype('float64') - - @physical_conversion('proper-motion_masyr') - def pmra(self,*args,**kwargs): - """ - NAME: - pmra - PURPOSE: - return proper motion in right ascension (in mas/yr) - INPUT: - t - (optional) time at which to get pmra - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - pm_ra(t) in mas / yr - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'pmra') - _check_voSet(self,kwargs,'pmra') - pmrapmdec= self._pmrapmdec(*args,**kwargs) - return pmrapmdec[:,0] - - @physical_conversion('proper-motion_masyr') - def pmdec(self,*args,**kwargs): - """ - NAME: - pmdec - PURPOSE: - return proper motion in declination (in mas/yr) - INPUT: - t - (optional) time at which to get pmdec - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - pm_dec(t) in mas/yr - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'pmdec') - _check_voSet(self,kwargs,'pmdec') - pmrapmdec= self._pmrapmdec(*args,**kwargs) - return pmrapmdec[:,1] - - @physical_conversion('proper-motion_masyr') - def pmll(self,*args,**kwargs): - """ - NAME: - pmll - PURPOSE: - return proper motion in Galactic longitude (in mas/yr) - INPUT: - t - (optional) time at which to get pmll - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - pm_l(t) in mas/yr - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'pmll') - _check_voSet(self,kwargs,'pmll') - lbdvrpmllpmbb= self._lbdvrpmllpmbb(*args,**kwargs) - return lbdvrpmllpmbb[:,4] - - @physical_conversion('proper-motion_masyr') - def pmbb(self,*args,**kwargs): - """ - NAME: - pmbb - PURPOSE: - return proper motion in Galactic latitude (in mas/yr) - INPUT: - t - (optional) time at which to get pmbb - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - pm_b(t) in mas/yr - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'pmbb') - _check_voSet(self,kwargs,'pmbb') - lbdvrpmllpmbb= self._lbdvrpmllpmbb(*args,**kwargs) - return lbdvrpmllpmbb[:,5] - - @physical_conversion('velocity_kms') - def vlos(self,*args,**kwargs): - """ - NAME: - vlos - PURPOSE: - return the line-of-sight velocity (in km/s) - INPUT: - t - (optional) time at which to get vlos - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - vlos(t) in km/s - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'vlos') - _check_voSet(self,kwargs,'vlos') - lbdvrpmllpmbb= self._lbdvrpmllpmbb(*args,**kwargs) - return lbdvrpmllpmbb[:,3] - - @physical_conversion('position_kpc') - def helioX(self,*args,**kwargs): - """ - NAME: - helioX - PURPOSE: - return Heliocentric Galactic rectangular x-coordinate (aka "X") - INPUT: - t - (optional) time at which to get X - obs=[X,Y,Z] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - helioX(t) in kpc - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'helioX') - X, Y, Z= self._helioXYZ(*args,**kwargs) - return X - - @physical_conversion('position_kpc') - def helioY(self,*args,**kwargs): - """ - NAME: - helioY - PURPOSE: - return Heliocentric Galactic rectangular y-coordinate (aka "Y") - INPUT: - t - (optional) time at which to get Y - obs=[X,Y,Z] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - helioY(t) in kpc - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'helioY') - X, Y, Z= self._helioXYZ(*args,**kwargs) - return Y - - @physical_conversion('position_kpc') - def helioZ(self,*args,**kwargs): - """ - NAME: - helioZ - PURPOSE: - return Heliocentric Galactic rectangular z-coordinate (aka "Z") - INPUT: - t - (optional) time at which to get Z - obs=[X,Y,Z] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - OUTPUT: - helioZ(t) in kpc - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'helioZ') - X, Y, Z= self._helioXYZ(*args,**kwargs) - return Z - - @physical_conversion('velocity_kms') - def U(self,*args,**kwargs): - """ - NAME: - U - PURPOSE: - return Heliocentric Galactic rectangular x-velocity (aka "U") - INPUT: - t - (optional) time at which to get U - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - U(t) in km/s - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'U') - _check_voSet(self,kwargs,'U') - X, Y, Z, U, V, W= self._XYZvxvyvz(*args,**kwargs) - return U - - @physical_conversion('velocity_kms') - def V(self,*args,**kwargs): - """ - NAME: - V - PURPOSE: - return Heliocentric Galactic rectangular y-velocity (aka "V") - INPUT: - t - (optional) time at which to get U - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - V(t) in km/s - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'V') - _check_voSet(self,kwargs,'V') - X, Y, Z, U, V, W= self._XYZvxvyvz(*args,**kwargs) - return V - - @physical_conversion('velocity_kms') - def W(self,*args,**kwargs): - """ - NAME: - W - PURPOSE: - return Heliocentric Galactic rectangular z-velocity (aka "W") - INPUT: - t - (optional) time at which to get W - obs=[X,Y,Z,vx,vy,vz] - (optional) position and velocity of observer - (in kpc and km/s) (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - W(t) in km/s - HISTORY: - 2011-02-24 - Written - Bovy (NYU) - """ - _check_roSet(self,kwargs,'W') - _check_voSet(self,kwargs,'W') - X, Y, Z, U, V, W= self._XYZvxvyvz(*args,**kwargs) - return W - - def SkyCoord(self,*args,**kwargs): - """ - NAME: - SkyCoord - PURPOSE: - return the position and velocity as an astropy SkyCoord - INPUT: - t - (optional) time at which to get the position - obs=[X,Y,Z] - (optional) position of observer (in kpc) - (default=Object-wide default) - OR Orbit object that corresponds to the orbit - of the observer - Y is ignored and always assumed to be zero - ro= distance in kpc corresponding to R=1. (default=Object-wide default) - vo= velocity in km/s corresponding to v=1. (default=Object-wide default) - OUTPUT: - SkyCoord(t) - HISTORY: - 2015-06-02 - Written - Bovy (IAS) - """ - kwargs.pop('quantity',None) # rm useless keyword to no conflict later - _check_roSet(self,kwargs,'SkyCoord') - radec= self._radec(*args,**kwargs) - tdist= self.dist(quantity=False,*args,**kwargs) - if not _APY3: # pragma: no cover - return coordinates.SkyCoord(radec[:,0]*units.degree, - radec[:,1]*units.degree, - distance=tdist*units.kpc, - frame='icrs') - _check_voSet(self,kwargs,'SkyCoord') - pmrapmdec= self._pmrapmdec(*args,**kwargs) - vlos= self._lbdvrpmllpmbb(*args,**kwargs)[:,3] - # Also return the Galactocentric frame used - v_sun= coordinates.CartesianDifferential(\ - nu.array([-self._solarmotion[0], - self._solarmotion[1]+self._vo, - self._solarmotion[2]])*units.km/units.s) - return coordinates.SkyCoord(radec[:,0]*units.degree, - radec[:,1]*units.degree, - distance=tdist*units.kpc, - pm_ra_cosdec=pmrapmdec[:,0]\ - *units.mas/units.yr, - pm_dec=pmrapmdec[:,1]*units.mas/units.yr, - radial_velocity=vlos*units.km/units.s, - frame='icrs', - galcen_distance=\ - nu.sqrt(self._ro**2.+self._zo**2.)\ - *units.kpc, - z_sun=self._zo*units.kpc, - galcen_v_sun=v_sun) - - def _radec(self,*args,**kwargs): - """Calculate ra and dec""" - return _radec(self,self(*args,**kwargs),*args,**kwargs) - - def _pmrapmdec(self,*args,**kwargs): - """Calculate pmra and pmdec""" - return _pmrapmdec(self,self(*args,**kwargs),*args,**kwargs) - - def _lbd(self,*args,**kwargs): - """Calculate l,b, and d""" - return _lbd(self,self(*args,**kwargs),*args,**kwargs) - - def _helioXYZ(self,*args,**kwargs): - """Calculate heliocentric rectangular coordinates""" - return _helioXYZ(self,self(*args,**kwargs),*args,**kwargs) - - def _lbdvrpmllpmbb(self,*args,**kwargs): - """Calculate l,b,d,vr,pmll,pmbb""" - return _lbdvrpmllpmbb(self,self(*args,**kwargs),*args,**kwargs) - - def _XYZvxvyvz(self,*args,**kwargs): - """Calculate X,Y,Z,U,V,W""" - return _XYZvxvyvz(self,self(*args,**kwargs),*args,**kwargs) - - def Jacobi(self,Omega,t=0.,pot=None): - """ - NAME: - Jacobi - PURPOSE: - calculate the Jacobi integral E - Omega L - INPUT: - Omega - pattern speed - t= time at which to evaluate the Jacobi integral - Pot= potential instance or list of such instances - OUTPUT: - Jacobi integral - HISTORY: - 2011-04-18 - Written - Bovy (NYU) - """ - raise NotImplementedError("'Jacobi' for this Orbit type is not implemented yet") - - @physical_conversion('action') - def L(self,*args,**kwargs): - """ - NAME: - L - PURPOSE: - calculate the angular momentum - INPUT: - (none) - OUTPUT: - angular momentum - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - #Make sure you are not using physical coordinates - old_physical= kwargs.get('use_physical',None) - kwargs['use_physical']= False - Omega= kwargs.pop('Omega',None) - t= kwargs.pop('t',None) - thiso= self(*args,**kwargs) - if not len(thiso.shape) == 2: thiso= thiso.reshape((thiso.shape[0],1)) - if len(thiso[:,0]) < 3: - raise AttributeError("'linearOrbit has no angular momentum") - elif len(thiso[:,0]) == 3 or len(thiso[:,0]) == 4: - if Omega is None: - out= thiso[0,:]*thiso[2,:] - else: - out= thiso[0,:]*(thiso[2,:]-Omega*t*thiso[0,:]) - elif len(thiso[:,0]) == 5: - raise AttributeError("You must track the azimuth to get the angular momentum of a 3D Orbit") - else: #len(thiso[:,0]) == 6 - vx= self.vx(*args,**kwargs) - vy= self.vy(*args,**kwargs) - vz= self.vz(*args,**kwargs) - x= self.x(*args,**kwargs) - y= self.y(*args,**kwargs) - z= self.z(*args,**kwargs) - out= nu.zeros((len(x),3)) - out[:,0]= y*vz-z*vy - out[:,1]= z*vx-x*vz - out[:,2]= x*vy-y*vx - if not old_physical is None: - kwargs['use_physical']= old_physical - else: - kwargs.pop('use_physical') - return out - - @physical_conversion('action') - def Lz(self,*args,**kwargs): - """ - NAME: - Lz - PURPOSE: - calculate the z component of the angular momentum - INPUT: - (none) - OUTPUT: - z-component of the angular momentum - HISTORY: - 2018-08-29 - Written - Bovy (UofT) - """ - thiso= self(*args,**kwargs) - return thiso[0]*thiso[2] - - def _resetaA(self,pot=None,type=None): - """ - NAME: - _resetaA - PURPOSE: - re-set up an actionAngle module for this Orbit - ONLY TO BE CALLED FROM WITHIN SETUPAA - INPUT: - pot - potential - OUTPUT: - True if reset happened, False otherwise - HISTORY: - 2012-06-01 - Written - Bovy (IAS) - """ - if (not pot is None and pot != self._aAPot) \ - or (not type is None and type != self._aAType): - delattr(self,'_aA') - return True - else: - pass #Already set up - - def _setupaA(self,pot=None,type='staeckel',**kwargs): - """ - NAME: - _setupaA - PURPOSE: - set up an actionAngle module for this Orbit - INPUT: - pot - potential - type= ('staeckel') type of actionAngle module to use - 1) 'adiabatic' - 2) 'staeckel' - 3) 'isochroneApprox' - 4) 'spherical' - OUTPUT: - HISTORY: - 2010-11-30 - Written - Bovy (NYU) - 2013-11-27 - Re-written in terms of new actionAngle modules - Bovy (IAS) - 2017-12-25 - Changed default method to 'staeckel' and automatic delta estimation - Bovy (UofT) - """ - if hasattr(self,'_aA'): - if not self._resetaA(pot=pot,type=type): return None - if (len(self.vxvv) == 3 or len(self.vxvv) == 4) \ - and (type == 'staeckel'): - # No reason to do Staeckel, don't include adiabatic here, bc - # rperi, rap, e analytic calc. requires special adiabatic funcs - type= 'spherical' - if pot is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - self._aAPot= pot - self._aAType= type - #Setup - if self._aAType.lower() == 'adiabatic': - self._aA= actionAngle.actionAngleAdiabatic(pot=self._aAPot, - **kwargs) - elif self._aAType.lower() == 'staeckel': - try: - delta= \ - kwargs.pop('delta', - actionAngle.estimateDeltaStaeckel(\ - self._aAPot,self.R(use_physical=False), - self.z(use_physical=False)+(nu.fabs(self.z(use_physical=False)) < 1e-8) * (2.*(self.z(use_physical=False) >= 0)-1.)*1e-10)) # try to make sure this is not 0 - except PotentialError as e: - if 'deriv' in str(e): - raise PotentialError('Automagic calculation of delta parameter for Staeckel approximation failed because the necessary second derivatives of the given potential are not implemented; set delta= explicitly') - elif 'non-axi' in str(e): - raise PotentialError('Automagic calculation of delta parameter for Staeckel approximation failed because the given potential is not axisymmetric; pass an axisymmetric potential instead') - else: #pragma: no cover - raise - if delta < 1e-6: - self._setupaA(pot=pot,type='spherical') - else: - self._aA= actionAngle.actionAngleStaeckel(pot=self._aAPot, - delta=delta, - **kwargs) - elif self._aAType.lower() == 'isochroneapprox': - from galpy.actionAngle import actionAngleIsochroneApprox - self._aA= actionAngleIsochroneApprox(pot=self._aAPot, - **kwargs) - elif self._aAType.lower() == 'spherical': - self._aA= actionAngle.actionAngleSpherical(pot=self._aAPot, - **kwargs) - return None - - def _xw(self,*args,**kwargs): #pragma: no cover - """ - NAME: - xw - PURPOSE: - return the Fourier transform of xx - INPUT: - t - (optional) time at which to get xw - OUTPUT: - xw(t) - HISTORY: - 2011-01-04 - Written - Bovy (NYU) - """ - #BOVY: REPLACE WITH CALCULATION FUNCTION - x= self.x(self.t) - xw= nu.fft.fft(x)#-nu.mean(x)) - xw= nu.abs(xw[0:len(xw)/2])*(self.t[1]-self.t[0])/(self.t[-1]-self.t[0]) - return xw - - def _plotxw(self,*args,**kwargs): #pragma: no cover - """ - NAME: - plotxw - PURPOSE: - plot the spectrum of x - INPUT: - bovy_plot.bovy_plot args and kwargs - OUTPUT: - x(t) - HISTORY: - 2010-09-21 - Written - Bovy (NYU) - """ - xw= self.xw() - #BOVY: CHECK THAT THIS IS CORRECT - plot.bovy_plot(2.*m.pi*nu.fft.fftfreq(len(self.t), - d=(self.t[1]-self.t[0]))\ - [0:len(xw)], - xw,*args,**kwargs) - - def __call__(self,*args,**kwargs): - """ - NAME: - __call__ - PURPOSE: - return the orbit vector at time t - INPUT: - t - desired time - OUTPUT: - [R,vR,vT,z,vz(,phi)] or [R,vR,vT(,phi)] depending on the orbit - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - if len(args) == 0: - return nu.array(self.vxvv) - else: - t= args[0] - # Parse t - if _APY_LOADED and isinstance(t,units.Quantity): - t= t.to(units.Gyr).value\ - /bovy_conversion.time_in_Gyr(self._vo,self._ro) - elif hasattr(self,'_integrate_t_asQuantity') \ - and self._integrate_t_asQuantity \ - and not nu.all(t == self.t): - warnings.warn("You specified integration times as a Quantity, but are evaluating at times not specified as a Quantity; assuming that time given is in natural (internal) units (multiply time by unit to get output at physical time)",galpyWarning) - if hasattr(self,'t') and nu.all(t == self.t): # Common case where one wants all integrated times - return self.orbit.T - elif isinstance(t,(int,float)) and hasattr(self,'t') \ - and t in list(self.t): - return self.orbit[list(self.t).index(t),:] - else: - if isinstance(t,(int,float)): - nt= 1 - t= [t] - else: - nt= len(t) - dim= len(self.vxvv) - try: - self._setupOrbitInterp() - except: - out= nu.zeros((dim,nt)) - for jj in range(nt): - try: - indx= list(self.t).index(t[jj]) - except ValueError: - raise LookupError("Orbit interpolaton failed; integrate on finer grid") - for ii in range(dim): - out[ii,jj]= self.orbit[indx,ii] - return out #should always have nt > 1, bc otherwise covered by above - out= [] - if _OLD_SCIPY and not isinstance(self._orbInterp[0],_fakeInterp) \ - and nu.any((nu.array(t) < self._orbInterp[0]._data[3])\ - +(nu.array(t) > self._orbInterp[0]._data[4])): #pragma: no cover - raise ValueError("One or more requested time is not within the integrated range") - if dim == 4 or dim == 6: - #Unpack interpolated x and y to R and phi - x= self._orbInterp[0](t) - y= self._orbInterp[-1](t) - R= nu.sqrt(x*x+y*y) - phi= nu.arctan2(y,x) % (2.*nu.pi) - for ii in range(dim): - if ii == 0: - out.append(R) - elif ii == dim-1: - out.append(phi) - else: - out.append(self._orbInterp[ii](t)) - else: - for ii in range(dim): - out.append(self._orbInterp[ii](t)) - if nt == 1: - return nu.array(out).reshape(dim) - else: - return nu.array(out).reshape((dim,nt)) - - def plot(self,*args,**kwargs): - """ - NAME: - plot - PURPOSE: - plot aspects of an Orbit - INPUT: - bovy_plot args and kwargs - ro= (Object-wide default) physical scale for distances to use to convert - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - - +kwargs for ra,dec,ll,bb, etc. functions - OUTPUT: - plot - HISTORY: - 2010-07-26 - Written - Bovy (NYU) - 2010-09-22 - Adapted to more general framework - Bovy (NYU) - 2013-11-29 - added ra,dec kwargs and other derived quantities - Bovy (IAS) - 2014-06-11 - Support for plotting in physical coordinates - Bovy (IAS) - 2017-11-28 - Allow arbitrary functions of time to be plotted - Bovy (UofT) - """ - if (kwargs.get('use_physical',False) \ - and kwargs.get('ro',self._roSet)) or \ - (not 'use_physical' in kwargs \ - and kwargs.get('ro',self._roSet)): - labeldict= {'t':r'$t\ (\mathrm{Gyr})$','R':r'$R\ (\mathrm{kpc})$', - 'vR':r'$v_R\ (\mathrm{km\,s}^{-1})$', - 'vT':r'$v_T\ (\mathrm{km\,s}^{-1})$', - 'z':r'$z\ (\mathrm{kpc})$', - 'vz':r'$v_z\ (\mathrm{km\,s}^{-1})$','phi':r'$\phi$', - 'r':r'$r\ (\mathrm{kpc})$', - 'x':r'$x\ (\mathrm{kpc})$','y':r'$y\ (\mathrm{kpc})$', - 'vx':r'$v_x\ (\mathrm{km\,s}^{-1})$', - 'vy':r'$v_y\ (\mathrm{km\,s}^{-1})$', - 'E':r'$E\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Ez':r'$E_z\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'ER':r'$E_R\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Enorm':r'$E(t)/E(0.)$', - 'Eznorm':r'$E_z(t)/E_z(0.)$', - 'ERnorm':r'$E_R(t)/E_R(0.)$', - 'Jacobi':r'$E-\Omega_p\,L\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$'} - else: - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'r':r'$r$', - 'x':r'$x$','y':r'$y$','vx':r'$v_x$','vy':r'$v_y$', - 'E':r'$E$','Enorm':r'$E(t)/E(0.)$', - 'Ez':r'$E_z$','Eznorm':r'$E_z(t)/E_z(0.)$', - 'ER':r'$E_R$','ERnorm':r'$E_R(t)/E_R(0.)$', - 'Jacobi':r'$E-\Omega_p\,L$', - 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$'} - labeldict.update({'ra':r'$\alpha\ (\mathrm{deg})$', - 'dec':r'$\delta\ (\mathrm{deg})$', - 'll':r'$l\ (\mathrm{deg})$', - 'bb':r'$b\ (\mathrm{deg})$', - 'dist':r'$d\ (\mathrm{kpc})$', - 'pmra':r'$\mu_\alpha\ (\mathrm{mas\,yr}^{-1})$', - 'pmdec':r'$\mu_\delta\ (\mathrm{mas\,yr}^{-1})$', - 'pmll':r'$\mu_l\ (\mathrm{mas\,yr}^{-1})$', - 'pmbb':r'$\mu_b\ (\mathrm{mas\,yr}^{-1})$', - 'vlos':r'$v_\mathrm{los}\ (\mathrm{km\,s}^{-1})$', - 'helioX':r'$X\ (\mathrm{kpc})$', - 'helioY':r'$Y\ (\mathrm{kpc})$', - 'helioZ':r'$Z\ (\mathrm{kpc})$', - 'U':r'$U\ (\mathrm{km\,s}^{-1})$', - 'V':r'$V\ (\mathrm{km\,s}^{-1})$', - 'W':r'$W\ (\mathrm{km\,s}^{-1})$'}) - # Cannot be using Quantity output - kwargs['quantity']= False - #Defaults - if not 'd1' in kwargs and not 'd2' in kwargs: - if len(self.vxvv) == 3: - d1= 'R' - d2= 'vR' - elif len(self.vxvv) == 4: - d1= 'x' - d2= 'y' - elif len(self.vxvv) == 2: - d1= 'x' - d2= 'vx' - elif len(self.vxvv) == 5 or len(self.vxvv) == 6: - d1= 'R' - d2= 'z' - elif not 'd1' in kwargs: - d2= kwargs.pop('d2') - d1= 't' - elif not 'd2' in kwargs: - d1= kwargs.pop('d1') - d2= 't' - else: - d1= kwargs.pop('d1') - d2= kwargs.pop('d2') - x= self._parse_plot_quantity(d1,**kwargs) - y= self._parse_plot_quantity(d2,**kwargs) - kwargs.pop('ro',None) - kwargs.pop('vo',None) - kwargs.pop('obs',None) - kwargs.pop('use_physical',None) - kwargs.pop('pot',None) - kwargs.pop('OmegaP',None) - kwargs.pop('quantity',None) - #Plot - if not 'xlabel' in kwargs: - kwargs['xlabel']= labeldict.get(d1,'\mathrm{No\ xlabel\ specified}') - if not 'ylabel' in kwargs: - kwargs['ylabel']= labeldict.get(d2,'\mathrm{No\ ylabel\ specified}') - return plot.bovy_plot(x,y,*args,**kwargs) - - def plot3d(self,*args,**kwargs): - """ - NAME: - plot3d - PURPOSE: - plot 3D aspects of an Orbit - INPUT: - ro= (Object-wide default) physical scale for distances to use to convert - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - - bovy_plot args and kwargs - OUTPUT: - plot - HISTORY: - 2010-07-26 - Written - Bovy (NYU) - 2010-09-22 - Adapted to more general framework - Bovy (NYU) - 2010-01-08 - Adapted to 3D - Bovy (NYU) - 2013-11-29 - added ra,dec kwargs and other derived quantities - Bovy (IAS) - 2014-06-11 - Support for plotting in physical coordinates - Bovy (IAS) - 2017-11-28 - Allow arbitrary functions of time to be plotted - Bovy (UofT) - """ - if (kwargs.get('use_physical',False) \ - and kwargs.get('ro',self._roSet)) or \ - (not 'use_physical' in kwargs \ - and kwargs.get('ro',self._roSet)): - labeldict= {'t':r'$t\ (\mathrm{Gyr})$','R':r'$R\ (\mathrm{kpc})$', - 'vR':r'$v_R\ (\mathrm{km\,s}^{-1})$', - 'vT':r'$v_T\ (\mathrm{km\,s}^{-1})$', - 'z':r'$z\ (\mathrm{kpc})$', - 'vz':r'$v_z\ (\mathrm{km\,s}^{-1})$','phi':r'$\phi$', - 'r':r'$r\ (\mathrm{kpc})$', - 'x':r'$x\ (\mathrm{kpc})$','y':r'$y\ (\mathrm{kpc})$', - 'vx':r'$v_x\ (\mathrm{km\,s}^{-1})$', - 'vy':r'$v_y\ (\mathrm{km\,s}^{-1})$'} - else: - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'r':r'$r$','x':r'$x$','y':r'$y$', - 'vx':r'$v_x$','vy':r'$v_y$'} - labeldict.update({'ra':r'$\alpha\ (\mathrm{deg})$', - 'dec':r'$\delta\ (\mathrm{deg})$', - 'll':r'$l\ (\mathrm{deg})$', - 'bb':r'$b\ (\mathrm{deg})$', - 'dist':r'$d\ (\mathrm{kpc})$', - 'pmra':r'$\mu_\alpha\ (\mathrm{mas\,yr}^{-1})$', - 'pmdec':r'$\mu_\delta\ (\mathrm{mas\,yr}^{-1})$', - 'pmll':r'$\mu_l\ (\mathrm{mas\,yr}^{-1})$', - 'pmbb':r'$\mu_b\ (\mathrm{mas\,yr}^{-1})$', - 'vlos':r'$v_\mathrm{los}\ (\mathrm{km\,s}^{-1})$', - 'helioX':r'$X\ (\mathrm{kpc})$', - 'helioY':r'$Y\ (\mathrm{kpc})$', - 'helioZ':r'$Z\ (\mathrm{kpc})$', - 'U':r'$U\ (\mathrm{km\,s}^{-1})$', - 'V':r'$V\ (\mathrm{km\,s}^{-1})$', - 'W':r'$W\ (\mathrm{km\,s}^{-1})$'}) - # Cannot be using Quantity output - kwargs['quantity']= False - #Defaults - if not 'd1' in kwargs and not 'd2' in kwargs and not 'd3' in kwargs: - if len(self.vxvv) == 3: - d1= 'R' - d2= 'vR' - d3= 'vT' - elif len(self.vxvv) == 4: - d1= 'x' - d2= 'y' - d3= 'vR' - elif len(self.vxvv) == 2: - raise AttributeError("Cannot plot 3D aspects of 1D orbits") - elif len(self.vxvv) == 5: - d1= 'R' - d2= 'vR' - d3= 'z' - elif len(self.vxvv) == 6: - d1= 'x' - d2= 'y' - d3= 'z' - elif not ('d1' in kwargs and 'd2' in kwargs and 'd3' in kwargs): - raise AttributeError("Please provide 'd1', 'd2', and 'd3'") - else: - d1= kwargs.pop('d1') - d2= kwargs.pop('d2') - d3= kwargs.pop('d3') - x= self._parse_plot_quantity(d1,**kwargs) - y= self._parse_plot_quantity(d2,**kwargs) - z= self._parse_plot_quantity(d3,**kwargs) - kwargs.pop('ro',None) - kwargs.pop('vo',None) - kwargs.pop('obs',None) - kwargs.pop('use_physical',None) - kwargs.pop('quantity',None) - #Plot - if not 'xlabel' in kwargs: - kwargs['xlabel']= labeldict.get(d1,'\mathrm{No\ xlabel\ specified}') - if not 'ylabel' in kwargs: - kwargs['ylabel']= labeldict.get(d2,'\mathrm{No\ ylabel\ specified}') - if not 'zlabel' in kwargs: - kwargs['zlabel']= labeldict.get(d3,'\mathrm{No\ zlabel\ specified}') - return plot.bovy_plot3d(x,y,z,*args,**kwargs) - - def _parse_plot_quantity(self,quant,**kwargs): - """Internal function to parse a quantity to be plotted based on input data""" - # Cannot be using Quantity output - kwargs['quantity']= False - if callable(quant): - return quant(self.t) - def _eval(q): - # Check those that don't have the exact name of the function - if q == 't': - return self.time(self.t,**kwargs) - elif q == 'Enorm': - return self.E(self.t,**kwargs)/self.E(0.,**kwargs) - elif q == 'Eznorm': - return self.Ez(self.t,**kwargs)/self.Ez(0.,**kwargs) - elif q == 'ERnorm': - return self.ER(self.t,**kwargs)/self.ER(0.,**kwargs) - elif q == 'Jacobinorm': - return self.Jacobi(self.t,**kwargs)/self.Jacobi(0.,**kwargs) - else: # these are exact, e.g., 'x' for self.x - return self.__getattribute__(q)(self.t,**kwargs) - try: - return _eval(quant) - except AttributeError: pass - try: - import numexpr - except ImportError: #pragma: no cover - raise ImportError('Parsing the quantity to be plotted failed; if you are trying to plot an expression, please make sure to install numexpr first') - # Figure out the variables in the expression to be computed to plot - try: - vars= numexpr.NumExpr(quant).input_names - except TypeError as err: - raise TypeError('Parsing the expression {} failed, with error message:\n"{}"'.format(quant,err)) - # Construct dictionary of necessary parameters - vars_dict= {} - for var in vars: - vars_dict[var]= _eval(var) - return numexpr.evaluate(quant,local_dict=vars_dict) - - def plotR(self,*args,**kwargs): - """ - NAME: - plotR - PURPOSE: - plot R(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'R' - return self.plot(*args,**kwargs) - - def plotz(self,*args,**kwargs): - """ - NAME: - plotz - PURPOSE: - plot z(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'z' - return self.plot(*args,**kwargs) - - def plotx(self,*args,**kwargs): - """ - NAME: - plotx - PURPOSE: - plot x(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'x' - return self.plot(*args,**kwargs) - - def plotvx(self,*args,**kwargs): - """ - NAME: - plotvx - PURPOSE: - plot vx(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'vx' - return self.plot(*args,**kwargs) - - def ploty(self,*args,**kwargs): - """ - NAME: - ploty - PURPOSE: - plot y(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'y' - return self.plot(*args,**kwargs) - - def plotvy(self,*args,**kwargs): - """ - NAME: - plotvy - PURPOSE: - plot vy(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'vy' - return self.plot(*args,**kwargs) - - def plotvR(self,*args,**kwargs): - """ - NAME: - plotvR - PURPOSE: - plot vR(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'vR' - return self.plot(*args,**kwargs) - - def plotvT(self,*args,**kwargs): - """ - NAME: - plotvT - PURPOSE: - plot vT(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'vT' - return self.plot(*args,**kwargs) - - def plotphi(self,*args,**kwargs): - """ - NAME: - plotphi - PURPOSE: - plot \phi(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'phi' - return self.plot(*args,**kwargs) - - def plotvz(self,*args,**kwargs): - """ - NAME: - plotvz - PURPOSE: - plot vz(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-07-10 - Written - Bovy (NYU) - """ - kwargs['d2']= 'vz' - return self.plot(*args,**kwargs) - - def plotE(self,*args,**kwargs): - """ - NAME: - plotE - PURPOSE: - plot E(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'Enorm' - else: - kwargs['d2']= 'E' - return self.plot(*args,**kwargs) - - def plotJacobi(self,*args,**kwargs): - """ - NAME: - plotE - PURPOSE: - plot Jacobi(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'Jacobinorm' - else: - kwargs['d2']= 'Jacobi' - return self.plot(*args,**kwargs) - - def _setupOrbitInterp(self): - if not hasattr(self,"_orbInterp"): - # First check that times increase - if hasattr(self,"t"): #Orbit has been integrated - if self.t[1] < self.t[0]: #must be backward - sindx= nu.argsort(self.t) - # sort - self.t= self.t[sindx] - self.orbit= self.orbit[sindx] - usindx= nu.argsort(sindx) # to later unsort - orbInterp= [] - for ii in range(len(self.vxvv)): - if (len(self.vxvv) == 4 or len(self.vxvv) == 6) and ii == 0: - #Interpolate x and y rather than R and phi to avoid issues w/ phase wrapping - if not hasattr(self,"t"): #Orbit has not been integrated - orbInterp.append(_fakeInterp(self.vxvv[0]*nu.cos(self.vxvv[-1]))) - else: - orbInterp.append(interpolate.InterpolatedUnivariateSpline(\ - self.t,self.orbit[:,0]*nu.cos(self.orbit[:,-1]), - **_KWINTERP)) - elif (len(self.vxvv) == 4 or len(self.vxvv) == 6) and \ - ii == len(self.vxvv)-1: - if not hasattr(self,"t"): #Orbit has not been integrated - orbInterp.append(_fakeInterp(self.vxvv[0]*nu.sin(self.vxvv[-1]))) - else: - orbInterp.append(interpolate.InterpolatedUnivariateSpline(\ - self.t,self.orbit[:,0]*nu.sin(self.orbit[:,-1]),**_KWINTERP)) - else: - if not hasattr(self,"t"): #Orbit has not been integrated - orbInterp.append(_fakeInterp(self.vxvv[ii])) - else: - orbInterp.append(interpolate.InterpolatedUnivariateSpline(\ - self.t,self.orbit[:,ii],**_KWINTERP)) - self._orbInterp= orbInterp - try: #unsort - self.t= self.t[usindx] - self.orbit= self.orbit[usindx] - except: pass - return None - - def animate(self,*args,**kwargs): #pragma: no cover - """ - NAME: - animate - PURPOSE: - animate an Orbit - INPUT: - d1= first dimension to plot ('x', 'y', 'R', 'vR', 'vT', 'z', 'vz', ...); can be list with up to three entries for three subplots - d2= second dimension to plot; can be list with up to three entries for three subplots - width= (600) width of output div in px - height= (400) height of output div in px - json_filename= (None) if set, save the data necessary for the figure in this filename (e.g., json_filename= 'orbit_data/orbit.json'); this path is also used in the output HTML, so needs to be accessible - load_jslibs= (True) if True, load the require and jQuery Javascript libraries (necessary in Jupyterlab, not necessary but harmless in notebooks; if embedding on a webpage one typically wants to load these libraries in the header) - ro= (Object-wide default) physical scale for distances to use to convert - vo= (Object-wide default) physical scale for velocities to use to convert - use_physical= use to override Object-wide default for using a physical scale for output - - +kwargs for ra,dec,ll,bb, etc. functions - OUTPUT: - IPython.display.HTML object with code to animate the orbit; can be directly shown in jupyter notebook or embedded in HTML pages; get a text version of the HTML using the _repr_html_() function - HISTORY: - 2017-09-17-24 - Written - Bovy (UofT) - 2017-11-28 - Allow arbitrary functions of time to be plotted - Bovy (UofT) - """ - try: - from IPython.display import HTML - except ImportError: - raise ImportError("Orbit.animate requires ipython/jupyter to be installed") - if (kwargs.get('use_physical',False) \ - and kwargs.get('ro',self._roSet)) or \ - (not 'use_physical' in kwargs \ - and kwargs.get('ro',self._roSet)): - labeldict= {'t':'t (Gyr)', - 'R':'R (kpc)', - 'vR':'v_R (km/s)', - 'vT':'v_T (km/s)', - 'z':'z (kpc)', - 'vz':'v_z (km/s)', - 'phi':'azimuthal angle', - 'r':'r (kpc)', - 'x':'x (kpc)', - 'y':'y (kpc)', - 'vx':'v_x (km/s)', - 'vy':'v_y (km/s)', - 'E':'E (km^2/s^2)', - 'Ez':'E_z (km^2/s^2)', - 'ER':'E_R (km^2/s^2)', - 'Enorm':'E(t)/E(0.)', - 'Eznorm':'E_z(t)/E_z(0.)', - 'ERnorm':'E_R(t)/E_R(0.)', - 'Jacobi':'E-Omega_p L (km^2/s^2)', - 'Jacobinorm':'(E-Omega_p L)(t)/(E-Omega_p L)(0)'} - else: - labeldict= {'t':'t','R':'R','vR':'v_R','vT':'v_T', - 'z':'z','vz':'v_z','phi':r'azimuthal angle', - 'r':'r', - 'x':'x','y':'y','vx':'v_x','vy':'v_y', - 'E':'E','Enorm':'E(t)/E(0.)', - 'Ez':'E_z','Eznorm':'E_z(t)/E_z(0.)', - 'ER':r'E_R','ERnorm':r'E_R(t)/E_R(0.)', - 'Jacobi':r'E-Omega_p L', - 'Jacobinorm':r'(E-Omega_p L)(t)/(E-Omega_p L)(0)'} - labeldict.update({'ra':'RA (deg)', - 'dec':'Dec (deg)', - 'll':'Galactic lon (deg)', - 'bb':'Galactic lat (deg)', - 'dist':'distance (kpc)', - 'pmra':'pmRA (mas/yr)', - 'pmdec':'pmDec (mas/yr)', - 'pmll':'pmGlon (mas/yr)', - 'pmbb':'pmGlat (mas/yr)', - 'vlos':'line-of-sight vel (km/s)', - 'helioX':'X (kpc)', - 'helioY':'Y (kpc)', - 'helioZ':'Z (kpc)', - 'U':'U (km/s)', - 'V':'V (km/s)', - 'W':'W (km/s)'}) - # Cannot be using Quantity output - kwargs['quantity']= False - #Defaults - if not 'd1' in kwargs and not 'd2' in kwargs: - if len(self.vxvv) == 3: - d1= 'R' - d2= 'vR' - elif len(self.vxvv) == 4: - d1= 'x' - d2= 'y' - elif len(self.vxvv) == 2: - d1= 'x' - d2= 'vx' - elif len(self.vxvv) == 5 or len(self.vxvv) == 6: - d1= 'R' - d2= 'z' - elif not 'd1' in kwargs: - d2= kwargs.pop('d2') - d1= 't' - elif not 'd2' in kwargs: - d1= kwargs.pop('d1') - d2= 't' - else: - d1= kwargs.pop('d1') - d2= kwargs.pop('d2') - xs= [] - ys= [] - xlabels= [] - ylabels= [] - if isinstance(d1,str) or callable(d1): - d1s= [d1] - d2s= [d2] - else: - d1s= d1 - d2s= d2 - if len(d1s) > 3: - raise ValueError('Orbit.animate only works for up to three subplots') - all_xlabel= kwargs.get('xlabel',[None for d in d1]) - all_ylabel= kwargs.get('ylabel',[None for d in d2]) - for d1,d2, xlabel, ylabel in zip(d1s,d2s,all_xlabel,all_ylabel): - #Get x and y for each subplot - x= self._parse_plot_quantity(d1,**kwargs) - y= self._parse_plot_quantity(d2,**kwargs) - xs.append(x) - ys.append(y) - if xlabel is None: - xlabels.append(labeldict.get(d1,'\mathrm{No\ xlabel\ specified}')) - else: - xlabels.append(xlabel) - if ylabel is None: - ylabels.append(labeldict.get(d2,'\mathrm{No\ ylabel\ specified}')) - else: - ylabels.append(ylabel) - kwargs.pop('ro',None) - kwargs.pop('vo',None) - kwargs.pop('obs',None) - kwargs.pop('use_physical',None) - kwargs.pop('pot',None) - kwargs.pop('OmegaP',None) - kwargs.pop('quantity',None) - width= kwargs.pop('width',600) - height= kwargs.pop('height',400) - load_jslibs= kwargs.pop('load_jslibs',True) - if load_jslibs: - load_jslibs_code= """ - - -""".format(json_code=json_code,close_json_code=close_json_code, - divid=self.divid,width=width,height=height, - button_margin_left=button_margin_left, - layout=layout,load_jslibs_code=load_jslibs_code, - setup_trace2=setup_trace2,setup_trace3=setup_trace3, - delete_trace4=delete_trace4,delete_trace6=delete_trace6, - delete_trace3=delete_trace3,delete_trace5=delete_trace5, - update_trace34=update_trace34, - update_trace56=update_trace56)) - -class _fakeInterp(object): - """Fake class to simulate interpolation when orbit was not integrated""" - def __init__(self,x): - self.x= x - def __call__(self,t): - if nu.any(nu.array(t) != 0.): - raise ValueError("Integrate instance before evaluating it at non-zero time") - else: - return nu.array([self.x for i in t]) - diff --git a/galpy/orbit/RZOrbit.py b/galpy/orbit/RZOrbit.py deleted file mode 100644 index 43dcafb34..000000000 --- a/galpy/orbit/RZOrbit.py +++ /dev/null @@ -1,535 +0,0 @@ -import warnings -import math as m -import numpy as nu -from scipy import integrate -from galpy.potential.Potential import _evaluateRforces, _evaluatezforces,\ - evaluatePotentials, evaluateDensities, _check_c -from galpy.util import galpyWarning -import galpy.util.bovy_plot as plot -import galpy.util.bovy_symplecticode as symplecticode -from .FullOrbit import _integrateFullOrbit -from .integrateFullOrbit import _ext_loaded as ext_loaded -from galpy.util.bovy_conversion import physical_conversion -from galpy.util.leung_dop853 import dop853 -from .OrbitTop import OrbitTop -class RZOrbit(OrbitTop): - """Class that holds and integrates orbits in axisymetric potentials - in the (R,z) plane""" - def __init__(self,vxvv=[1.,0.,0.9,0.,0.1],vo=220.,ro=8.0,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): - """ - NAME: - - __init__ - - PURPOSE: - - intialize an RZ-orbit - - INPUT: - - vxvv - initial condition [R,vR,vT,z,vz] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - (none) - - HISTORY: - - 2010-07-10 - Written - Bovy (NYU) - - 2014-06-11 - Added conversion kwargs to physical coordinates - Bovy (IAS) - - """ - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=zo,vo=vo,solarmotion=solarmotion) - return None - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - (none) (get the actual orbit using getOrbit() - HISTORY: - 2010-07-10 - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - self.t= nu.array(t) - self._pot= pot - self.orbit= _integrateRZOrbit(self.vxvv,pot,t,method,dt) - - @physical_conversion('energy') - def E(self,*args,**kwargs): - """ - NAME: - E - PURPOSE: - calculate the energy - INPUT: - t - (optional) time at which to get the radius - pot= RZPotential instance or list thereof - OUTPUT: - energy - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],thiso[3], - t=t,use_physical=False)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2.\ - +thiso[4]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],thiso[3,ii], - t=t[ii],use_physical=False)\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2.\ - +thiso[4,ii]**2./2. for ii in range(len(t))]) - - @physical_conversion('energy') - def ER(self,*args,**kwargs): - """ - NAME: - ER - PURPOSE: - calculate the radial energy - INPUT: - t - (optional) time at which to get the energy - pot= potential instance or list of such instances - OUTPUT: - radial energy - HISTORY: - 2013-11-30 - Written - Bovy (IAS) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],0., - t=t,use_physical=False)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],0., - t=t[ii],use_physical=False)\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2. for ii in range(len(t))]) - - @physical_conversion('energy') - def Ez(self,*args,**kwargs): - """ - NAME: - Ez - PURPOSE: - calculate the vertical energy - INPUT: - t - (optional) time at which to get the energy - pot= potential instance or list of such instances - OUTPUT: - vertical energy - HISTORY: - 2013-11-30 - Written - Bovy (IAS) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatePotentials(pot,thiso[0],thiso[3], - t=t,use_physical=False)\ - -evaluatePotentials(pot,thiso[0],0., - t=t, - use_physical=False)\ - +thiso[4]**2./2. - else: - return nu.array([evaluatePotentials(pot,thiso[0,ii],thiso[3,ii], - t=t[ii],use_physical=False)\ - -evaluatePotentials(pot,thiso[0,ii],0., - t=t[ii],use_physical=False)\ - +thiso[4,ii]**2./2. for ii in range(len(t))]) - - @physical_conversion('energy') - def Jacobi(self,*args,**kwargs): - """ - NAME: - Jacobi - PURPOSE: - calculate the Jacobi integral of the motion - INPUT: - t - (optional) time at which to get the radius - OmegaP= pattern speed of rotating frame (scalar) - pot= potential instance or list of such instances - OUTPUT: - Jacobi integral - HISTORY: - 2011-04-18 - Written - Bovy (NYU) - """ - if not 'OmegaP' in kwargs or kwargs['OmegaP'] is None: - OmegaP= 1. - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - else: - pot= kwargs['pot'] - if isinstance(pot,list): - for p in pot: - if hasattr(p,'OmegaP'): - OmegaP= p.OmegaP() - break - else: - if hasattr(pot,'OmegaP'): - OmegaP= pot.OmegaP() - kwargs.pop('OmegaP',None) - else: - OmegaP= kwargs.pop('OmegaP') - #Make sure you are not using physical coordinates - old_physical= kwargs.get('use_physical',None) - kwargs['use_physical']= False - thiso= self(*args,**kwargs) - out= self.E(*args,**kwargs)-OmegaP*thiso[0]*thiso[2] - if not old_physical is None: - kwargs['use_physical']= old_physical - else: - kwargs.pop('use_physical') - return out - - def e(self,analytic=False,pot=None,**kwargs): - """ - NAME: - e - PURPOSE: - calculate the eccentricity - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - eccentricity - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[0]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate eccentricity") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs)) - - @physical_conversion('position') - def rap(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rap - PURPOSE: - return the apocenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_ap - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[3]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate rap") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return nu.amax(self.rs) - - @physical_conversion('position') - def rperi(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rperi - PURPOSE: - return the pericenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_peri - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[2]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate rperi") - if not hasattr(self,'rs'): - self.rs= nu.sqrt(self.orbit[:,0]**2.+self.orbit[:,3]**2.) - return nu.amin(self.rs) - - @physical_conversion('position') - def zmax(self,analytic=False,pot=None,**kwargs): - """ - NAME: - zmax - PURPOSE: - return the maximum vertical height - INPUT: - OUTPUT: - Z_max - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,**kwargs) - return float(self._aA.EccZmaxRperiRap(self)[1]) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate zmax") - return nu.amax(nu.fabs(self.orbit[:,3])) - - def plotEz(self,*args,**kwargs): - """ - NAME: - plotEz - PURPOSE: - plot Ez(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'Eznorm' - else: - kwargs['d2']= 'Ez' - return self.plot(*args,**kwargs) - - def plotER(self,*args,**kwargs): - """ - NAME: - plotER - PURPOSE: - plot ER(.) along the orbit - INPUT: - bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2014-06-16 - Written - Bovy (IAS) - """ - if kwargs.pop('normed',False): - kwargs['d2']= 'ERnorm' - else: - kwargs['d2']= 'ER' - return self.plot(*args,**kwargs) - - def plotEzJz(self,*args,**kwargs): - """ - NAME: - plotEzJz - PURPOSE: - plot E_z(.)/sqrt(dens(R)) along the orbit - INPUT: - pot= Potential instance or list of instances in which the orbit was - integrated - d1= - plot Ez vs d1: e.g., 't', 'z', 'R', 'vR', 'vT', 'vz' - +bovy_plot.bovy_plot inputs - OUTPUT: - figure to output device - HISTORY: - 2010-08-08 - Written - Bovy (NYU) - """ - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'x':r'$x$','y':r'$y$','vx':r'$v_x$','vy':r'$v_y$'} - if not 'pot' in kwargs: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit first or specify pot=") - else: - pot= kwargs.pop('pot') - d1= kwargs.pop('d1','t') - self.EzJz= [(evaluatePotentials(pot,self.orbit[ii,0],self.orbit[ii,3], - t=self.t[ii],use_physical=False)- - evaluatePotentials(pot,self.orbit[ii,0],0.,t=self.t[ii], - use_physical=False)+ - self.orbit[ii,4]**2./2.)/\ - nu.sqrt(evaluateDensities(pot,self.orbit[ii,0],0., - t=self.t[ii], - use_physical=False))\ - for ii in range(len(self.t))] - if not 'xlabel' in kwargs: - kwargs['xlabel']= labeldict[d1] - if not 'ylabel' in kwargs: - kwargs['ylabel']= r'$E_z/\sqrt{\rho}$' - if d1 == 't': - return plot.bovy_plot(nu.array(self.t), - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'z': - return plot.bovy_plot(self.orbit[:,3], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'R': - return plot.bovy_plot(self.orbit[:,0], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vR': - return plot.bovy_plot(self.orbit[:,1], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vT': - return plot.bovy_plot(self.orbit[:,2], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - elif d1 == 'vz': - return plot.bovy_plot(self.orbit[:,4], - nu.array(self.EzJz)/self.EzJz[0], - *args,**kwargs) - -def _integrateRZOrbit(vxvv,pot,t,method,dt): - """ - NAME: - _integrateRZOrbit - PURPOSE: - integrate an orbit in a Phi(R,z) potential in the (R,z) plane - INPUT: - vxvv - array with the initial conditions stacked like - [R,vR,vT,z,vz]; vR outward! - pot - Potential instance - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - dt - if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - [:,5] array of [R,vR,vT,z,vz] at each t - HISTORY: - 2010-04-16 - Written - Bovy (NYU) - """ - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - if method.lower() == 'leapfrog' \ - or method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \ - or method.lower() == 'rk6_c' or method.lower() == 'symplec4_c' \ - or method.lower() == 'symplec6_c' or method.lower() == 'dopr54_c' or method.lower() == 'dop853_c': - #We hack this by upgrading to a FullOrbit - this_vxvv= nu.zeros(len(vxvv)+1) - this_vxvv[0:len(vxvv)]= vxvv - tmp_out= _integrateFullOrbit(this_vxvv,pot,t,method,dt) - #tmp_out is (nt,6) - out= tmp_out[:,0:5] - elif method.lower() == 'odeint' or method.lower() == 'dop853': - l= vxvv[0]*vxvv[2] - l2= l**2. - init= [vxvv[0],vxvv[1],vxvv[3],vxvv[4]] - if method.lower() == "dop853": - intOut = dop853(_RZEOM, init, t, args=(pot, l2)) - else: - intOut = integrate.odeint(_RZEOM, init, t, args=(pot, l2), - rtol=10. ** -8.) # ,mxstep=100000000) - out= nu.zeros((len(t),5)) - out[:,0]= intOut[:,0] - out[:,1]= intOut[:,1] - out[:,3]= intOut[:,2] - out[:,4]= intOut[:,3] - out[:,2]= l/out[:,0] - #post-process to remove negative radii - neg_radii= (out[:,0] < 0.) - out[neg_radii,0]= -out[neg_radii,0] - return out - -def _RZEOM(y,t,pot,l2): - """ - NAME: - _RZEOM - PURPOSE: - implements the EOM, i.e., the right-hand side of the differential - equation - INPUT: - y - current phase-space position - t - current time - pot - (list of) Potential instance(s) - l2 - angular momentum squared - OUTPUT: - dy/dt - HISTORY: - 2010-04-16 - Written - Bovy (NYU) - """ - return [y[1], - l2/y[0]**3.+_evaluateRforces(pot,y[0],y[2],t=t), - y[3], - _evaluatezforces(pot,y[0],y[2],t=t)] - diff --git a/galpy/orbit/linearOrbit.py b/galpy/orbit/linearOrbit.py deleted file mode 100644 index 382131d2c..000000000 --- a/galpy/orbit/linearOrbit.py +++ /dev/null @@ -1,195 +0,0 @@ -import warnings -import numpy as nu -from scipy import integrate -from .OrbitTop import OrbitTop -from ..potential.linearPotential import \ - _evaluatelinearForces, evaluatelinearPotentials -from ..potential.Potential import _check_c -from galpy.util.leung_dop853 import dop853 -import galpy.util.bovy_symplecticode as symplecticode -from galpy.util.bovy_conversion import physical_conversion -from galpy.util import galpyWarning, galpyWarningVerbose -from .integrateLinearOrbit import integrateLinearOrbit_c, _ext_loaded -ext_loaded= _ext_loaded -class linearOrbit(OrbitTop): - """Class that represents an orbit in a (effectively) one-dimensional potential""" - def __init__(self,vxvv=[1.,0.],vo=220.,ro=8.0): - """ - NAME: - - __init__ - - PURPOSE: - - Initialize a linear orbit - - INPUT: - - vxvv - [x,vx] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - OUTPUT: - - (none) - - HISTORY: - - 2010-07-13 - Written - Bovy (NYU) - - """ - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=None,vo=vo,solarmotion=None) - return None - - def integrate(self,t,pot,method='odeint',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint'= scipy's odeint, or 'leapfrog' - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize (NOT USED FOR LINEAR ORBIT SO FAR) - OUTPUT: - (none) (get the actual orbit using getOrbit() - HISTORY: - 2010-07-13 - Written - Bovy (NYU) - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - self.t= nu.array(t) - self._pot= pot - self.orbit= _integrateLinearOrbit(self.vxvv,pot,t,method,dt) - - @physical_conversion('energy') - def E(self,*args,**kwargs): - """ - NAME: - E - PURPOSE: - calculate the energy - INPUT: - t - (optional) time at which to get the radius - pot= linearPotential instance or list thereof - OUTPUT: - energy - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return evaluatelinearPotentials(pot,thiso[0], - t=t,use_physical=False)\ - +thiso[1]**2./2. - else: - return nu.array([evaluatelinearPotentials(pot,thiso[0,ii], - t=t[ii], - use_physical=False)\ - +thiso[1,ii]**2./2.\ - for ii in range(len(t))]) - - def e(self,analytic=False,pot=None): #pragma: no cover - """ - NAME: - e - PURPOSE: - calculate the eccentricity - INPUT: - OUTPUT: - eccentricity - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - raise AttributeError("linearOrbit does not have an eccentricity") - - def rap(self,analytic=False,pot=None): #pragma: no cover - raise AttributeError("linearOrbit does not have an apocenter") - - def rperi(self,analytic=False,pot=None): #pragma: no cover - raise AttributeError("linearOrbit does not have a pericenter") - - def zmax(self): #pragma: no cover - raise AttributeError("linearOrbit does not have a zmax") - -def _integrateLinearOrbit(vxvv,pot,t,method,dt): - """ - NAME: - integrateLinearOrbit - PURPOSE: - integrate a one-dimensional orbit - INPUT: - vxvv - initial condition [x,vx] - pot - linearPotential or list of linearPotentials - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - OUTPUT: - [:,2] array of [x,vx] at each t - HISTORY: - 2010-07-13- Written - Bovy (NYU) - 2018-10-05- Added support for C integration - Bovy (UofT) - """ - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - if method.lower() == 'leapfrog': - return symplecticode.leapfrog(lambda x,t=t: _evaluatelinearForces(pot,x, - t=t), - nu.array(vxvv), - t,rtol=10.**-8) - elif method.lower() == 'dop853': - return dop853(func=_linearEOM, x=vxvv, t=t, args=(pot,)) - elif ext_loaded and \ - (method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \ - or method.lower() == 'rk6_c' or method.lower() == 'symplec4_c' \ - or method.lower() == 'symplec6_c' or method.lower() == 'dopr54_c' \ - or method.lower() == 'dop853_c'): - warnings.warn("Using C implementation to integrate orbits", - galpyWarningVerbose) - out, msg= integrateLinearOrbit_c(pot,nu.array(vxvv),t,method,dt=dt) - return out - elif method.lower() == 'odeint' or not ext_loaded: - return integrate.odeint(_linearEOM,vxvv,t,args=(pot,),rtol=10.**-8.) - -def _linearEOM(y,t,pot): - """ - NAME: - linearEOM - PURPOSE: - the one-dimensional equation-of-motion - INPUT: - y - current phase-space position - t - current time - pot - (list of) linearPotential instance(s) - OUTPUT: - dy/dt - HISTORY: - 2010-07-13 - Bovy (NYU) - """ - return [y[1],_evaluatelinearForces(pot,y[0],t=t)] diff --git a/galpy/orbit/planarOrbit.py b/galpy/orbit/planarOrbit.py deleted file mode 100644 index ead8d4acb..000000000 --- a/galpy/orbit/planarOrbit.py +++ /dev/null @@ -1,855 +0,0 @@ -import math as m -import warnings -import numpy as nu -from scipy import integrate -from galpy.util.leung_dop853 import dop853 -import galpy.util.bovy_symplecticode as symplecticode -from galpy.util.bovy_conversion import physical_conversion -from .OrbitTop import OrbitTop -from galpy.potential.planarPotential import _evaluateplanarRforces,\ - RZToplanarPotential, toPlanarPotential, _evaluateplanarphiforces,\ - _evaluateplanarPotentials -from galpy.potential.Potential import Potential, _check_c -from galpy.util import galpyWarning, galpyWarningVerbose -#try: -from .integratePlanarOrbit import integratePlanarOrbit_c,\ - integratePlanarOrbit_dxdv_c, _ext_loaded -ext_loaded= _ext_loaded -class planarOrbitTop(OrbitTop): - """Top-level class representing a planar orbit (i.e., one in the plane - of a galaxy)""" - def __init__(self,vxvv=None,vo=220.,ro=8.0,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): #pragma: no cover (never used) - """ - NAME: - - __init__ - - PURPOSE: - - Initialize a planar orbit - - INPUT: - - vxvv - [R,vR,vT(,phi)] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - HISTORY: - - 2010-07-12 - Written - Bovy (NYU) - - 2014-06-11 - Added conversion kwargs to physical coordinates - Bovy (IAS) - - """ - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=zo,vo=vo,solarmotion=solarmotion) - return None - - def e(self,analytic=False,pot=None): - """ - NAME: - e - PURPOSE: - calculate the eccentricity - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - eccentricity - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,type='adiabatic') - (rperi,rap)= self._aA.calcRapRperi(self) - return (rap-rperi)/(rap+rperi) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate eccentricity") - if not hasattr(self,'rs'): - self.rs= self.orbit[:,0] - return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs)) - - @physical_conversion('energy') - def Jacobi(self,*args,**kwargs): - """ - NAME: - Jacobi - PURPOSE: - calculate the Jacobi integral of the motion - INPUT: - t - (optional) time at which to get the radius - OmegaP= pattern speed of rotating frame - pot= potential instance or list of such instances - OUTPUT: - Jacobi integral - HISTORY: - 2011-04-18 - Written - Bovy (NYU) - """ - if not 'OmegaP' in kwargs or kwargs['OmegaP'] is None: - OmegaP= 1. - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - else: - pot= kwargs['pot'] - if isinstance(pot,list): - for p in pot: - if hasattr(p,'OmegaP'): - OmegaP= p.OmegaP() - break - else: - if hasattr(pot,'OmegaP'): - OmegaP= pot.OmegaP() - kwargs.pop('OmegaP',None) - else: - OmegaP= kwargs.pop('OmegaP') - #Make sure you are not using physical coordinates - old_physical= kwargs.get('use_physical',None) - kwargs['use_physical']= False - out= self.E(*args,**kwargs)-OmegaP*self.L(*args,**kwargs) - if not old_physical is None: - kwargs['use_physical']= old_physical - else: - kwargs.pop('use_physical') - return out - - @physical_conversion('position') - def rap(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rap - PURPOSE: - return the apocenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_ap - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,type='adiabatic') - (rperi,rap)= self._aA.calcRapRperi(self) - return rap - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first") - if not hasattr(self,'rs'): - self.rs= self.orbit[:,0] - return nu.amax(self.rs) - - @physical_conversion('position') - def rperi(self,analytic=False,pot=None,**kwargs): - """ - NAME: - rperi - PURPOSE: - return the pericenter radius - INPUT: - analytic - compute this analytically - pot - potential to use for analytical calculation - OUTPUT: - R_peri - HISTORY: - 2010-09-20 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,type='adiabatic') - (rperi,rap)= self._aA.calcRapRperi(self) - return rperi - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first") - if not hasattr(self,'rs'): - self.rs= self.orbit[:,0] - return nu.amin(self.rs) - - @physical_conversion('position') - def zmax(self,pot=None,analytic=False,**kwargs): - raise AttributeError("planarOrbit does not have a zmax") - -class planarROrbit(planarOrbitTop): - """Class representing a planar orbit, without \phi. Useful for - orbit-integration in axisymmetric potentials when you don't care about the - azimuth""" - def __init__(self,vxvv=[1.,0.,1.],vo=220.,ro=8.0,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): - """ - NAME: - - __init__ - - PURPOSE: - - Initialize a planarROrbit - - INPUT: - - vxvv - [R,vR,vT] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - HISTORY: - - 2010-07-12 - Written - Bovy (NYU) - - 2014-06-11 - Added conversion kwargs to physical coordinates - Bovy (IAS) - - """ - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=zo,vo=vo,solarmotion=solarmotion) - return None - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - error message number (get the actual orbit using getOrbit() - HISTORY: - 2010-07-20 - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - thispot= RZToplanarPotential(pot) - self.t= nu.array(t) - self._pot= thispot - self.orbit, msg= _integrateROrbit(self.vxvv,thispot,t,method,dt) - return msg - - @physical_conversion('energy') - def E(self,*args,**kwargs): - """ - NAME: - E - PURPOSE: - calculate the energy - INPUT: - t - (optional) time at which to get the radius - pot= potential instance or list of such instances - OUTPUT: - energy - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - 2011-04-18 - Added t - Bovy (NYU) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if isinstance(pot,Potential): - thispot= RZToplanarPotential(pot) - elif isinstance(pot,list): - thispot= [] - for p in pot: - if isinstance(p,Potential): thispot.append(RZToplanarPotential(p)) - else: thispot.append(p) - else: - thispot= pot - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return _evaluateplanarPotentials(thispot,thiso[0], - t=t)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2. - else: - return nu.array([_evaluateplanarPotentials(thispot,thiso[0,ii], - t=t[ii])\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2. for ii in range(len(t))]) - -class planarOrbit(planarOrbitTop): - """Class representing a full planar orbit (R,vR,vT,phi)""" - def __init__(self,vxvv=[1.,0.,1.,0.],vo=220.,ro=8.0,zo=0.025, - solarmotion=nu.array([-10.1,4.0,6.7])): - """ - NAME: - - __init__ - - PURPOSE: - - Initialize a planarOrbit - - INPUT: - - vxvv - [R,vR,vT,phi] - - vo - circular velocity at ro (km/s) - - ro - distance from vantage point to GC (kpc) - - zo - offset toward the NGP of the Sun wrt the plane (kpc) - - solarmotion - value in [-U,V,W] (km/s) - - OUTPUT: - - HISTORY: - - 2010-07-12 - Written - Bovy (NYU) - - 2014-06-11 - Added conversion kwargs to physical coordinates - Bovy (IAS) - - """ - if len(vxvv) == 3: #pragma: no cover - raise ValueError("You only provided R,vR, & vT, but not phi; you probably want planarROrbit") - OrbitTop.__init__(self,vxvv=vxvv, - ro=ro,zo=zo,vo=vo,solarmotion=solarmotion) - return None - - def integrate(self,t,pot,method='symplec4_c',dt=None): - """ - NAME: - integrate - PURPOSE: - integrate the orbit - INPUT: - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - dt= (None) if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - (none) (get the actual orbit using getOrbit() - HISTORY: - 2010-07-20 - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - thispot= toPlanarPotential(pot) - self.t= nu.array(t) - self._pot= thispot - self.orbit, msg= _integrateOrbit(self.vxvv,thispot,t,method,dt) - return msg - - def integrate_dxdv(self,dxdv,t,pot,method='dopr54_c', - rectIn=False,rectOut=False): - """ - NAME: - integrate_dxdv - PURPOSE: - integrate the orbit and a small area of phase space - INPUT: - dxdv - [dR,dvR,dvT,dphi] - t - list of times at which to output (0 has to be in this!) - pot - potential instance or list of instances - method= 'odeint' for scipy's odeint - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a Dormand-Prince integrator in C (generally the fastest) - rectIn= (False) if True, input dxdv is in rectangular coordinates - rectOut= (False) if True, output dxdv (that in orbit_dxdv) is in rectangular coordinates - OUTPUT: - (none) (get the actual orbit using getOrbit_dxdv() - HISTORY: - 2010-10-17 - Written - Bovy (NYU) - 2014-06-29 - Added rectIn and rectOut - Bovy (IAS) - """ - if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') - if hasattr(self,'rs'): delattr(self,'rs') - thispot= toPlanarPotential(pot) - self.t= nu.array(t) - self._pot_dxdv= thispot - self._pot= thispot - self.orbit_dxdv, msg= _integrateOrbit_dxdv(self.vxvv,dxdv,thispot,t, - method,rectIn,rectOut) - self.orbit= self.orbit_dxdv[:,:4] - return msg - - @physical_conversion('energy') - def E(self,*args,**kwargs): - """ - NAME: - E - PURPOSE: - calculate the energy - INPUT: - pot= - t= time at which to evaluate E - OUTPUT: - energy - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if not 'pot' in kwargs or kwargs['pot'] is None: - try: - pot= self._pot - except AttributeError: - raise AttributeError("Integrate orbit or specify pot=") - if 'pot' in kwargs and kwargs['pot'] is None: - kwargs.pop('pot') - else: - pot= kwargs.pop('pot') - if isinstance(pot,Potential): - thispot= toPlanarPotential(pot) - elif isinstance(pot,list): - thispot= [] - for p in pot: - if isinstance(p,Potential): thispot.append(toPlanarPotential(p)) - else: thispot.append(p) - else: - thispot= pot - if len(args) > 0: - t= args[0] - else: - t= 0. - #Get orbit - thiso= self(*args,**kwargs) - onet= (len(thiso.shape) == 1) - if onet: - return _evaluateplanarPotentials(thispot,thiso[0], - phi=thiso[3],t=t)\ - +thiso[1]**2./2.\ - +thiso[2]**2./2. - else: - return nu.array([_evaluateplanarPotentials(thispot,thiso[0,ii], - phi=thiso[3,ii], - t=t[ii])\ - +thiso[1,ii]**2./2.\ - +thiso[2,ii]**2./2. for ii in range(len(t))]) - - def e(self,analytic=False,pot=None): - """ - NAME: - e - PURPOSE: - calculate the eccentricity - INPUT: - analytic - calculate e analytically - pot - potential used to analytically calculate e - OUTPUT: - eccentricity - HISTORY: - 2010-09-15 - Written - Bovy (NYU) - """ - if analytic: - self._setupaA(pot=pot,type='adiabatic') - (rperi,rap)= self._aA.calcRapRperi(self) - return (rap-rperi)/(rap+rperi) - if not hasattr(self,'orbit'): - raise AttributeError("Integrate the orbit first or use analytic=True for approximate eccentricity") - if not hasattr(self,'rs'): - self.rs= self.orbit[:,0] - return (nu.amax(self.rs)-nu.amin(self.rs))/(nu.amax(self.rs)+nu.amin(self.rs)) - -def _integrateROrbit(vxvv,pot,t,method,dt): - """ - NAME: - _integrateROrbit - PURPOSE: - integrate an orbit in a Phi(R) potential in the R-plane - INPUT: - vxvv - array with the initial conditions stacked like - [R,vR,vT]; vR outward! - pot - Potential instance - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - dt - if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - [:,3] array of [R,vR,vT] at each t - HISTORY: - 2010-07-20 - Written - Bovy (NYU) - """ - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - if method.lower() == 'leapfrog': - #We hack this by putting in a dummy phi - this_vxvv= nu.zeros(len(vxvv)+1) - this_vxvv[0:len(vxvv)]= vxvv - tmp_out, msg= _integrateOrbit(this_vxvv,pot,t,method,dt) - #tmp_out is (nt,4) - out= tmp_out[:,0:3] - elif ext_loaded and \ - (method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \ - or method.lower() == 'rk6_c' or method.lower() == 'symplec4_c' \ - or method.lower() == 'symplec6_c' or method.lower() == 'dopr54_c' or method.lower() == 'dop853_c'): - #We hack this by putting in a dummy phi - this_vxvv= nu.zeros(len(vxvv)+1) - this_vxvv[0:len(vxvv)]= vxvv - tmp_out, msg= _integrateOrbit(this_vxvv,pot,t,method,dt) - #tmp_out is (nt,4) - out= tmp_out[:,0:3] - elif method.lower() == 'odeint' or method.lower() == 'dop853' or not ext_loaded: - l= vxvv[0]*vxvv[2] - l2= l**2. - init= [vxvv[0],vxvv[1]] - if method.lower() == 'dop853': - intOut = dop853(_REOM, init, t, args=(pot, l2)) - else: - intOut= integrate.odeint(_REOM,init,t,args=(pot,l2), - rtol=10.**-8.)#,mxstep=100000000) - out= nu.zeros((len(t),3)) - out[:,0]= intOut[:,0] - out[:,1]= intOut[:,1] - out[:,2]= l/out[:,0] - msg= 0 - #post-process to remove negative radii - neg_radii= (out[:,0] < 0.) - out[neg_radii,0]= -out[neg_radii,0] - _parse_warnmessage(msg) - return (out,msg) - -def _REOM(y,t,pot,l2): - """ - NAME: - _REOM - PURPOSE: - implements the EOM, i.e., the right-hand side of the differential - equation - INPUT: - y - current phase-space position - t - current time - pot - (list of) Potential instance(s) - l2 - angular momentum squared - OUTPUT: - dy/dt - HISTORY: - 2010-07-20 - Written - Bovy (NYU) - """ - return [y[1], - l2/y[0]**3.+_evaluateplanarRforces(pot,y[0],t=t)] - -def _integrateOrbit(vxvv,pot,t,method,dt): - """ - NAME: - _integrateOrbit - PURPOSE: - integrate an orbit in a Phi(R) potential in the (R,phi)-plane - INPUT: - vxvv - array with the initial conditions stacked like - [R,vR,vT,phi]; vR outward! - pot - Potential instance - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - dt- if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize - OUTPUT: - [:,4] array of [R,vR,vT,phi] at each t - HISTORY: - 2010-07-20 - Written - Bovy (NYU) - """ - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - if method.lower() == 'leapfrog': - #go to the rectangular frame - this_vxvv= nu.array([vxvv[0]*nu.cos(vxvv[3]), - vxvv[0]*nu.sin(vxvv[3]), - vxvv[1]*nu.cos(vxvv[3])-vxvv[2]*nu.sin(vxvv[3]), - vxvv[2]*nu.cos(vxvv[3])+vxvv[1]*nu.sin(vxvv[3])]) - #integrate - tmp_out= symplecticode.leapfrog(_rectForce,this_vxvv, - t,args=(pot,),rtol=10.**-8) - #go back to the cylindrical frame - R= nu.sqrt(tmp_out[:,0]**2.+tmp_out[:,1]**2.) - phi= nu.arccos(tmp_out[:,0]/R) - phi[(tmp_out[:,1] < 0.)]= 2.*nu.pi-phi[(tmp_out[:,1] < 0.)] - vR= tmp_out[:,2]*nu.cos(phi)+tmp_out[:,3]*nu.sin(phi) - vT= tmp_out[:,3]*nu.cos(phi)-tmp_out[:,2]*nu.sin(phi) - out= nu.zeros((len(t),4)) - out[:,0]= R - out[:,1]= vR - out[:,2]= vT - out[:,3]= phi - msg= 0 - elif ext_loaded and \ - (method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \ - or method.lower() == 'rk6_c' or method.lower() == 'symplec4_c' \ - or method.lower() == 'symplec6_c' or method.lower() == 'dopr54_c' \ - or method.lower() == 'dop853_c'): - warnings.warn("Using C implementation to integrate orbits",galpyWarningVerbose) - out, msg= integratePlanarOrbit_c(pot,nu.copy(vxvv),t,method,dt=dt) - elif method.lower() == 'odeint' or method.lower() == 'dop853' or not ext_loaded: - vphi= vxvv[2]/vxvv[0] - init= [vxvv[0],vxvv[1],vxvv[3],vphi] - if method == 'dop853': - intOut = dop853(_EOM, init, t, args=(pot,)) - else: - intOut= integrate.odeint(_EOM,init,t,args=(pot,), - rtol=10.**-8.)#,mxstep=100000000) - out= nu.zeros((len(t),4)) - out[:,0]= intOut[:,0] - out[:,1]= intOut[:,1] - out[:,3]= intOut[:,2] - out[:,2]= out[:,0]*intOut[:,3] - msg= 0 - else: - raise NotImplementedError("requested integration method does not exist") - #post-process to remove negative radii - neg_radii= (out[:,0] < 0.) - out[neg_radii,0]= -out[neg_radii,0] - out[neg_radii,3]+= m.pi - _parse_warnmessage(msg) - return (out,msg) - -def _integrateOrbit_dxdv(vxvv,dxdv,pot,t,method,rectIn,rectOut): - """ - NAME: - _integrateOrbit_dxdv - PURPOSE: - integrate an orbit and area of phase space in a Phi(R) potential - in the (R,phi)-plane - INPUT: - vxvv - array with the initial conditions stacked like - [R,vR,vT,phi]; vR outward! - dxdv - difference to integrate [dR,dvR,dvT,dphi] - pot - Potential instance - t - list of times at which to output (0 has to be in this!) - method - 'odeint' or 'leapfrog' - rectIn= (False) if True, input dxdv is in rectangular coordinates - rectOut= (False) if True, output dxdv (that in orbit_dxdv) is in rectangular coordinates - OUTPUT: - [:,8] array of [R,vR,vT,phi,dR,dvR,dvT,dphi] at each t - error message from integrator - HISTORY: - 2010-10-17 - Written - Bovy (IAS) - """ - #First check that the potential has C - if '_c' in method: - allHasC= _check_c(pot) and _check_c(pot,dxdv=True) - if not ext_loaded or \ - (not allHasC and not 'leapfrog' in method and not 'symplec' in method): - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Using odeint because C extension not loaded",galpyWarning) - else: - warnings.warn("Using odeint because not all used potential have adequate C implementations to integrate phase-space volumes",galpyWarning) - #go to the rectangular frame - this_vxvv= nu.array([vxvv[0]*nu.cos(vxvv[3]), - vxvv[0]*nu.sin(vxvv[3]), - vxvv[1]*nu.cos(vxvv[3])-vxvv[2]*nu.sin(vxvv[3]), - vxvv[2]*nu.cos(vxvv[3])+vxvv[1]*nu.sin(vxvv[3])]) - if not rectIn: - this_dxdv= nu.array([nu.cos(vxvv[3])*dxdv[0] - -vxvv[0]*nu.sin(vxvv[3])*dxdv[3], - nu.sin(vxvv[3])*dxdv[0] - +vxvv[0]*nu.cos(vxvv[3])*dxdv[3], - -(vxvv[1]*nu.sin(vxvv[3]) - +vxvv[2]*nu.cos(vxvv[3]))*dxdv[3] - +nu.cos(vxvv[3])*dxdv[1]-nu.sin(vxvv[3])*dxdv[2], - (vxvv[1]*nu.cos(vxvv[3]) - -vxvv[2]*nu.sin(vxvv[3]))*dxdv[3] - +nu.sin(vxvv[3])*dxdv[1]+nu.cos(vxvv[3])*dxdv[2]]) - else: - this_dxdv= dxdv - if 'leapfrog' in method.lower() or 'symplec' in method.lower(): - raise TypeError('Symplectic integration for phase-space volume is not possible') - elif ext_loaded and \ - (method.lower() == 'rk4_c' or method.lower() == 'rk6_c' \ - or method.lower() == 'dopr54_c' or method.lower() == 'dop853_c'): - warnings.warn("Using C implementation to integrate orbits",galpyWarningVerbose) - #integrate - tmp_out, msg= integratePlanarOrbit_dxdv_c(pot,this_vxvv,this_dxdv, - t,method) - elif method.lower() == 'odeint' or not ext_loaded or method.lower() == 'dop853': - init= [this_vxvv[0],this_vxvv[1],this_vxvv[2],this_vxvv[3], - this_dxdv[0],this_dxdv[1],this_dxdv[2],this_dxdv[3]] - #integrate - if method.lower() == "dop853": - tmp_out = dop853(_EOM_dxdv, init, t, args=(pot,)) - else: - tmp_out= integrate.odeint(_EOM_dxdv,init,t,args=(pot,), - rtol=10.**-8.)#,mxstep=100000000) - msg= 0 - else: - raise NotImplementedError("requested integration method does not exist") - #go back to the cylindrical frame - R= nu.sqrt(tmp_out[:,0]**2.+tmp_out[:,1]**2.) - phi= nu.arccos(tmp_out[:,0]/R) - phi[(tmp_out[:,1] < 0.)]= 2.*nu.pi-phi[(tmp_out[:,1] < 0.)] - vR= tmp_out[:,2]*nu.cos(phi)+tmp_out[:,3]*nu.sin(phi) - vT= tmp_out[:,3]*nu.cos(phi)-tmp_out[:,2]*nu.sin(phi) - cp= nu.cos(phi) - sp= nu.sin(phi) - dR= cp*tmp_out[:,4]+sp*tmp_out[:,5] - dphi= (cp*tmp_out[:,5]-sp*tmp_out[:,4])/R - dvR= cp*tmp_out[:,6]+sp*tmp_out[:,7]+vT*dphi - dvT= cp*tmp_out[:,7]-sp*tmp_out[:,6]-vR*dphi - out= nu.zeros((len(t),8)) - out[:,0]= R - out[:,1]= vR - out[:,2]= vT - out[:,3]= phi - if rectOut: - out[:,4:]= tmp_out[:,4:] - else: - out[:,4]= dR - out[:,7]= dphi - out[:,5]= dvR - out[:,6]= dvT - _parse_warnmessage(msg) - return (out,msg) - -def _EOM_dxdv(x,t,pot): - """ - NAME: - _EOM_dxdv - PURPOSE: - implements the EOM, i.e., the right-hand side of the differential - equation, for integrating phase space differences, rectangular - INPUT: - x - current phase-space position - t - current time - pot - (list of) Potential instance(s) - OUTPUT: - dy/dt - HISTORY: - 2011-10-18 - Written - Bovy (NYU) - """ - #x is rectangular so calculate R and phi - R= nu.sqrt(x[0]**2.+x[1]**2.) - phi= nu.arccos(x[0]/R) - sinphi= x[1]/R - cosphi= x[0]/R - if x[1] < 0.: phi= 2.*nu.pi-phi - #calculate forces - Rforce= _evaluateplanarRforces(pot,R,phi=phi,t=t) - phiforce= _evaluateplanarphiforces(pot,R,phi=phi,t=t) - R2deriv= _evaluateplanarPotentials(pot,R,phi=phi,t=t,dR=2) - phi2deriv= _evaluateplanarPotentials(pot,R,phi=phi,t=t,dphi=2) - Rphideriv= _evaluateplanarPotentials(pot,R,phi=phi,t=t,dR=1,dphi=1) - #Calculate derivatives and derivatives+time derivatives - dFxdx= -cosphi**2.*R2deriv\ - +2.*cosphi*sinphi/R**2.*phiforce\ - +sinphi**2./R*Rforce\ - +2.*sinphi*cosphi/R*Rphideriv\ - -sinphi**2./R**2.*phi2deriv - dFxdy= -sinphi*cosphi*R2deriv\ - +(sinphi**2.-cosphi**2.)/R**2.*phiforce\ - -cosphi*sinphi/R*Rforce\ - -(cosphi**2.-sinphi**2.)/R*Rphideriv\ - +cosphi*sinphi/R**2.*phi2deriv - dFydx= -cosphi*sinphi*R2deriv\ - +(sinphi**2.-cosphi**2.)/R**2.*phiforce\ - +(sinphi**2.-cosphi**2.)/R*Rphideriv\ - -sinphi*cosphi/R*Rforce\ - +sinphi*cosphi/R**2.*phi2deriv - dFydy= -sinphi**2.*R2deriv\ - -2.*sinphi*cosphi/R**2.*phiforce\ - -2.*sinphi*cosphi/R*Rphideriv\ - +cosphi**2./R*Rforce\ - -cosphi**2./R**2.*phi2deriv - return nu.array([x[2],x[3], - cosphi*Rforce-1./R*sinphi*phiforce, - sinphi*Rforce+1./R*cosphi*phiforce, - x[6],x[7], - dFxdx*x[4]+dFxdy*x[5], - dFydx*x[4]+dFydy*x[5]]) - -def _EOM(y,t,pot): - """ - NAME: - _EOM - PURPOSE: - implements the EOM, i.e., the right-hand side of the differential - equation - INPUT: - y - current phase-space position - t - current time - pot - (list of) Potential instance(s) - l2 - angular momentum squared - OUTPUT: - dy/dt - HISTORY: - 2010-07-20 - Written - Bovy (NYU) - """ - l2= (y[0]**2.*y[3])**2. - return [y[1], - l2/y[0]**3.+_evaluateplanarRforces(pot,y[0],phi=y[2],t=t), - y[3], - 1./y[0]**2.*(_evaluateplanarphiforces(pot,y[0],phi=y[2],t=t)- - 2.*y[0]*y[1]*y[3])] - -def _rectForce(x,pot,t=0.): - """ - NAME: - _rectForce - PURPOSE: - returns the force in the rectangular frame - INPUT: - x - current position - t - current time - pot - (list of) Potential instance(s) - OUTPUT: - force - HISTORY: - 2011-02-02 - Written - Bovy (NYU) - """ - #x is rectangular so calculate R and phi - R= nu.sqrt(x[0]**2.+x[1]**2.) - phi= nu.arccos(x[0]/R) - sinphi= x[1]/R - cosphi= x[0]/R - if x[1] < 0.: phi= 2.*nu.pi-phi - #calculate forces - Rforce= _evaluateplanarRforces(pot,R,phi=phi,t=t) - phiforce= _evaluateplanarphiforces(pot,R,phi=phi,t=t) - return nu.array([cosphi*Rforce-1./R*sinphi*phiforce, - sinphi*Rforce+1./R*cosphi*phiforce]) - -def _parse_warnmessage(msg): - if msg == 1: #pragma: no cover - warnings.warn("During numerical integration, steps smaller than the smallest step were requested; integration might not be accurate",galpyWarning) -