diff --git a/CHANGES.rst b/CHANGES.rst index bd12bd68..464dd76b 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -11,6 +11,8 @@ New Features objects, or to create Gala potential objects from a pre-existing Galpy potential. +- Added a ``plot_3d()`` method for ``Orbit`` objects to make 3D plots of the + orbital trajectories. Bug fixes --------- @@ -33,6 +35,8 @@ Bug fixes API changes ----------- +- Removed the deprecated ``gala.coordinates.get_galactocentric2019()`` function. + 1.3 (2020-10-27) ================ diff --git a/docs/conventions.rst b/docs/conventions.rst index c0b4b009..fc430d23 100644 --- a/docs/conventions.rst +++ b/docs/conventions.rst @@ -23,16 +23,16 @@ Array shapes ============ The arrays and :class:`~astropy.units.Quantity` objects expected as input and -returned as ouput throughout ``Gala`` have shapes that follow a particular +returned as output throughout ``Gala`` have shapes that follow a particular convention, unless otherwise specified in function or class docstrings. For arrays containing coordinate or kinematic information, ``axis=0`` is assumed -to be the coordinate dimension. For example, for representing 128 Cartesian 3D -positions, the object would have shape ``(3,128)``. +to be the coordinate dimension. For example, for representing 128 different 3D +Cartesian positions, the object should have shape ``(3, 128)``. -For collections of orbits, arrays have 3 axes. As above, ``axis=0`` is assumed -to be the coordinate dimension, but now ``axis=1`` is interpreted as the time -axis, and ``axis=2`` are the individual orbits. +For collections of orbits, arrays or array-valued objects have three axes: As +above, ``axis=0`` is assumed to be the coordinate dimension, ``axis=1`` is +interpreted as the time axis, and ``axis=2`` are the different orbits. .. _energy-momentum: @@ -40,6 +40,6 @@ Energy and momentum =================== The `gala` documentation and functions often refer to energy and angular -momentum and the respective quantities *per unit mass* interchangibly. Unless +momentum and the respective quantities *per unit mass* interchangeably. Unless otherwise specified, all such quantities -- e.g., energy, angular momentum, momenta, conjugate momenta -- are in fact used and returned *per unit mass*. diff --git a/docs/coordinates/greatcircle.rst b/docs/coordinates/greatcircle.rst index b329e220..bd20990d 100644 --- a/docs/coordinates/greatcircle.rst +++ b/docs/coordinates/greatcircle.rst @@ -21,14 +21,14 @@ coordinate system, such as the ICRS. The great circle system is defined by a specified north pole, with a additional (optional) specification of the longitude zero point of the final system. -Currently, ``gala`` supports great circle frames that are defined as a rotation +`gala` currentlt supports great circle frames that are defined as a rotation away from the ICRS (ra, dec) through the `~gala.coordinates.GreatCircleICRSFrame` class. To create a new great circle -frame, you must specify a pole using the ``pole=`` keywrod, and optionally +frame, you must specify a pole using the ``pole=`` keyword, and optionally specify the longitude zero point either by specifying the right ascension of the longitude zero point, ``ra0``, or by specifying a final rotation to be applied to the transformation, ``rotation``. For example, to define a great circle -system with the pole at (RA, Dec) = (32.5, 19.8), we first have to create a +system with the pole at (RA, Dec) = (32.5, 19.8)ยบ, we first have to create a coordinate object for the pole:: >>> pole = coord.SkyCoord(ra=32.5*u.deg, dec=19.8*u.deg) @@ -39,7 +39,7 @@ define our coordinate frame:: >>> fr = gc.GreatCircleICRSFrame(pole=pole) We can then use this frame like any other Astropy coordinate frame. For example, -we can transform other coordinate to this new coordinate system using:: +we can transform other coordinates to this new coordinate system using:: >>> c = coord.SkyCoord(ra=[160, 53]*u.deg, dec=[-11, 9]*u.deg) >>> c_fr = c.transform_to(fr) @@ -61,7 +61,8 @@ The transformation also works for velocity components. For example, if we have a sky position and proper motions, we can transform to the great circle frame in the same way:: - >>> c = coord.SkyCoord(ra=160*u.deg, dec=-11*u.deg, + >>> c = coord.SkyCoord(ra=160*u.deg, + ... dec=-11*u.deg, ... pm_ra_cosdec=5*u.mas/u.yr, ... pm_dec=0.3*u.mas/u.yr) >>> c_fr = c.transform_to(fr) diff --git a/docs/coordinates/index.rst b/docs/coordinates/index.rst index 492fadaa..b3a0d959 100644 --- a/docs/coordinates/index.rst +++ b/docs/coordinates/index.rst @@ -30,7 +30,7 @@ values adopted in Astropy v4.0: Stellar stream coordinate frames ================================ -`Gala` provides Astropy coordinate frame classes for transforming to several +`gala` provides Astropy coordinate frame classes for transforming to several built-in stellar stream stream coordinate frames (as defined in the references below), and for transforming positions and velocities to and from coordinate systems defined by great circles or poles. These classes behave like the @@ -42,7 +42,7 @@ with the Sagittarius stream with the `~gala.coordinates.SagittariusLaw10` frame:: >>> c = coord.ICRS(ra=100.68458*u.degree, dec=41.26917*u.degree) - >>> sgr = c.transform_to(gc.SagittariusLaw10) + >>> sgr = c.transform_to(gc.SagittariusLaw10()) >>> (sgr.Lambda, sgr.Beta) # doctest: +FLOAT_CMP (, ) @@ -50,7 +50,7 @@ Or, to transform from `~gala.coordinates.SagittariusLaw10` coordinates to the `~astropy.coordinates.Galactic` frame:: >>> sgr = gc.SagittariusLaw10(Lambda=156.342*u.degree, Beta=1.1*u.degree) - >>> c = sgr.transform_to(coord.Galactic) + >>> c = sgr.transform_to(coord.Galactic()) >>> (c.l, c.b) # doctest: +FLOAT_CMP (, ) @@ -59,10 +59,10 @@ can be transformed between the systems. For example, to transform from `~gala.coordinates.GD1Koposov10` proper motions to `~astropy.coordinates.Galactic` proper motions:: - >>> gd1 = gc.GD1Koposov10(phi1=-35.00*u.degree, phi2=0*u.degree, + >>> gd1 = gc.GD1Koposov10(phi1=-35*u.degree, phi2=0*u.degree, ... pm_phi1_cosphi2=-12.20*u.mas/u.yr, ... pm_phi2=-3.10*u.mas/u.yr) - >>> gd1.transform_to(coord.Galactic) # doctest: +FLOAT_CMP + >>> gd1.transform_to(coord.Galactic()) # doctest: +FLOAT_CMP >> gd1.transform_to(coord.Galactocentric) # doctest: +FLOAT_CMP + >>> gd1.transform_to(coord.Galactocentric()) # doctest: +FLOAT_CMP , galcen_distance=8.122 kpc, galcen_v_sun=(12.9, 245.6, 7.78) km / s, z_sun=20.8 pc, roll=0.0 deg): (x, y, z) in kpc (-12.61622659, -0.09870921, 6.43179403) @@ -97,7 +97,8 @@ and returns a coordinate object with the solar motion added back in to the velocity components. This is useful for computing velocities in a Galactocentric reference frame, rather than a solar system barycentric frame. -To use, you pass in a coordinate object with scalar or array values:: +The `~gala.coordinates.reflex_correct` function accepts a coordinate object with +scalar or array values:: >>> c = coord.SkyCoord(ra=[180.323, 1.523]*u.deg, ... dec=[-17, 29]*u.deg, @@ -112,10 +113,12 @@ To use, you pass in a coordinate object with scalar or array values:: [(139.47001884, 175.45769809, -47.09032586), (-61.01738781, 61.51055793, 163.36721898)]> -By default, this uses the default solar location and velocity from the -`astropy.coordinates.Galactocentric` frame class. To modify these, for example, -to change the solar velocity, or the sun's height above the Galactic midplane, -use the arguments of that class:: +By default, this uses the solar location and velocity from the +`astropy.coordinates.Galactocentric` frame class. To modify these parameters, +for example, to change the solar velocity, or the sun's height above the +Galactic midplane, use the arguments of the `astropy.coordinates.Galactocentric` +class and pass in an instance of the `astropy.coordinates.Galactocentric` +frame:: >>> vsun = coord.CartesianDifferential([11., 245., 7.]*u.km/u.s) >>> gc_frame = coord.Galactocentric(galcen_v_sun=vsun, z_sun=0*u.pc) @@ -162,19 +165,18 @@ values of the proper motions) -- this is sometimes called "v_GSR":: Transforming a proper motion covariance matrix to a new coordinate frame ------------------------------------------------------------------------ -When working with Gaia or other astrometric data sets, we often need to -transform the reported covariance matrix between proper motion components into a -new coordinate system. For example, Gaia data are provided in the +When working with Gaia or other astrometric data sets, you may need to transform +the reported covariance matrix between proper motion components into a new +coordinate system. For example, Gaia data are provided in the `~astropy.coordinates.ICRS` (equatorial) coordinate frame, but for Galactic -science, we often want to instead work in the -`~astropy.coordinates.Galactic` coordinate system. -For this and other transformations that only require a rotation (i.e. the origin -doesn't change), the astrometric covariance matrix can be transformed exactly -through a projection of the rotation onto the tangent plane at a given location. -The details of this procedure are explained in `this document from the Gaia data -processing team +science, we often want to instead work in the `~astropy.coordinates.Galactic` +coordinate system. For this and other transformations that only require a +rotation (i.e. the origin doesn't change), the astrometric covariance matrix can +be transformed exactly through a projection of the rotation onto the tangent +plane at a given location. The details of this procedure are explained in `this +document from the Gaia data processing team `_, -and this functionality is implemented in gala. Let's first create a coordinate +and this functionality is implemented in `gala`. Let's first create a coordinate object to transform:: >>> c = coord.SkyCoord(ra=62*u.deg, @@ -190,23 +192,22 @@ be constructed from a single row from a Gaia data release source catalog:: This matrix specifies the 2D error distribution for the proper motion measurement *in the ICRS frame*. To transform this matrix to, e.g., the Galactic -coordinate system, we can use the function +coordinate system, use the function `~gala.coordinates.transform_pm_cov`:: - >>> gc.transform_pm_cov(c, cov, coord.Galactic) # doctest: +FLOAT_CMP + >>> gc.transform_pm_cov(c, cov, coord.Galactic()) # doctest: +FLOAT_CMP array([[ 0.69450047, -0.309945 ], [-0.309945 , 0.96413005]]) -Note that this works for all of the great circle or stellar stream coordinate -frames implemented in gala:: +Note that this also works for all of the great circle or stellar stream +coordinate frames implemented in `gala`:: - >>> gc.transform_pm_cov(c, cov, gc.GD1Koposov10) # doctest: +FLOAT_CMP + >>> gc.transform_pm_cov(c, cov, gc.GD1Koposov10()) # doctest: +FLOAT_CMP array([[1.10838914, 0.19067958], [0.19067958, 0.55024138]]) -Also note that this works for multiple coordinates (i.e. array coordinates) as -well, so try to avoid looping over this function and instead apply it to -array-valued coordinate objects. +This works for array-valued coordinates as well, so try to avoid looping over +this function and instead apply it to array-valued coordinate objects. References diff --git a/docs/dynamics/actionangle.rst b/docs/dynamics/actionangle.rst index 74c11145..b6c12fab 100644 --- a/docs/dynamics/actionangle.rst +++ b/docs/dynamics/actionangle.rst @@ -89,6 +89,7 @@ We will now integrate the orbit and plot it in the meridional plane:: .. plot:: :align: center :context: close-figs + :width: 60% import astropy.coordinates as coord import astropy.units as u @@ -110,10 +111,10 @@ We will now integrate the orbit and plot it in the meridional plane:: To solve for the actions in the true potential, we first compute the actions in a "toy" potential -- a potential in which we can compute the actions and angles analytically. The two simplest potentials for which this is possible are the -`~gala.potential.IsochronePotential` and -`~gala.potential.HarmonicOscillatorPotential`. We will use the Isochrone -potential as our toy potential for tube orbits and the harmonic oscillator for -box orbits. +`~gala.potential.potential.IsochronePotential` and +`~gala.potential.potential.HarmonicOscillatorPotential`. We will use the +Isochrone potential as our toy potential for tube orbits and the harmonic +oscillator for box orbits. We start by finding the parameters of the toy potential (Isochrone in this case) by minimizing the dispersion in energy for the orbit:: @@ -138,6 +139,7 @@ Instead, the orbit is wobbly in the toy potential angles:: .. plot:: :align: center :context: close-figs + :width: 60% toy_potential = gd.fit_isochrone(w) toy_actions,toy_angles,toy_freqs = toy_potential.action_angle(w) @@ -160,6 +162,7 @@ time-independent in the toy potential:: .. plot:: :align: center :context: close-figs + :width: 60% fig,ax = plt.subplots(1,1) ax.plot(w.t, toy_actions[0].to(u.km/u.s*u.kpc), marker='') @@ -168,7 +171,7 @@ time-independent in the toy potential:: fig.tight_layout() We can now find approximations to the actions in the true potential. We have to -choose the maximum integer vector norm, `N_max`, which here we arbitrarilty set +choose the maximum integer vector norm, `N_max`, which here we arbitrarily set to 8. This will change depending on the convergence of the action correction (the properties of the orbit and potential) and the accuracy desired:: @@ -200,6 +203,7 @@ actions computed using this machinery:: .. plot:: :align: center :context: close-figs + :width: 60% import warnings with warnings.catch_warnings(record=True): @@ -252,6 +256,7 @@ and the same initial conditions as above: .. plot:: :align: center :include-source: + :width: 60% import astropy.coordinates as coord import astropy.units as u diff --git a/docs/dynamics/index.rst b/docs/dynamics/index.rst index 1f8d8a72..2a5f049f 100644 --- a/docs/dynamics/index.rst +++ b/docs/dynamics/index.rst @@ -6,8 +6,7 @@ Dynamics (`gala.dynamics`) ******************************** -For code blocks below and any pages linked below, we'll assume the following -imports have already been executed:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np @@ -73,6 +72,7 @@ coordinate representation, for example, cylindrical radius :math:`\rho` and .. plot:: :align: center + :width: 60% import astropy.units as u import gala.potential as gp diff --git a/docs/dynamics/mockstreams.rst b/docs/dynamics/mockstreams.rst index f39ecc33..67a3989e 100644 --- a/docs/dynamics/mockstreams.rst +++ b/docs/dynamics/mockstreams.rst @@ -45,8 +45,8 @@ km/s, and a scale radius of 15 kpc:: The mock stream generation supports any of the reference frames implemented in ``gala`` (e.g., non-static / rotating reference frames), so we must create a -`~gala.potential.Hamiltonian` object to use when generating streams. By default, -this will use a static reference frame:: +`~gala.potential.hamiltonian.Hamiltonian` object to use when generating streams. +By default, this will use a static reference frame:: >>> H = gp.Hamiltonian(pot) @@ -99,6 +99,7 @@ Let's plot the stream:: .. plot:: :align: center :context: close-figs + :width: 60% import astropy.units as u import numpy as np @@ -143,7 +144,7 @@ Hamiltonian. It is possible to include a potential object for the progenitor system to account for the self-gravity of the progenitor as stream star particles are released. We can use any of the ``gala.potential`` potential objects to represent the progenitor system, but here we will use a simple -`~gala.potential.PlummerPotential`. We pass this in to the +`~gala.potential.potential.PlummerPotential`. We pass this in to the `~gala.dynamics.mockstream.MockStreamGenerator` - let's see what the stream looks like when generated including self-gravity:: @@ -156,6 +157,7 @@ looks like when generated including self-gravity:: .. plot:: :align: center :context: close-figs + :width: 60% prog_pot = gp.PlummerPotential(m=prog_mass, b=2*u.pc, units=galactic) gen2 = ms.MockStreamGenerator(df, H, progenitor_potential=prog_pot) diff --git a/docs/dynamics/nbody.rst b/docs/dynamics/nbody.rst index c30a0c23..083491cb 100644 --- a/docs/dynamics/nbody.rst +++ b/docs/dynamics/nbody.rst @@ -7,19 +7,18 @@ N-body (`gala.dynamics.nbody`) Introduction ============ -With the `~gala.potential.Hamiltonian` and potential classes (:ref:`potential`), -Gala contains functionality for integrating test particle orbits in background -gravitational fields. To supplement this, Gala also now contains some limited -functionality for performing N-body orbit integrations through direct N-body -force calculations between particles. With the `gala.dynamics.nbody` subpackage, -gravitational fields (i.e., any potential class from ``gala.potential``) can be -sourced by particles that interact, and optionally feel a background/external -potential. To use this functionality, the core class is -`~gala.dynamics.DirectNBody`. Below, we'll go through a few examples of using -this class to perform orbit integrations - -For code blocks below, we'll assume the following imports have already been -executed:: +With the `~gala.potential.hamiltonian.Hamiltonian` and potential classes +(:ref:`potential`), Gala contains functionality for integrating test particle +orbits in background gravitational fields. To supplement this, Gala also now +contains some limited functionality for performing N-body orbit integrations +through direct N-body force calculations between particles. With the +`gala.dynamics.nbody` subpackage, gravitational fields (i.e., any potential +class from :mod:`gala.potential`) can be sourced by particles that interact, and +optionally feel a background/external potential. To use this functionality, the +core class is `~gala.dynamics.nbody.DirectNBody`. Below, we'll go through a few +examples of using this class to perform orbit integrations + +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np @@ -32,25 +31,25 @@ executed:: Getting started =============== -The `~gala.dynamics.DirectNBody`, at minimum, must be instantiated with a set of -particle orbital initial conditions along with a specification of the +The `~gala.dynamics.nbody.DirectNBody`, at minimum, must be instantiated with a +set of particle orbital initial conditions along with a specification of the gravitational fields sourced by each particle --- that is, the number of initial conditions must match the input list of gravitational potential objects that specify the particle mass distributions. Other optional arguments to -`~gala.dynamics.DirectNBody` allow you to set the unit system (i.e., to improve -numerical precision when time-stepping the orbit integration), or to specify a -background gravitational potential. Let's now go through a few examples of using -this class in practice. +`~gala.dynamics.nbody.DirectNBody` allow you to set the unit system (i.e., to +improve numerical precision when time-stepping the orbit integration), or to +specify a background gravitational potential. Let's now go through a few +examples of using this class in practice. Example: Mixed test particle and massive particle orbit integration =================================================================== -Like with `~gala.potential.Hamiltonian` orbit integration, orbital initial -conditions are passed in to `~gala.dynamics.DirectNBody` by passing in a -single `~gala.dynamics.PhaseSpacePosition` object. Let's create two initial -conditions by specifying the position and velocity of two particles, then -combine them into a single `~gala.dynamics.PhaseSpacePosition` object:: +Like with `~gala.potential.hamiltonian.Hamiltonian` orbit integration, orbital +initial conditions are passed in to `~gala.dynamics.nbody.DirectNBody` by +passing in a single `~gala.dynamics.PhaseSpacePosition` object. Let's create two +initial conditions by specifying the position and velocity of two particles, +then combine them into a single `~gala.dynamics.PhaseSpacePosition` object:: >>> w0_1 = gd.PhaseSpacePosition(pos=[0, 0, 0] * u.pc, ... vel=[0, 1.5, 0] * u.km/u.s) @@ -61,10 +60,10 @@ combine them into a single `~gala.dynamics.PhaseSpacePosition` object:: (2,) We'll then treat particle 1 as a massive object by sourcing a -`~gala.potential.HernquistPotential` at the location of the particle, and -particle 2 as a test particle: To treat some particles as test particles, you -can pass ``None`` or a `~gala.potential.NullPotential` instance for the -corresponding particle potential:: +`~gala.potential.potential.HernquistPotential` at the location of the particle, +and particle 2 as a test particle: To treat some particles as test particles, +you can pass ``None`` or a `~gala.potential.potential.NullPotential` instance +for the corresponding particle potential:: >>> pot1 = gp.HernquistPotential(m=1e7*u.Msun, c=0.5*u.kpc, units=galactic) >>> particle_pot = [pot1, None] @@ -85,6 +84,7 @@ orbits:: .. plot:: :align: center :context: close-figs + :width: 60% import astropy.units as u import numpy as np @@ -113,13 +113,13 @@ orbits:: Example: N-body integration with a background potential ======================================================= -With `~gala.dynamics.DirectNBody`, we can also specify a background or external -potential to integrate all orbits in. To do this, you can optionally pass in an -external potential as a potential object to `~gala.dynamics.DirectNBody`. Here, -as an example, we'll repeat a similar integration as above, but (1) add a -positional offset of the initial conditions from the origin, and (2) specify an -external potential using the `~gala.potential.MilkyWayPotential` class as an -external potential:: +With `~gala.dynamics.nbody.DirectNBody`, we can also specify a background or +external potential to integrate all orbits in. To do this, you can optionally +pass in an external potential as a potential object to +`~gala.dynamics.nbody.DirectNBody`. Here, as an example, we'll repeat a similar +integration as above, but (1) add a positional offset of the initial conditions +from the origin, and (2) specify an external potential using the +`~gala.potential.potential.MilkyWayPotential` class as an external potential:: >>> external_pot = gp.MilkyWayPotential() >>> w0_1 = gd.PhaseSpacePosition(pos=[10, 0, 0] * u.kpc, @@ -138,6 +138,7 @@ external potential:: .. plot:: :align: center :context: close-figs + :width: 60% external_pot = gp.MilkyWayPotential() w0_1 = gd.PhaseSpacePosition(pos=[10, 0, 0] * u.kpc, @@ -167,6 +168,7 @@ should give us a sense of whether particle 2 is bound or unbound to this mass:: .. plot:: :align: center :context: close-figs + :width: 60% dxyz = orbits[:, 0].xyz - orbits[:, 1].xyz diff --git a/docs/dynamics/nd-representations.rst b/docs/dynamics/nd-representations.rst index 6c4998f8..2b8480e7 100644 --- a/docs/dynamics/nd-representations.rst +++ b/docs/dynamics/nd-representations.rst @@ -6,8 +6,7 @@ N-dimensional representation classes ************************************ -For code blocks below and any pages linked below, we'll assume the following -imports have already been executed:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np @@ -17,10 +16,10 @@ Introduction ============ The Astropy |astropyrep|_ presently only support 3D positions and differential -objects. The `~gala.dynamics.NDCartesianRepresentation` and -`~gala.dynamics.NDCartesianDifferential` classes add Cartesian representation -classes that can handle arbitrary numbers of dimensions. For example, 2D -coordinates:: +objects. The `~gala.dynamics.representation_nd.NDCartesianRepresentation` and +`~gala.dynamics.representation_nd.NDCartesianDifferential` classes add Cartesian +representation classes that can handle arbitrary numbers of dimensions. For +example, 2D coordinates:: >>> xy = np.arange(16).reshape(2, 8) * u.kpc >>> rep = gd.NDCartesianRepresentation(xy) @@ -48,6 +47,7 @@ core representation objects:: .. plot:: :align: center + :width: 60% import astropy.units as u import numpy as np diff --git a/docs/dynamics/nonlinear.rst b/docs/dynamics/nonlinear.rst index 503a6dbd..2b2bf087 100644 --- a/docs/dynamics/nonlinear.rst +++ b/docs/dynamics/nonlinear.rst @@ -1,7 +1,7 @@ .. _gala-nonlinear-dynamics: ****************** -Nonlinear dynamics +Nonlinear Dynamics ****************** Introduction @@ -27,11 +27,12 @@ Chaotic orbit ------------- There are two ways to compute Lyapunov exponents implemented in `gala.dynamics`. -In most cases, you'll want to use the `~gala.dynamics.fast_lyapunov_max` function -because the integration is implemented in C and is quite fast. This function only -works if the potential you are working with is implemented in C (e.g., it is a -`~gala.potential.CPotentialBase` subclass). With a potential object and -a set of initial conditions:: +In most cases, you'll want to use the +`~gala.dynamics.nonlinear.fast_lyapunov_max` function because the integration is +implemented in C and is quite fast. This function only works if the potential +you are working with is implemented in C (e.g., it is a +`~gala.potential.potential.CPotentialBase` subclass). With a potential object +and a set of initial conditions:: >>> pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc, ... q1=1., q2=0.8, q3=0.6, units=galactic) @@ -45,13 +46,13 @@ contains the maximum Lyapunov exponent estimate for each offset orbit, argument) and an `~gala.dynamics.Orbit` object that contains the parent orbit and each offset orbit. Let's plot the parent orbit:: - >>> fig = orbit[:,0].plot(marker=',', alpha=0.1, linestyle='none') # doctest: +SKIP + >>> fig = orbit[:,0].plot(marker=',', alpha=0.25, linestyle='none') # doctest: +SKIP .. plot:: :align: center import astropy.units as u - import matplotlib.pyplot as pl + import matplotlib.pyplot as plt import gala.potential as gp import gala.dynamics as gd from gala.units import galactic @@ -61,7 +62,7 @@ the parent orbit and each offset orbit. Let's plot the parent orbit:: w0 = gd.PhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc, vel=[0.,100.,0]*u.km/u.s) lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., n_steps=100000) - fig = orbit[:,0].plot(marker=',', linestyle='none', alpha=0.1) + fig = orbit[:,0].plot(marker=',', linestyle='none', alpha=0.25) Visually, this looks like a chaotic orbit. This means the Lyapunov exponent should saturate to some value. We'll now plot the estimate of the Lyapunov @@ -70,17 +71,18 @@ several time-steps (controllable with the ``n_steps_per_pullback`` argument), we have to down-sample the time array to align it with the Lyapunov exponent array. This plots one line per offset orbit:: - >>> pl.figure() # doctest: +SKIP - >>> pl.loglog(orbit.t[11::10], lyap, marker='') # doctest: +SKIP - >>> pl.xlabel("Time [{}]".format(orbit.t.unit)) # doctest: +SKIP - >>> pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) # doctest: +SKIP - >>> pl.tight_layout() # doctest: +SKIP + >>> plt.figure() # doctest: +SKIP + >>> plt.loglog(orbit.t[11::10], lyap, marker='') # doctest: +SKIP + >>> plt.xlabel("Time [{}]".format(orbit.t.unit)) # doctest: +SKIP + >>> plt.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) # doctest: +SKIP + >>> plt.tight_layout() # doctest: +SKIP .. plot:: :align: center + :width: 60% import astropy.units as u - import matplotlib.pyplot as pl + import matplotlib.pyplot as plt import gala.potential as gp import gala.dynamics as gd from gala.units import galactic @@ -91,11 +93,11 @@ array. This plots one line per offset orbit:: vel=[0.,100.,0]*u.km/u.s) lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., n_steps=100000) - pl.figure() - pl.loglog(orbit.t[11::10], lyap, marker='') - pl.xlabel("Time [{}]".format(orbit.t.unit)) - pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) - pl.tight_layout() + plt.figure() + plt.loglog(orbit.t[11::10], lyap, marker='') + plt.xlabel("Time [{}]".format(orbit.t.unit)) + plt.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) + plt.tight_layout() The estimate is clearly starting to diverge from a simple power law decay. @@ -113,7 +115,6 @@ To compare, we will compute the estimate for a regular orbit as well:: :align: center import astropy.units as u - import matplotlib.pyplot as pl import gala.potential as gp import gala.dynamics as gd from gala.units import galactic @@ -136,6 +137,7 @@ following a characteristic power-law (a straight line in a log-log plot):: .. plot:: :align: center + :width: 60% import astropy.units as u import matplotlib.pyplot as pl diff --git a/docs/dynamics/orbits-in-detail.rst b/docs/dynamics/orbits-in-detail.rst index 89ee63f3..38b98c64 100644 --- a/docs/dynamics/orbits-in-detail.rst +++ b/docs/dynamics/orbits-in-detail.rst @@ -6,7 +6,7 @@ Orbit and phase-space position objects in more detail ***************************************************** -We'll assume the following imports have already been executed:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np @@ -40,7 +40,7 @@ reading below: .. _phase-space-position: -Phase-space positions +Phase-space Positions ===================== The |psp| class provides an interface for representing full phase-space @@ -52,16 +52,16 @@ The easiest way to create a |psp| object is to pass in a pair of `~astropy.units.Quantity` objects that represent the Cartesian position and velocity vectors:: - >>> gd.PhaseSpacePosition(pos=[4.,8.,15.]*u.kpc, - ... vel=[-150.,50.,15.]*u.km/u.s) + >>> gd.PhaseSpacePosition(pos=[4., 8., 15.] * u.kpc, + ... vel=[-150., 50., 15.] * u.km/u.s) By default, passing in `~astropy.units.Quantity`'s are interpreted as Cartesian coordinates and velocities. This works with arrays of positions and velocities as well:: - >>> x = np.arange(24).reshape(3,8) - >>> v = np.arange(24).reshape(3,8) + >>> x = np.arange(24).reshape(3, 8) + >>> v = np.arange(24).reshape(3, 8) >>> w = gd.PhaseSpacePosition(pos=x * u.kpc, ... vel=v * u.km/u.s) >>> w @@ -90,8 +90,10 @@ These mappings are defined in the class definition of `~gala.dynamics.PhaseSpacePosition`. For example, to access the ``x`` component of the position and the ``v_x`` component of the velocity:: - >>> w.x, w.v_x # doctest: +FLOAT_CMP - (, ) + >>> w.x # doctest: +FLOAT_CMP + + >>> w.v_x # doctest: +FLOAT_CMP + The default representation is Cartesian, but the class can also be instantiated with representation objects instead of `~astropy.units.Quantity`'s -- this is @@ -174,6 +176,7 @@ function and any keyword arguments are passed through to that function:: .. plot:: :align: center :context: close-figs + :width: 60% fig = w.plot(components=['x', 'v_z'], color='r', facecolor='none', marker='o', s=20, alpha=0.5) @@ -212,7 +215,7 @@ orbits, each with the same 128 times:: >>> t = np.linspace(0, 100, 128) * u.Myr >>> Om = np.random.uniform(size=16) * u.rad / u.Myr - >>> angle = Om[None] * t[:,None] + >>> angle = Om[None] * t[:, None] >>> pos = np.stack((5*np.cos(angle), np.sin(angle))).value * u.kpc >>> vel = np.stack((-5*np.sin(angle), np.cos(angle))).value * u.kpc/u.Myr >>> orbit = gd.Orbit(pos=pos, vel=vel) @@ -220,8 +223,8 @@ orbits, each with the same 128 times:: To make full use of the orbit functionality, you must also pass in an array with -the time values and an instance of a `~gala.potential.PotentialBase` subclass -that represents the potential that the orbit was integrated in:: +the time values and an instance of a `~gala.potential.potential.PotentialBase` +subclass that represents the potential that the orbit was integrated in:: >>> pot = gp.PlummerPotential(m=1E10, b=1., units=galactic) >>> orbit = gd.Orbit(pos=pos*u.kpc, vel=vel*u.km/u.s, @@ -230,10 +233,10 @@ that represents the potential that the orbit was integrated in:: (note, in this case ``pos`` and ``vel`` were not generated from integrating an orbit in the potential ``pot``!). However, most of the time you won't need to create |orb| objects from scratch! They are returned from any of the numerical -intergration routines provided in `gala`. For example, they are returned by the -`~gala.potential.PotentialBase.integrate_orbit` method of potential objects and -will automatically contain the ``time`` array and ``potential`` object. For -example:: +integration routines provided in `gala`. For example, they are returned by the +`~gala.potential.potential.PotentialBase.integrate_orbit` method of potential +objects and will automatically contain the ``time`` array and ``potential`` +object. For example:: >>> pot = gp.PlummerPotential(m=1E10 * u.Msun, b=1. * u.kpc, units=galactic) >>> w0 = gd.PhaseSpacePosition(pos=[10.,0,0] * u.kpc, diff --git a/docs/integrate/index.rst b/docs/integrate/index.rst index ef09f03b..89ac8363 100644 --- a/docs/integrate/index.rst +++ b/docs/integrate/index.rst @@ -9,34 +9,35 @@ Integration (`gala.integrate`) Introduction ============ -Scipy provides numerical ODE integration functions but they aren't very -user friendly or object oriented. This subpackage implements the Leapfrog -integration scheme (not available in Scipy) and provides wrappers to -higher order integration schemes such as a 5th order Runge-Kutta and the -Dormand-Prince 85(3) method. +:mod:`scipy` provides numerical ODE integration functions (e.g., +:func:`scipy.integrate.odeint` and :func:`scipy.integrate.solve_ivp`), but these +functions are not object-oriented or accessible from C. The +:mod:`gala.integrate` subpackage implements the Leapfrog integration scheme (not +available in Scipy) and provides C wrappers for higher order integration schemes +such as a 5th order Runge-Kutta and the Dormand-Prince 85(3) method. -For code blocks below and any pages linked below, I assume the following -imports have already been executed:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np + >>> import gala.dynamics as gd >>> import gala.integrate as gi - >>> from gala.units import galactic + >>> from gala.units import galactic, UnitSystem Getting Started =============== -All of the integrator classes have the same basic call structure. To create an -integrator object, you pass in a function that evaluates derivatives of, for -example, phase-space coordinates, then you call the +All of the integrator classes in :mod:`gala.integrate` have the same basic call +structure. To create an integrator object, you pass in a function that evaluates +derivatives of, for example, phase-space coordinates, then you call the `~gala.integrate.Integrator.run` method while specifying timestep information. The integration function must accept, at minimum, two arguments: the current time, ``t``, and the current position in phase-space, ``w``. The time is a -single floating-point number and the position will have shape ``(ndim, -norbits)`` where ``ndim`` is the full dimensionality of the phase-space (e.g., 6 -for a 3D coordinate system) and ``norbits`` is the number of orbits. These -inputs will *not* have units associated with them (e.g., they are not +single floating-point number and the phase-space position will have shape +``(ndim, norbits)`` where ``ndim`` is the full dimensionality of the phase-space +(e.g., 6 for a 3D coordinate system) and ``norbits`` is the number of orbits. +These inputs will *not* have units associated with them (e.g., they are not :class:`astropy.units.Quantity` objects). An example of such a function (that represents a simple harmonic oscillator) is:: @@ -45,20 +46,20 @@ represents a simple harmonic oscillator) is:: ... return np.array([x_dot, -x]) Even though time does not explicitly enter into the equation, the function must -still accept a time argument. We will now create an instance of -`~gala.integrate.LeapfrogIntegrator` to integrate an orbit:: +still accept a time argument. We can now create an instance of +`~gala.integrate.LeapfrogIntegrator` to integrate an orbit in a harmonic +oscillator potential:: >>> integrator = gi.LeapfrogIntegrator(F) -To actually run the integrator, we need to specify a set of initial conditions. -The simplest way to do this is to specify an array:: +To run the integrator, we need to specify a set of initial conditions. The +simplest way to do this is to specify an array:: - >>> w0 = np.array([1.,0.]) + >>> w0 = np.array([1., 0.]) -However this then causes the integrator to work without units. The orbit object -returned by the integrator will then also have no associated units. For example, -to integrate from these initial conditions with a time step of 0.5 for 100 -steps:: +This causes the integrator to work without units, so the orbit object returned +by the integrator will then also have no associated units. For example, to +integrate from these initial conditions with a time step of 0.5 for 100 steps:: >>> orbit = integrator.run(w0, dt=0.5, n_steps=100) >>> orbit.t.unit @@ -66,59 +67,63 @@ steps:: >>> orbit.pos.xyz.unit Unit(dimensionless) -We can instead specify the unit system that the function (``F``) expects, and +We could instead specify the unit system that the function (``F``) expects, and then pass in a `~gala.dynamics.PhaseSpacePosition` object with arbitrary units -in as initial conditions:: +as initial conditions:: - >>> import gala.dynamics as gd - >>> from gala.units import UnitSystem >>> usys = UnitSystem(u.m, u.s, u.kg, u.radian) >>> integrator = gi.LeapfrogIntegrator(F, func_units=usys) >>> w0 = gd.PhaseSpacePosition(pos=[100.]*u.cm, vel=[0]*u.cm/u.yr) >>> orbit = integrator.run(w0, dt=0.5, n_steps=100) + +The returned orbit object has quantities in the specified unit system, for +example:: + >>> orbit.t.unit Unit("s") + >>> orbit.x1.unit + Unit("m") -Now the orbit object will have quantities in the specified unit system. Example: Forced pendulum ------------------------- -Here we demonstrate how to use the Dormand-Prince integrator to compute the +Here we will demonstrate how to use the Dormand-Prince integrator to compute the orbit of a forced pendulum. We will use the variable ``q`` as the angle of the pendulum with respect to the vertical and ``p`` as the conjugate momentum. Our Hamiltonian is .. math:: - H(q,p) = \frac{1}{2}p^2 + \cos(q) + A\sin(\omega_D t) + H(q, p) = \frac{1}{2} \, p^2 + \cos(q) + A \, \sin(\omega_D \, t) so that .. math:: \dot{q} &= p\\ - \dot{p} &= -\sin(q) + A\omega_D\cos(\omega_D t) + \dot{p} &= -\sin(q) + A\, \omega_D \, \cos(\omega_D \, t) -Our derivative function is then:: +For numerical integration, the function to compute the time derivatives of our +phase-space coordinates is then:: >>> def F(t, w, A, omega_D): ... q, p = w ... wdot = np.zeros_like(w) ... wdot[0] = p - ... wdot[1] = -np.sin(q) + A*omega_D*np.cos(omega_D*t) + ... wdot[1] = -np.sin(q) + A * omega_D * np.cos(omega_D * t) ... return wdot -This function has two arguments -- :math:`A` (``A``), the amplitude of the -forcing, and :math:`\omega_D` (``omega_D``), the driving frequency. We define an -integrator object by specifying this function, along with values for the -function arguments:: +This function has two arguments: :math:`A` (``A``), the amplitude of the +forcing,and :math:`\omega_D` (``omega_D``), the driving frequency. We define an +integrator object by specifying this function along with values for the function +arguments:: >>> integrator = gi.DOPRI853Integrator(F, func_args=(0.07, 0.75)) To integrate an orbit, we use the `~gala.integrate.Integrator.run` method. We have to specify the initial conditions along with information about how long to -integrate and with what stepsize. There are several options for how to specify +integrate and with what step size. There are several options for how to specify the time step information. We could pre-generate an array of times and pass that in, or pass in an initial time, end time, and timestep. Or, we could simply pass in the number of steps to run for and a timestep. For this example, we will use @@ -129,7 +134,7 @@ information.:: We can plot the integrated (chaotic) orbit:: - >>> fig = orbit.plot(subplots_kwargs=dict(figsize=(8,4))) # doctest: +SKIP + >>> fig = orbit.plot(subplots_kwargs=dict(figsize=(8, 4))) # doctest: +SKIP .. plot:: :align: center @@ -154,7 +159,7 @@ We can plot the integrated (chaotic) orbit:: Example: Lorenz equations ------------------------- -Here's another example of integrating the +Here's another example of numerical ODE integration using the `Lorenz equations `_, a 3D nonlinear system:: @@ -185,7 +190,7 @@ nonlinear system:: sigma, rho, beta = 10., 28., 8/3. integrator = gi.DOPRI853Integrator(F, func_args=(sigma, rho, beta)) - orbit = integrator.run([0.5,0.5,0.5,0,0,0], dt=1E-2, n_steps=1E4) + orbit = integrator.run([0.5, 0.5, 0.5, 0, 0, 0], dt=1E-2, n_steps=1E4) fig = orbit.plot() API diff --git a/docs/interop.rst b/docs/interop.rst index c319a714..f0a3c71a 100644 --- a/docs/interop.rst +++ b/docs/interop.rst @@ -63,8 +63,8 @@ Similarly, a Galpy ``Orbit`` can be converted to a Gala Gala also provides tools for converting potential objects to `galpy` potential objects, or creating Gala potential objects from existing `galpy` potentials. To convert a Gala potential to a Galpy potential, use the -:meth:`~gala.potential.PotentialBase.to_galpy_potential()` method on any Gala -potential object. For example: +:meth:`~gala.potential.potential.PotentialBase.to_galpy_potential()` method on +any Gala potential object. For example: .. doctest-requires:: galpy diff --git a/docs/potential/compositepotential.rst b/docs/potential/compositepotential.rst index bad3056e..0dc20129 100644 --- a/docs/potential/compositepotential.rst +++ b/docs/potential/compositepotential.rst @@ -16,9 +16,10 @@ implemented in C, it is always faster to use :class:`~gala.potential.potential.CCompositePotential`, where the composition is done at the C layer rather than in Python. -With either class, interaction with the class is identical. To compose -potentials with unique but arbitrary names, you can simply add pre-defined -potential class instances:: +With either class, interaction with the class (e.g., by calling methods) is +identical to the individual potential classes. To compose potentials with unique +but arbitrary names, you can also simply add pre-defined potential class +instances:: >>> import numpy as np >>> import gala.potential as gp @@ -56,27 +57,25 @@ is equivalent to:: >>> pot['disk'] = disk >>> pot['bulge'] = bulge -In detail, the composite potential classes subclass -:class:`~collections.OrderedDict`, so in this sense there is a slight difference -between the two examples above (if you are using a Python version < 3.6). By -defining components after creating the instance, the order is preserved. In the -above example, the disk potential would always be called first and the bulge -would always be called second. +The order of insertion is preserved, and sets the order that the potentials are +called. In the above example, the disk potential would always be called first +and the bulge would always be called second. The resulting potential object has all of the same properties as individual potential objects:: - >>> pot.energy([1.,-1.,0.]) # doctest: +FLOAT_CMP + >>> pot.energy([1., -1., 0.]) # doctest: +FLOAT_CMP - >>> pot.acceleration([1.,-1.,0.]) # doctest: +FLOAT_CMP + >>> pot.acceleration([1., -1., 0.]) # doctest: +FLOAT_CMP - >>> grid = np.linspace(-3.,3.,100) - >>> fig = pot.plot_contours(grid=(grid,0,grid)) # doctest: +SKIP + >>> grid = np.linspace(-3., 3., 100) + >>> fig = pot.plot_contours(grid=(grid, 0, grid)) # doctest: +SKIP .. plot:: :align: center + :width: 60% import numpy as np import gala.dynamics as gd diff --git a/docs/potential/define-new-potential.rst b/docs/potential/define-new-potential.rst index 43874953..4514f338 100644 --- a/docs/potential/define-new-potential.rst +++ b/docs/potential/define-new-potential.rst @@ -12,13 +12,12 @@ and Cython. The advantage to writing a new class in Cython is that the computations can execute with C-like speeds, however only certain integrators support using this functionality (Leapfrog and DOP853) and it is a bit more complicated to set up the code to build the C+Cython code properly. If you are -not familiar with Cython, you probably want to stick to a pure-Python class for +not familiar with Cython, you probably want to stick to a pure Python class for initial testing. If there is a potential class that you think should be -included, feel free to suggest the new addition as a `GitHub issue -`_! +included as a built-in Cython potential, feel free to suggest the new addition +as a `GitHub issue `_! -For code blocks below and any pages linked below, we assume the following -imports have already been executed:: +For the examples below the following imports have already been executed:: >>> import numpy as np >>> import gala.potential as gp @@ -43,24 +42,25 @@ The expression for the potential is: With this parametrization, there is only one free parameter (``A``), and the potential is two-dimensional. -At minimum, the subclass must implement: +At minimum, the subclass must implement the following methods: - ``__init__()`` - ``_energy()`` - ``_gradient()`` -The ``_energy()`` method will compute the potential energy at a given position -and time. The ``_gradient()`` computes the gradient of the potential. Both of -these methods must accept two arguments: a position, and a time. These internal -methods are then called by the :class:`~gala.potential.potential.PotentialBase` -superclass methods :meth:`~gala.potential.potential.PotentialBase.energy` and +The ``_energy()`` method should compute the potential energy at a given position +and time. The ``_gradient()`` method should compute the gradient of the +potential. Both of these methods must accept two arguments: a position, and a +time. These internal methods are then called by the +:class:`~gala.potential.potential.PotentialBase` superclass methods +:meth:`~gala.potential.potential.PotentialBase.energy` and :meth:`~gala.potential.potential.PotentialBase.gradient`. The superclass methods -convert the input position to an array in the unit system of the potetial for -fast evalutaion. The input to these superclass methods can be +convert the input position to an array in the unit system of the potential for +fast evaluation. The input to these superclass methods can be :class:`~astropy.units.Quantity` objects, :class:`~gala.dynamics.PhaseSpacePosition` objects, or :class:`~numpy.ndarray`. -Because this potential has a free parameter, the ``__init__`` method must accept +Because this potential has a parameter, the ``__init__`` method must accept a parameter argument and store this in the ``parameters`` dictionary attribute (a required attribute of any subclass). Let's write it out, then work through what each piece means in detail:: @@ -85,12 +85,12 @@ what each piece means in detail:: The internal energy and gradient methods compute the numerical value and gradient of the potential. The ``__init__`` method must take a single argument, -``A``, and store this to a paremeter dictionary. The expected shape of the +``A``, and store this to a parameter dictionary. The expected shape of the position array (``xy``) passed to the internal ``_energy()`` and ``_gradient()`` methods is always 2-dimensional with shape ``(n_points, n_dim)`` where ``n_points >= 1`` and ``n_dim`` must match the dimensionality of the potential specified in the initializer. Note that this is different from the shape -expected when calling the actual methods ``energy()`` and ``gradient()``! +expected when calling the public methods ``energy()`` and ``gradient()``! Let's now create an instance of the class and see how it works. For now, let's pass in ``None`` for the unit system to designate that we'll work with @@ -98,17 +98,20 @@ dimensionless quantities:: >>> pot = CustomHenonHeilesPotential(A=1., units=None) -That's it! Now we have a fully-fledged potential object. For example, we -can integrate an orbit in this potential:: +That's it! We now have a potential object with all of the same functionality as +the built-in potential classes. For example, we can integrate an orbit in this +potential (but note that this potential is two-dimensional, so we only have to +specify four coordinate values):: - >>> w0 = gd.PhaseSpacePosition(pos=[0.,0.3], - ... vel=[0.38,0.]) + >>> w0 = gd.PhaseSpacePosition(pos=[0., 0.3], + ... vel=[0.38, 0.]) >>> orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.05, n_steps=10000) >>> fig = orbit.plot(marker=',', linestyle='none', alpha=0.5) # doctest: +SKIP .. plot:: :align: center :context: close-figs + :width: 60% import matplotlib.pyplot as pl import numpy as np @@ -136,13 +139,13 @@ can integrate an orbit in this potential:: orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.05, n_steps=10000) fig = orbit.plot(marker=',', linestyle='none', alpha=0.5) -Or, we could create a contour plot of equipotentials:: +We could also, for example, create a contour plot of equipotentials:: >>> grid = np.linspace(-1., 1., 100) >>> from matplotlib import colors >>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(1, 1, figsize=(5,5)) - >>> fig = pot.plot_contours(grid=(grid,grid), + >>> fig = pot.plot_contours(grid=(grid, grid), ... levels=np.logspace(-3, 1, 10), ... norm=colors.LogNorm(), ... cmap='Blues', ax=ax) @@ -150,6 +153,7 @@ Or, we could create a contour plot of equipotentials:: .. plot:: :align: center :context: close-figs + :width: 60% from matplotlib import colors import matplotlib.pyplot as plt diff --git a/docs/potential/hamiltonian-reference-frames.rst b/docs/potential/hamiltonian-reference-frames.rst index 03e17b0f..80f118a6 100644 --- a/docs/potential/hamiltonian-reference-frames.rst +++ b/docs/potential/hamiltonian-reference-frames.rst @@ -4,8 +4,7 @@ Hamiltonian objects and reference frames **************************************** -For code blocks below, I assume the following imports have already been -excuted:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np @@ -28,7 +27,7 @@ When :ref:`integrating orbits using the potential classes directly >>> orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=1000) it is implicitly assumed that the initial conditions and orbit are in an -intertial (static) reference frame. In this case, the total energy or value +inertial (static) reference frame. In this case, the total energy or value of the Hamiltonian (per unit mass) is simply .. math:: @@ -62,7 +61,7 @@ Using the potential objects and :meth:`~gala.potential.potential.PotentialBase.integrate_orbit()` to integrate an orbit is equivalent to defining a :class:`~gala.potential.hamiltonian.Hamiltonian` object with the potential -object and a :class:`~gala.potential.frame.Staticframe` instance:: +object and a :class:`~gala.potential.frame.builtin.Staticframe` instance:: >>> pot = gp.HernquistPotential(m=1E10*u.Msun, c=1.*u.kpc, ... units=galactic) diff --git a/docs/potential/index.rst b/docs/potential/index.rst index 4fe16423..b5f943f5 100644 --- a/docs/potential/index.rst +++ b/docs/potential/index.rst @@ -11,22 +11,22 @@ Introduction This subpackage provides a number of classes for working with parametric models of gravitational potentials. There are a number of built-in potentials -implemented in C and Cython for speed and there are base classes that allow for -easy creation of `new custom potential classes `_ in -pure-Python. The potential objects have methods for computing, for example, the -potential energy, gradient, density, or mass profiles. These are particularly -useful in combination with the `~gala.integrate` and `~gala.dynamics` -subpackages. +implemented in C and Cython (for speed), and there are base classes that allow +for easy creation of `new custom potential classes `_ +in pure Python or by writing custom C/Cython extensions. The ``Potential`` +objects have convenience methods for computing common dynamical quantities, for +example: potential energy, spatial gradient, density, or mass profiles. These +are particularly useful in combination with the `~gala.integrate` and +`~gala.dynamics` subpackages. Also defined in this subpackage are a set of reference frames which can be used for numerical integration of orbits in non-static reference frames. See the page -on :ref:`hamiltonian-reference-frames` for more information. Potential objects -can be combined with a reference frame and stored in a -`~gala.potential.Hamiltonian` object that provides an easy interface to -numerical orbit integration. +on :ref:`hamiltonian-reference-frames` for more information. ``Potential`` +objects can be combined with a reference frame and stored in a +`~gala.potential.hamiltonian.Hamiltonian` object that provides an easy interface +to numerical orbit integration. -For code blocks below and any pages linked below, I assume the following imports -have already been excuted:: +For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import matplotlib.pyplot as plt @@ -34,22 +34,20 @@ have already been excuted:: >>> import gala.potential as gp >>> from gala.units import galactic, solarsystem, dimensionless -Getting started: built-in methods of potential classes +Getting Started: Built-in Methods of Potential Classes ====================================================== -The built-in potentials are all initialized by passing in keyword argument -parameter values as :class:`~astropy.units.Quantity` objects or as numeric -values in a specified unit system. To see what parameters are available for a -given potential, check the documentation for the individual classes below. You -must also specify a `~gala.units.UnitSystem` when initializing a potential. A -unit system is a set of non-reducible units that define (at minimum) the length, -mass, time, and angle units. A few common unit systems are built in to the -package (e.g., ``galactic``, ``solarsystem``, ``dimensionless``). - -All of the built-in potential objects have defined methods to evaluate the -potential energy and the gradient/acceleration at a given position(s). -For example, here we will create a potential object for a 2D point mass located -at the origin with unit mass:: +Any of the built-in ``Potential`` classes are initialized by passing in keyword +argument parameter values as :class:`~astropy.units.Quantity` objects or as +numeric values in a specified unit system. To see what parameters are available +for a given potential, check the documentation for the individual classes below. +You must also specify a `~gala.units.UnitSystem` when initializing a potential. +A unit system is a set of non-reducible units that define (at minimum) the +length, mass, time, and angle units. A few common unit systems are built in to +the package (e.g., ``galactic``, ``solarsystem``, ``dimensionless``). For +example, to create an object to represent a Kepler potential (point mass) at the +origin with mass = 1 solar mass, we would instantiate a +:class:`~gala.potential.potential.KeplerPotential` object: >>> ptmass = gp.KeplerPotential(m=1.*u.Msun, units=solarsystem) >>> ptmass @@ -61,60 +59,64 @@ specified unit system:: >>> gp.KeplerPotential(m=1047.6115*u.Mjup, units=solarsystem) -If no units are specified for a parameter, it is assumed to already be -consistent with the `~gala.units.UnitSystem` passed in:: +If no units are specified for a parameter (i.e. a parameter value is passed in +as a Python numeric value or array), it is assumed to be in the specified +`~gala.units.UnitSystem`:: >>> gp.KeplerPotential(m=1., units=solarsystem) The potential classes work well with the :mod:`astropy.units` framework, but to ignore units you can use the `~gala.units.DimensionlessUnitSystem` or pass -``None`` as the unit system:: +`None` as the unit system:: >>> gp.KeplerPotential(m=1., units=None) -With an instantiated potential object, we can evaluate the potential energy at -some position or a set of positions. Here, we'll evaluate a 3D point mass -potential at the 3D position (1,-1,0):: +All of the built-in potential objects have defined methods to evaluate the +potential energy and the gradient/acceleration at a given position or array of +positions. For example, to evaluate the potential energy at the 3D position +``(x, y, z) = (1, -1, 0) AU``:: - >>> ptmass.energy([1.,-1.,0.]*u.au) + >>> ptmass.energy([1., -1., 0.] * u.au) These functions also accept both :class:`~astropy.units.Quantity` objects or plain :class:`~numpy.ndarray`-like objects (in which case the position is assumed to be in the unit system of the potential):: - >>> ptmass.energy([1.,-1.,0.]) + >>> ptmass.energy([1., -1., 0.]) This also works for multiple positions by passing in a 2D position (but see :ref:`conventions` for a description of the interpretation of different axes):: - >>> pos = np.array([[1.,-1.,0], - ... [2.,3.,0]]).T - >>> ptmass.energy(pos*u.au) + >>> pos = np.array([[1., -1. ,0], + ... [2., 3., 0]]).T + >>> ptmass.energy(pos * u.au) -We may also compute the gradient or acceleration:: +We can also compute the gradient or acceleration:: - >>> ptmass.gradient([1.,-1.,0]*u.au) # doctest: +FLOAT_CMP + >>> ptmass.gradient([1., -1., 0] * u.au) # doctest: +FLOAT_CMP - >>> ptmass.acceleration([1.,-1.,0]*u.au) # doctest: +FLOAT_CMP + >>> ptmass.acceleration([1., -1., 0] * u.au) # doctest: +FLOAT_CMP -Some of the potential objects also have methods implemented for computing the -corresponding mass density and the Hessian of the potential (matrix of 2nd -derivatives) at given locations. For example:: +Most of the potential objects also have methods implemented for computing the +corresponding mass density and the Hessian of the potential (the matrix of 2nd +derivatives) at given locations. For example, with the +:class:`~gala.potential.potential.HernquistPotential`, we can evaluate both the +mass density and Hessian at the position ``(x, y, z) = (1, -1, 0) kpc``:: >>> pot = gp.HernquistPotential(m=1E9*u.Msun, c=1.*u.kpc, units=galactic) - >>> pot.density([1.,-1.,0]*u.kpc) # doctest: +FLOAT_CMP + >>> pot.density([1., -1., 0] * u.kpc) # doctest: +FLOAT_CMP - >>> pot.hessian([1.,-1.,0]*u.kpc) # doctest: +SKIP + >>> pot.hessian([1., -1., 0] * u.kpc) # doctest: +SKIP -Another useful method is :meth:`~gala.potential.PotentialBase.mass_enclosed`, -which numerically estimates the mass enclosed within a spherical shell defined -by the specified position. This numerically estimates -:math:`\frac{d \Phi}{d r}` along the vector pointing at the specified position -and estimates the enclosed mass simply as -:math:`M(>> pot = gp.NFWPotential(m=1E11*u.Msun, r_s=20.*u.kpc, units=galactic) >>> pos = np.zeros((3,100)) * u.kpc @@ -146,6 +148,7 @@ be used to compute, for example, a mass profile:: .. plot:: :align: center :context: close-figs + :width: 60% import astropy.units as u import numpy as np @@ -164,18 +167,20 @@ be used to compute, for example, a mass profile:: plt.ylabel("$M(>> p = gp.MiyamotoNagaiPotential(m=1E11, a=6.5, b=0.27, units=galactic) - >>> fig, ax = plt.subplots(1,1) # doctest: +SKIP + >>> fig, ax = plt.subplots() # doctest: +SKIP >>> p.plot_contours(grid=(np.linspace(-15,15,100), 0., 1.), marker='', ax=ax) # doctest: +SKIP >>> E_unit = p.units['energy'] / p.units['mass'] >>> ax.set_xlabel("$x$ [{}]".format(p.units['length'].to_string(format='latex'))) # doctest: +SKIP @@ -184,6 +189,7 @@ the potential value as a function of :math:`x` position at :math:`y=0, z=1`:: .. plot:: :align: center :context: close-figs + :width: 90% pot = gp.MiyamotoNagaiPotential(m=1E11, a=6.5, b=0.27, units=galactic) @@ -200,9 +206,9 @@ for :math:`z` (the meshgridding is taken care of internally). Here, we choose to draw on a pre-defined `matplotlib` axes object so we can set the labels and aspect ratio of the plot:: - >>> fig,ax = plt.subplots(1, 1, figsize=(12,4)) - >>> x = np.linspace(-15,15,100) - >>> z = np.linspace(-5,5,100) + >>> fig,ax = plt.subplots(1, 1, figsize=(12, 4)) + >>> x = np.linspace(-15, 15, 100) + >>> z = np.linspace(-5, 5, 100) >>> p.plot_contours(grid=(x, 1., z), ax=ax) # doctest: +SKIP >>> ax.set_xlabel("$x$ [kpc]") # doctest: +SKIP >>> ax.set_ylabel("$z$ [kpc]") # doctest: +SKIP @@ -210,6 +216,7 @@ aspect ratio of the plot:: .. plot:: :align: center :context: close-figs + :width: 60% x = np.linspace(-15,15,100) z = np.linspace(-5,5,100) @@ -223,12 +230,14 @@ aspect ratio of the plot:: Saving / loading potential objects ================================== -Potential objects can be `pickled `_ +Potential objects can be `pickled `_ and can therefore be stored for later use. However, pickles are saved as binary files. It may be useful to save to or load from text-based specifications of -Potential objects. This can be done with :func:`gala.potential.save` and -:func:`gala.potential.load`, or with the :meth:`~gala.potential.PotentialBase.save` -and method:: +Potential objects. This can be done with the +:meth:`~gala.potential.potential.PotentialBase.save` method and the +:func:`~gala.potential.potential.load` function, which serialize and +de-serialize (respectively) a ``Potential`` object to a `YAML +`_ file:: >>> from gala.potential import load >>> pot = gp.NFWPotential(m=6E11*u.Msun, r_s=20.*u.kpc, @@ -240,10 +249,11 @@ and method:: Exporting potentials as ``sympy`` expressions ============================================= -Most of the potential classes can be exported to a ``sympy`` expression that can +Most of the potential classes can be exported to a `sympy` expression that can be used to manipulate or evaluate the form of the potential. To access this -functionality, the potential classes have a ``.to_sympy()`` classmethod (note: -this requires that ``sympy`` is installed): +functionality, the potential classes have a +`~gala.potential.potential.PotentialBase.to_sympy` classmethod (note: this +requires `sympy` to be installed): .. doctest-requires:: sympy @@ -259,7 +269,7 @@ in the expression as ``sympy`` symbols, here defined as ``vars_``: >>> vars_ {'x': x, 'y': y, 'z': z} -A second dictionary containing the potential parameters as ``sympy`` symbols is +A second dictionary containing the potential parameters as `sympy` symbols is also returned, here defined as ``pars``: .. doctest-requires:: sympy @@ -268,7 +278,7 @@ also returned, here defined as ``pars``: {'v_c': v_c, 'r_h': r_h, 'q1': q1, 'q2': q2, 'q3': q3, 'phi': phi, 'G': G} The expressions and variables returned can be used to perform operations on the -potential expression. For example, to create a ``sympy`` expression for the +potential expression. For example, to create a `sympy` expression for the gradient of the potential: .. doctest-requires:: sympy diff --git a/docs/potential/origin-rotation.rst b/docs/potential/origin-rotation.rst index 6d678308..ff1ae0e5 100644 --- a/docs/potential/origin-rotation.rst +++ b/docs/potential/origin-rotation.rst @@ -1,23 +1,22 @@ .. _rotate-origin-potential: -For code blocks below, we assume the following imports have already been -executed:: - - >>> import astropy.units as u - >>> import numpy as np - >>> import gala.potential as gp - >>> from gala.units import galactic, solarsystem - ************************************************************** -Specifying rotations or origin shifts to ``Potential`` classes +Specifying rotations or origin shifts in ``Potential`` classes ************************************************************** -Most of the gravitational potential classes implemented in ``gala`` support +Most of the gravitational potential classes implemented in `gala` support shifting the origin of the potential relative to the coordinate system, and specifying a rotation of the potential relative to the coordinate system. By default, the origin is assumed to be at (0,0,0) or (0,0), and there is no rotation assumed. +For the examples below the following imports have already been executed:: + + >>> import astropy.units as u + >>> import numpy as np + >>> import gala.potential as gp + >>> from gala.units import galactic, solarsystem + Origin shifts ============= @@ -34,8 +33,8 @@ one another such that one potential is at ``(1, 0, 0)`` AU and the other is at ... units=solarsystem) To see that these are shifted from the coordinate system origin, let's combine -these two objects into a `~gala.potential.CCompositePotential` and visualize the -potential:: +these two objects into a `~gala.potential.potential.CCompositePotential` and +visualize the potential:: >>> pot = gp.CCompositePotential(p1=p1, p2=p2) >>> fig, ax = plt.subplots(1, 1, figsize=(5, 5)) # doctest: +SKIP @@ -47,6 +46,7 @@ potential:: .. plot:: :align: center :context: close-figs + :width: 60% import astropy.units as u import matplotlib.pyplot as plt @@ -72,10 +72,13 @@ Rotations ========= Rotations can be specified either by passing in a -`scipy.spatial.transform.Rotation` instance, or by passing in a 2D numpy array +`scipy.spatial.transform.Rotation` instance, or by passing in a 2D `numpy` array specifying a rotation matrix. For example, let's see what happens if we rotate a -bar potential using these two possible inputs. First, we'll define a rotation matrix specifying a 30 degree -rotation around the z axis (i.e. counter-clockwise) using `astropy.coordinates.matrix_utilities.rotation_matrix`. Next, we'll define a rotation using a scipy `~scipy.spatial.transform.Rotation` object:: +bar potential using these two possible inputs. First, we'll define a rotation +matrix specifying a 30 degree rotation around the z axis (i.e. +counter-clockwise) using `astropy.coordinates.matrix_utilities.rotation_matrix`. +Next, we'll define a rotation using a `scipy` +`~scipy.spatial.transform.Rotation` object:: >>> from astropy.coordinates.matrix_utilities import rotation_matrix >>> from scipy.spatial.transform import Rotation diff --git a/docs/potential/scf-examples.rst b/docs/potential/scf-examples.rst index 9853231f..e421eb3c 100644 --- a/docs/potential/scf-examples.rst +++ b/docs/potential/scf-examples.rst @@ -2,7 +2,7 @@ Examples ******** -For code blocks below, the following imports will be needed:: +For the examples below the following imports have already been executed:: import astropy.units as u import matplotlib as mpl diff --git a/docs/testing.rst b/docs/testing.rst index 0efa871d..a6adeb01 100644 --- a/docs/testing.rst +++ b/docs/testing.rst @@ -4,19 +4,33 @@ Running the tests ================= -The tests are written assuming they will be run with `pytest -`_ using the Astropy `custom test runner -`_. To set up a -Conda environment to run the full set of tests, see the `environment-dev.yml -`_ environment -file. Once the dependencies are installed, you can run the tests two ways: +The tests are written assuming they will be run with `tox +`_ or `pytest `_. -1. By importing ``gala``:: +To run the tests with tox, first make sure that tox is installed; - import gala - gala.test() + pip install tox -2. By cloning the ``gala`` repository and running:: +then run the basic test suite with: - python setup.py test + tox -e test +or run the test suite with all optional dependencies with: + + tox -e test-alldeps + +You can see a list of available test environments with: + + tox -l -v + +which will also explain what each of them does. + +You can also run the tests directly with pytest. To do this, make sure to +install the testing requirements (from the cloned ``gala`` repository +directory):: + + pip install -e ".[test]" + +Then you can run the tests with: + + pytest gala diff --git a/docs/units.rst b/docs/units.rst index b408a9a9..8ad2d456 100644 --- a/docs/units.rst +++ b/docs/units.rst @@ -39,7 +39,7 @@ unit system using :meth:`~astropy.units.Quantity.decompose`:: `~gala.units.UnitSystem` objects can also act as a dictionary to look up a unit for a given physical type. For example, if we want to know what a 'velocity' -unit is in a given unit system, pass the key ``'velocity'``:: +unit is in a given unit system, pass the key ``'speed'`` or ``'velocity'``:: >>> usys['speed'] Unit("cm / ms") diff --git a/gala/coordinates/__init__.py b/gala/coordinates/__init__.py index 7cab629d..6f6dbc5c 100644 --- a/gala/coordinates/__init__.py +++ b/gala/coordinates/__init__.py @@ -11,4 +11,3 @@ from .reflex import * from .greatcircle import * from .pm_cov_transform import * -from .galactocentric import get_galactocentric2019 diff --git a/gala/coordinates/galactocentric.py b/gala/coordinates/galactocentric.py deleted file mode 100644 index bf7c071b..00000000 --- a/gala/coordinates/galactocentric.py +++ /dev/null @@ -1,21 +0,0 @@ -# Third-party -import astropy.coordinates as coord - -__all__ = ['get_galactocentric2019'] - - -def get_galactocentric2019(): - """ - DEPRECATED! - - Use astropy v4.0 or later: - - >>> with coord.galactocentric_frame_defaults.set('v4.0'): - ... galcen = coord.Galactocentric() - - """ - - with coord.galactocentric_frame_defaults.set('v4.0'): - galcen = coord.Galactocentric() - - return galcen diff --git a/gala/coordinates/reflex.py b/gala/coordinates/reflex.py index 5ad44bc1..30301e23 100644 --- a/gala/coordinates/reflex.py +++ b/gala/coordinates/reflex.py @@ -1,5 +1,4 @@ import astropy.coordinates as coord -from .galactocentric import get_galactocentric2019 __all__ = ['reflex_correct'] @@ -33,7 +32,7 @@ def reflex_correct(coords, galactocentric_frame=None): # If not specified, use the Astropy default Galactocentric frame if galactocentric_frame is None: - galactocentric_frame = get_galactocentric2019() + galactocentric_frame = coord.Galactocentric() v_sun = galactocentric_frame.galcen_v_sun diff --git a/gala/coordinates/tests/test_galactocentric.py b/gala/coordinates/tests/test_galactocentric.py deleted file mode 100644 index f96c3ddb..00000000 --- a/gala/coordinates/tests/test_galactocentric.py +++ /dev/null @@ -1,11 +0,0 @@ -# Third-party -import astropy.coordinates as coord - -# This package -from ..galactocentric import get_galactocentric2019 - - -def test_simple(): - # Not meant to be a sensible test - that's done in Astropy! - fr = get_galactocentric2019() - assert isinstance(fr, coord.Galactocentric) diff --git a/gala/coordinates/tests/test_greatcircle.py b/gala/coordinates/tests/test_greatcircle.py index 9e16bda5..6535abfa 100644 --- a/gala/coordinates/tests/test_greatcircle.py +++ b/gala/coordinates/tests/test_greatcircle.py @@ -7,7 +7,6 @@ # This project from ..greatcircle import (GreatCircleICRSFrame, make_greatcircle_cls, pole_from_endpoints, sph_midpoint) -from ..galactocentric import get_galactocentric2019 def test_cls_init(): @@ -32,7 +31,7 @@ def test_cls_init(): def test_init_center(): - galcen = get_galactocentric2019() + galcen = coord.Galactocentric() stupid_gal = GreatCircleICRSFrame( pole=coord.Galactic._ngp_J2000.transform_to(coord.ICRS()), center=galcen.galcen_coord) diff --git a/gala/coordinates/velocity_frame_transforms.py b/gala/coordinates/velocity_frame_transforms.py index 00aed6b2..0d90596c 100644 --- a/gala/coordinates/velocity_frame_transforms.py +++ b/gala/coordinates/velocity_frame_transforms.py @@ -2,8 +2,6 @@ import astropy.coordinates as coord -from .galactocentric import get_galactocentric2019 - __all__ = ["vgsr_to_vhel", "vhel_to_vgsr"] @@ -39,7 +37,7 @@ def vgsr_to_vhel(coordinate, vgsr, vsun=None): """ if vsun is None: - galcen = get_galactocentric2019() + galcen = coord.Galactocentric() vsun = galcen.galcen_v_sun.to_cartesian().xyz return vgsr - _get_vproj(coordinate, vsun) @@ -70,7 +68,7 @@ def vhel_to_vgsr(coordinate, vhel, vsun): """ if vsun is None: - galcen = get_galactocentric2019() + galcen = coord.Galactocentric() vsun = galcen.galcen_v_sun.to_cartesian().xyz return vhel + _get_vproj(coordinate, vsun) diff --git a/gala/dynamics/core.py b/gala/dynamics/core.py index 8570f343..944e4a30 100644 --- a/gala/dynamics/core.py +++ b/gala/dynamics/core.py @@ -802,11 +802,10 @@ def plot(self, components=None, units=None, auto_aspect=True, **kwargs): """ - try: - import matplotlib.pyplot as plt - except ImportError: - msg = 'matplotlib is required for visualization.' - raise ImportError(msg) + from gala.tests.optional_deps import HAS_MATPLOTLIB + if not HAS_MATPLOTLIB: + raise ImportError('matplotlib is required for visualization.') + import matplotlib.pyplot as plt if components is None: components = self.pos.components @@ -814,15 +813,10 @@ def plot(self, components=None, units=None, auto_aspect=True, **kwargs): x, labels = self._plot_prepare(components=components, units=units) - default_kwargs = { - 'marker': '.', - 'labels': labels, - 'plot_function': plt.scatter, - 'autolim': False - } - - for k, v in default_kwargs.items(): - kwargs[k] = kwargs.get(k, v) + kwargs.setdefault('marker', '.') + kwargs.setdefault('labels', labels) + kwargs.setdefault('plot_function', plt.scatter) + kwargs.setdefault('autolim', False) fig = plot_projections(x, **kwargs) diff --git a/gala/dynamics/orbit.py b/gala/dynamics/orbit.py index 1c25091d..6ea42da8 100644 --- a/gala/dynamics/orbit.py +++ b/gala/dynamics/orbit.py @@ -857,12 +857,10 @@ def plot(self, components=None, units=None, auto_aspect=True, **kwargs): fig : `~matplotlib.Figure` """ - - try: - import matplotlib.pyplot as plt - except ImportError: - msg = 'matplotlib is required for visualization.' - raise ImportError(msg) + from gala.tests.optional_deps import HAS_MATPLOTLIB + if not HAS_MATPLOTLIB: + raise ImportError('matplotlib is required for visualization.') + import matplotlib.pyplot as plt if components is None: if self.ndim == 1: # only a 1D orbit, so just plot time series @@ -873,15 +871,10 @@ def plot(self, components=None, units=None, auto_aspect=True, **kwargs): x, labels = self._plot_prepare(components=components, units=units) - default_kwargs = { - 'marker': '', - 'linestyle': '-', - 'labels': labels, - 'plot_function': plt.plot - } - - for k, v in default_kwargs.items(): - kwargs[k] = kwargs.get(k, v) + kwargs.setdefault('marker', '') + kwargs.setdefault('linestyle', '-') + kwargs.setdefault('labels', labels) + kwargs.setdefault('plot_function', plt.plot) fig = plot_projections(x, **kwargs) @@ -894,6 +887,96 @@ def plot(self, components=None, units=None, auto_aspect=True, **kwargs): return fig + def plot_3d(self, components=None, units=None, auto_aspect=True, + subplots_kwargs=None, **kwargs): + """ + Plot the specified 3D components. + + Parameters + ---------- + components : iterable (optional) + A list of component names (strings) to plot. By default, this is the + Cartesian positions ``['x', 'y', 'z']``. To plot Cartesian + velocities, pass in the velocity component names + ``['v_x', 'v_y', 'v_z']``. If the representation is different, the + component names will be different. For example, for a Cylindrical + representation, the components are ``['rho', 'phi', 'z']`` and + ``['v_rho', 'pm_phi', 'v_z']``. + units : `~astropy.units.UnitBase`, iterable (optional) + A single unit or list of units to display the components in. + auto_aspect : bool (optional) + Automatically enforce an equal aspect ratio. + ax : `matplotlib.axes.Axes` + The matplotlib Axes object to draw on. + subplots_kwargs : dict (optional) + Dictionary of kwargs passed to :func:`~matplotlib.pyplot.subplots`. + labels : iterable (optional) + List or iterable of axis labels as strings. They should correspond + to the dimensions of the input orbit. + plot_function : str (optional) + The ``matplotlib`` plot function to use. By default, this is 'plot' + but can also be, e.g., 'scatter'. + **kwargs + All other keyword arguments are passed to the ``plot_function``. + You can pass in any of the usual style kwargs like ``color=...``, + ``marker=...``, etc. + + Returns + ------- + fig : `~matplotlib.Figure` + + """ + from gala.tests.optional_deps import HAS_MATPLOTLIB + if not HAS_MATPLOTLIB: + raise ImportError('matplotlib is required for visualization.') + import matplotlib.pyplot as plt + from mpl_toolkits import mplot3d # noqa + + if components is None: + components = self.pos.components + + if subplots_kwargs is None: + subplots_kwargs = dict() + + if len(components) != 3: + raise ValueError( + f"The number of components ({len(components)}) must be 3") + + x, labels = self._plot_prepare(components=components, + units=units) + + kwargs.setdefault('marker', '') + kwargs.setdefault('linestyle', kwargs.pop('ls', '-')) + plot_function_name = kwargs.pop('plot_function', 'plot') + + ax = kwargs.pop('ax', None) + subplots_kwargs.setdefault('constrained_layout', True) + if ax is None: + fig, ax = plt.subplots(figsize=(6, 6), + subplot_kw=dict(projection='3d'), + **subplots_kwargs) + else: + fig = ax.figure + + plot_function = getattr(ax, plot_function_name) + if x[0].ndim > 1: + for n in range(x[0].shape[1]): + plot_function(*[xx[:, n] for xx in x], **kwargs) + else: + plot_function(*x, **kwargs) + ax.set_xlabel(labels[0]) + ax.set_ylabel(labels[1]) + ax.set_zlabel(labels[2]) + + if self.pos.get_name() == 'cartesian' and \ + all([not c.startswith('d_') for c in components]) and \ + 't' not in components and \ + auto_aspect: + for ax in fig.axes: + ax.set(aspect='auto', adjustable='datalim') + + return fig, ax + def to_frame(self, frame, current_frame=None, **kwargs): """ TODO: diff --git a/gala/dynamics/plot.py b/gala/dynamics/plot.py index c576ddf3..6a916045 100644 --- a/gala/dynamics/plot.py +++ b/gala/dynamics/plot.py @@ -4,7 +4,7 @@ __all__ = ['plot_projections'] -def _get_axes(dim, subplots_kwargs=dict()): +def _get_axes(dim, subplots_kwargs=None): """ Parameters ---------- @@ -13,16 +13,23 @@ def _get_axes(dim, subplots_kwargs=dict()): subplots_kwargs : dict (optional) Dictionary of kwargs passed to :func:`~matplotlib.pyplot.subplots`. """ - + from gala.tests.optional_deps import HAS_MATPLOTLIB + if not HAS_MATPLOTLIB: + raise ImportError('matplotlib is required for visualization.') import matplotlib.pyplot as plt + if subplots_kwargs is None: + subplots_kwargs = dict() + if dim > 1: n_panels = int(dim * (dim - 1) / 2) else: n_panels = 1 - figsize = subplots_kwargs.pop('figsize', (4*n_panels, 4)) - fig, axes = plt.subplots(1, n_panels, figsize=figsize, - **subplots_kwargs) + + subplots_kwargs.setdefault('figsize', (4*n_panels, 4)) + subplots_kwargs.setdefault('constrained_layout', True) + + fig, axes = plt.subplots(1, n_panels, **subplots_kwargs) if n_panels == 1: axes = [axes] @@ -34,7 +41,7 @@ def _get_axes(dim, subplots_kwargs=dict()): def plot_projections(x, relative_to=None, autolim=True, axes=None, - subplots_kwargs=dict(), labels=None, plot_function=None, + subplots_kwargs=None, labels=None, plot_function=None, **kwargs): """ Given N-dimensional quantity, ``x``, make a figure containing 2D projections @@ -118,5 +125,4 @@ def plot_projections(x, relative_to=None, autolim=True, axes=None, k += 1 - axes[0].figure.tight_layout() return axes[0].figure diff --git a/gala/units.py b/gala/units.py index 4924cbf1..f23b7a1c 100644 --- a/gala/units.py +++ b/gala/units.py @@ -1,4 +1,5 @@ -__all__ = ['UnitSystem', 'galactic', 'dimensionless', 'solarsystem'] +__all__ = ['UnitSystem', 'DimensionlessUnitSystem', + 'galactic', 'dimensionless', 'solarsystem'] # Third-party import astropy.units as u