astropy:docs

Source code for astropy.coordinates.distances

# Licensed under a 3-clause BSD style license - see LICENSE.rst

"""
This module contains the classes and utility functions for distance and
cartesian coordinates.
"""

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import math

import numpy as np

from .. import units as u

__all__ = ['Distance', 'CartesianPoints', 'cartesian_to_spherical',
           'spherical_to_cartesian']


__doctest_requires__ = {'*': ['scipy.integrate']}


[docs]class Distance(u.Quantity): """ A one-dimensional distance. This can be initialized in one of four ways: * A distance ``value`` (array or float) and a ``unit`` * A `~astropy.units.Quantity` object * A redshift and (optionally) a cosmology. * Providing a distance modulus Parameters ---------- value : scalar or `~astropy.units.Quantity`. The value of this distance. unit : `~astropy.units.UnitBase` The units for this distance, *if* ``value`` is not a `~astropy.units.Quantity`. Must have dimensions of distance. z : float A redshift for this distance. It will be converted to a distance by computing the luminosity distance for this redshift given the cosmology specified by ``cosmology``. Must be given as a keyword argument. cosmology : ``Cosmology`` or `None` A cosmology that will be used to compute the distance from ``z``. If `None`, the current cosmology will be used (see `astropy.cosmology` for details). distmod : float or `~astropy.units.Quantity` The distance modulus for this distance. dtype : `~numpy.dtype`, optional See `~astropy.units.Quantity`. Must be given as a keyword argument. copy : bool, optional See `~astropy.units.Quantity`. Must be given as a keyword argument. Raises ------ `~astropy.units.UnitsError` If the ``unit`` is not a distance. ValueError If value specified is less than 0 and ``allow_negative=False``. If ``z`` is provided with a ``unit`` or ``cosmology`` is provided when ``z`` is *not* given, or ``value`` is given as well as ``z``. Examples -------- >>> from astropy import units as u >>> from astropy import cosmology >>> from astropy.cosmology import WMAP5, WMAP7 >>> cosmology.set_current(WMAP7) >>> d1 = Distance(10, u.Mpc) >>> d2 = Distance(40, unit=u.au) >>> d3 = Distance(value=5, unit=u.kpc) >>> d4 = Distance(z=0.23) >>> d5 = Distance(z=0.23, cosmology=WMAP5) >>> d6 = Distance(distmod=24.47) >>> d7 = Distance(Distance(10 * u.Mpc)) """ _include_easy_conversion_members = True def __new__(cls, value=None, unit=None, z=None, cosmology=None, distmod=None, dtype=None, copy=True, allow_negative=False): from ..cosmology import default_cosmology if isinstance(value, u.Quantity): # This includes Distances as well if z is not None or distmod is not None: raise ValueError('`value` was given along with `z` or `distmod`' ' in Quantity constructor.') if unit is not None: value = value.to(unit).value else: unit = value.unit value = value.value elif value is None: if z is not None: if distmod is not None: raise ValueError('Both `z` and `distmod` given in Distance ' 'constructor') if cosmology is None: cosmology = default_cosmology.get() ld = cosmology.luminosity_distance(z) if unit is not None: ld = ld.to(unit) value = ld.value unit = ld.unit elif distmod is not None: value = cls._distmod_to_pc(distmod) if unit is None: # choose unit based on most reasonable of Mpc, kpc, or pc if value > 1e6: value = value / 1e6 unit = u.megaparsec elif value > 1e3: value = value / 1e3 unit = u.kiloparsec else: unit = u.parsec else: value = u.Quantity(value, u.parsec).to(unit).value else: raise ValueError('None of `value`, `z`, or `distmod` were given' ' to Distance constructor') value = ld.value unit = ld.unit elif z is not None: # and value is not None based on above raise ValueError('Both `z` and a `value` were provided in Distance ' 'constructor') elif cosmology is not None: raise ValueError('A `cosmology` was given but `z` was not provided ' 'in Distance constructor') elif unit is None: raise u.UnitsError('No unit was provided for Distance') #"else" the baseline ``value`` + ``unit`` case unit = _convert_to_and_validate_length_unit(unit) try: value = np.asarray(value) except ValueError as e: raise TypeError(str(e)) if value.dtype.kind not in 'iuf': raise TypeError("Unsupported dtype '{0}'".format(value.dtype)) if np.any(value < 0) and not allow_negative: raise ValueError("Distance must be >= 0. Set the kwarg " "'allow_negative=True' to allow negative values.") return super(Distance, cls).__new__(cls, value, unit, dtype=dtype, copy=copy) def __quantity_subclass__(self, unit): if unit.is_equivalent(u.m): return Distance, True else: return super(Distance, self).__quantity_subclass__(unit)[0], False @property def z(self): """Short for ``self.compute_z()``""" return self.compute_z()
[docs] def compute_z(self, cosmology=None): """ The redshift for this distance assuming its physical distance is a luminosity distance. Parameters ---------- cosmology : ``Cosmology`` or `None` The cosmology to assume for this calculation, or `None` to use the current cosmology (see `astropy.cosmology` for details). Returns ------- z : float The redshift of this distance given the provided ``cosmology``. """ from ..cosmology import luminosity_distance from scipy import optimize # FIXME: array: need to make this calculation more vector-friendly f = lambda z, d, cos: (luminosity_distance(z, cos).value - d) ** 2 return optimize.brent(f, (self.Mpc, cosmology))
@property def distmod(self): """ The distance modulus of this distance as a `~astropy.units.Quantity` """ val = 5. * np.log10(self.to(u.pc).value) - 5. return u.Quantity(val, u.mag) @staticmethod def _distmod_to_pc(dm): return 10 ** ((dm + 5) / 5.)
[docs]class CartesianPoints(u.Quantity): """ A cartesian representation of a point in three-dimensional space. Parameters ---------- x : `~astropy.units.Quantity` or array-like The first cartesian coordinate or a single array or `~astropy.units.Quantity` where the first dimension is length-3. y : `~astropy.units.Quantity` or array-like, optional The second cartesian coordinate. z : `~astropy.units.Quantity` or array-like, optional The third cartesian coordinate. unit : `~astropy.units.UnitBase` object or `None` The physical unit of the coordinate values. If ``x``, ``y``, or ``z`` are quantities, they will be converted to this unit. dtype : `~numpy.dtype`, optional See `~astropy.units.Quantity`. Must be given as a keyword argument. copy : bool, optional See `~astropy.units.Quantity`. Must be given as a keyword argument. Raises ------ UnitsError If the units on ``x``, ``y``, and ``z`` do not match or an invalid unit is given. ValueError If ``y`` and ``z`` don't match ``x``'s shape or ``x`` is not length-3 TypeError If incompatible array types are passed into ``x``, ``y``, or ``z`` """ #this ensures that __array_wrap__ gets called for ufuncs even when #where a quantity is first, like ``3*u.m + c`` __array_priority__ = 10001 def __new__(cls, x, y=None, z=None, unit=None, dtype=None, copy=True): if y is None and z is None: if len(x) != 3: raise ValueError('Input to CartesianPoints is not length 3') qarr = x if unit is None and hasattr(qarr, 'unit'): unit = qarr.unit # for when a Quantity is given elif y is not None and z is not None: if unit is None: #they must all match units or this fails for coo in (x, y, z): if hasattr(coo, 'unit'): if unit is not None and coo.unit != unit: raise u.UnitsError('Units for `x`, `y`, and `z` do ' 'not match in CartesianPoints ') unit = coo.unit #if `unit` is still None at this point, it means none were #Quantties, which is fine, because it means the user wanted #the unit to be None else: #convert any that are like a Quantity to the given unit if hasattr(x, 'to'): x = x.to(unit) if hasattr(y, 'to'): y = y.to(unit) if hasattr(z, 'to'): z = z.to(unit) qarr = [np.asarray(coo) for coo in (x, y, z)] if not (qarr[0].shape == qarr[1].shape == qarr[2].shape): raise ValueError("Shapes for `x`, `y`, and `z` don't match in " "CartesianPoints") #let the unit be whatever it is else: raise TypeError('Must give all of `x`, `y`, and `z` or just array in ' 'CartesianPoints') try: unit = _convert_to_and_validate_length_unit(unit, True) except TypeError as e: raise u.UnitsError(str(e)) try: qarr = np.asarray(qarr) except ValueError as e: raise TypeError(str(e)) if qarr.dtype.kind not in 'iuf': raise TypeError("Unsupported dtype '{0}'".format(qarr.dtype)) return super(CartesianPoints, cls).__new__(cls, qarr, unit, dtype=dtype, copy=copy) def __quantity_subclass__(self, unit): if unit.is_equivalent(u.m): return CartesianPoints, True else: return super(CartesianPoints, self).__quantity_subclass__(unit)[0], False def __array_wrap__(self, obj, context=None): #always convert to CartesianPoints because all operations that would #screw up the units are killed by _convert_to_and_validate_length_unit obj = super(CartesianPoints, self).__array_wrap__(obj, context=context) #always prefer self's unit, if possible if obj.unit.is_equivalent(self.unit): return obj.to(self.unit) else: return obj @property def x(self): """ The second cartesian coordinate as a `~astropy.units.Quantity`. """ return self.view(u.Quantity)[0] @property def y(self): """ The second cartesian coordinate as a `~astropy.units.Quantity`. """ return self.view(u.Quantity)[1] @property def z(self): """ The third cartesian coordinate as a `~astropy.units.Quantity`. """ return self.view(u.Quantity)[2]
[docs] def to_spherical(self): """ Converts to the spherical representation of this point. Returns ------- r : `~astropy.units.Quantity` The radial coordinate (in the same units as this `CartesianPoints`). lat : `~astropy.units.Quantity` The spherical coordinates latitude. lon : `~astropy.units.Quantity` The spherical coordinates longitude. """ from .angles import Latitude, Longitude rarr, latarr, lonarr = cartesian_to_spherical(self.x, self.y, self.z) r = Distance(rarr, unit=self.unit) lat = Latitude(latarr, unit=u.radian) lon = Longitude(lonarr, unit=u.radian) return r, lat, lon
def _convert_to_and_validate_length_unit(unit, allow_dimensionless=False): """ raises UnitsError if not a length unit """ try: unit = u.Unit(unit) assert (unit.is_equivalent(u.kpc) or allow_dimensionless and unit == u.dimensionless_unscaled) except (TypeError, AssertionError): raise u.UnitsError('Unit "{0}" is not a length type'.format(unit)) return unit #<------------transformation-related utility functions----------------->
[docs]def cartesian_to_spherical(x, y, z): """ Converts 3D rectangular cartesian coordinates to spherical polar coordinates. Note that the resulting angles are latitude/longitude or elevation/azimuthal form. I.e., the origin is along the equator rather than at the north pole. .. note:: This is a low-level function used internally in `astropy.coordinates`. It is provided for users if they really want to use it, but it is recommended that you use the `astropy.coordinates` coordinate systems. Parameters ---------- x : scalar or array-like The first cartesian coordinate. y : scalar or array-like The second cartesian coordinate. z : scalar or array-like The third cartesian coordinate. Returns ------- r : float or array The radial coordinate (in the same units as the inputs). lat : float or array The latitude in radians lon : float or array The longitude in radians """ xsq = x ** 2 ysq = y ** 2 zsq = z ** 2 r = (xsq + ysq + zsq) ** 0.5 s = (xsq + ysq) ** 0.5 if np.isscalar(x) and np.isscalar(y) and np.isscalar(z): lon = math.atan2(y, x) lat = math.atan2(z, s) else: lon = np.arctan2(y, x) lat = np.arctan2(z, s) return r, lat, lon
[docs]def spherical_to_cartesian(r, lat, lon): """ Converts spherical polar coordinates to rectangular cartesian coordinates. Note that the input angles should be in latitude/longitude or elevation/azimuthal form. I.e., the origin is along the equator rather than at the north pole. .. note:: This is a low-level function used internally in `astropy.coordinates`. It is provided for users if they really want to use it, but it is recommended that you use the `astropy.coordinates` coordinate systems. Parameters ---------- r : scalar or array-like The radial coordinate (in the same units as the inputs). lat : scalar or array-like The latitude in radians lon : scalar or array-like The longitude in radians Returns ------- x : float or array The first cartesian coordinate. y : float or array The second cartesian coordinate. z : float or array The third cartesian coordinate. """ if np.isscalar(r) and np.isscalar(lat) and np.isscalar(lon): x = r * math.cos(lat) * math.cos(lon) y = r * math.cos(lat) * math.sin(lon) z = r * math.sin(lat) else: x = r * np.cos(lat) * np.cos(lon) y = r * np.cos(lat) * np.sin(lon) z = r * np.sin(lat) return x, y, z

Page Contents