# -*- coding: utf-8 -*- """Simple distance, velocity, and angle support for Skyfield. """ import numpy as np from numpy import abs, copysign, isnan from types import MethodType from .constants import AU_KM, AU_M, C, DAY_S, tau from .descriptorlib import reify from .functions import _to_array, length_of _dfmt = '{0}{1:02}deg {2:02}\' {3:02}.{4:0{5}}"' _dsgn = '{0:+>1}{1:02}deg {2:02}\' {3:02}.{4:0{5}}"' _hfmt = '{0}{1:02}h {2:02}m {3:02}.{4:0{5}}s' class UnpackingError(Exception): """You cannot iterate directly over a Skyfield measurement object.""" class Unit(object): """A measurement that can be expressed in several choices of unit.""" def __getitem__(self, *args): """Tell users to ask for a specific unit before indexing or slicing.""" cls = self.__class__ name = cls.__name__ s = 'to use this {0}, ask for its value in a particular unit:\n\n{1}' attrs = sorted(k for k, v in cls.__dict__.items() if k[0].islower() and isinstance(v, (getset, reify))) examples = [' {0}.{1}'.format(name.lower(), attr) for attr in attrs] raise UnpackingError(s.format(name, '\n'.join(examples))) __iter__ = __getitem__ # give advice about both foo[i] and "x,y,z = foo" class getset(object): """Unit name that serves as both a class constructor and instance attribute. This supports two use cases: * When called as a class method like ``Distance.km(5.0)``, we build and return an instance of ``Distance`` whose ``km`` has been set to 5.0 and whose base unit ``m``, using the appropriate conversion factor, has been set to 5000.0. * When invoked like ``d.km`` on a particular ``Distance`` that doesn't yet have a ``km`` attribute (which otherwise Python itself would have returned), we apply the conversion factor to ``d.m`` and return the result. """ def __init__(self, name, docstring, conversion_factor=None, core_unit=None): self.name = name self.__doc__ = docstring self.conversion_factor = conversion_factor self.core_unit = core_unit if conversion_factor is None: def _constructor(cls, value): value = _to_array(value) obj = cls.__new__(cls) setattr(obj, name, value) return obj else: def _constructor(cls, value): value = _to_array(value) obj = cls.__new__(cls) setattr(obj, name, value) setattr(obj, core_unit, value / conversion_factor) return obj _constructor.__doc__ = self.__doc__ self._constructor = _constructor def __get__(self, instance, objtype=None): if instance is None: # the class itself has been asked for this name return MethodType(self._constructor, objtype) value = getattr(instance, self.core_unit) * self.conversion_factor instance.__dict__[self.name] = value return value class Distance(Unit): """A distance, stored internally as au and available in other units. You can initialize a ``Distance`` by providing a single float or a float array as either an ``au=``, ``km=``, or ``m=`` parameter. You can access the magnitude of the distance with its three attributes ``.au``, ``.km``, and ``.m``. By default a distance prints itself in astronomical units (au), but you can take control of the formatting and choice of units yourself using standard Python numeric formatting: >>> d = Distance(au=1) >>> print(d) 1.0 au >>> print('{:.2f} km'.format(d.km)) 149597870.70 km """ def __init__(self, au=None, km=None, m=None): if au is not None: self.au = _to_array(au) elif km is not None: self.km = km = _to_array(km) self.au = km / AU_KM elif m is not None: self.m = m = _to_array(m) self.au = m / AU_M else: raise ValueError('to construct a Distance provide au, km, or m') @classmethod def from_au(cls, au): # deprecated and no longer used internally return cls.au(au) au = getset('au', 'Astronomical units' ' (the Earth-Sun distance of 149,597,870,700 m).') km = getset('km', 'Kilometers (1,000 meters).', AU_KM, 'au') m = getset('m', 'Meters.', AU_M, 'au') def __str__(self): n = self.au return ('{0} au' if getattr(n, 'shape', 0) else '{0:.6} au').format(n) def __repr__(self): return '<{0} {1}>'.format(type(self).__name__, self) def length(self): """Compute the length when this is an |xyz| vector. The Euclidean vector length of this vector is returned as a new :class:`~skyfield.units.Distance` object. >>> from skyfield.api import Distance >>> d = Distance(au=[1, 1, 0]) >>> d.length() """ return Distance(au=length_of(self.au)) def light_seconds(self): """Return the length of this vector in light seconds.""" return self.m / C def to(self, unit): """Convert this distance to the given AstroPy unit.""" from astropy.units import au return (self.au * au).to(unit) class Velocity(Unit): """A velocity, stored internally as au/day and available in other units. You can initialize a ``Velocity`` by providing a float or float array to its ``au_per_d=`` parameter. """ # TODO: consider reworking this class to return a Rate object. def __init__(self, au_per_d=None, km_per_s=None): if km_per_s is not None: self.km_per_s = km_per_s = _to_array(km_per_s) self.au_per_d = km_per_s * DAY_S / AU_KM elif au_per_d is not None: self.au_per_d = _to_array(au_per_d) else: raise ValueError('to construct a Velocity provide' ' au_per_d or km_per_s') au_per_d = getset('au_per_d', 'Astronomical units per day.') km_per_s = getset('km_per_s', 'Kilometers per second.', AU_KM / DAY_S, 'au_per_d') m_per_s = getset('m_per_s', 'Meters per second.', AU_M / DAY_S, 'au_per_d') def __str__(self): n = self.au_per_d fmt = '{0} au/day' if getattr(n, 'shape', 0) else '{0:.6} au/day' return fmt.format(n) def __repr__(self): return '<{0} {1}>'.format(type(self).__name__, self) def to(self, unit): """Convert this velocity to the given AstroPy unit.""" from astropy.units import au, d return (self.au_per_d * au / d).to(unit) class AngleRate(object): """The rate at which an angle is changing.""" # TODO: design and implement public constructor. @classmethod def _from_radians_per_day(cls, radians_per_day): ar = cls() ar._radians_per_day = radians_per_day return ar @reify def radians(self): """:class:`Rate` of change in radians.""" return Rate._from_per_day(self._radians_per_day) @reify def degrees(self): """:class:`Rate` of change in degrees.""" return Rate._from_per_day(self._radians_per_day / tau * 360.0) @reify def arcminutes(self): """:class:`Rate` of change in arcminutes.""" return Rate._from_per_day(self._radians_per_day / tau * 21600.0) @reify def arcseconds(self): """:class:`Rate` of change in arcseconds.""" return Rate._from_per_day(self._radians_per_day / tau * 1296000.0) @reify def mas(self): """:class:`Rate` of change in milliarcseconds.""" return Rate._from_per_day(self._radians_per_day / tau * 1.296e9) # TODO: str; repr; conversion to AstroPy units class Rate(object): """Measurement whose denominator is time.""" # TODO: design and implement public constructor. @classmethod def _from_per_day(cls, per_day): r = cls() r._per_day = per_day return r @reify def per_day(self): """Units per day of Terrestrial Time.""" return self._per_day @reify def per_hour(self): """Units per hour of Terrestrial Time.""" return self._per_day / 24.0 @reify def per_minute(self): """Units per minute of Terrestrial Time.""" return self._per_day / 1440.0 @reify def per_second(self): """Units per second of Terrestrial Time.""" return self._per_day / 86400.0 # Angle units. _instantiation_instructions = """to instantiate an Angle, try one of: Angle(angle=another_angle) Angle(radians=value) Angle(degrees=value) Angle(hours=value) where `value` can be either a Python float, a list of Python floats, or a NumPy array of floats""" class Angle(Unit): def __init__(self, angle=None, radians=None, degrees=None, hours=None, preference=None, signed=False): if angle is not None: if not isinstance(angle, Angle): raise ValueError(_instantiation_instructions) self.radians = angle.radians elif radians is not None: self.radians = _to_array(radians) elif degrees is not None: self.degrees = degrees = _to_array(_unsexagesimalize(degrees)) self.radians = degrees / 360.0 * tau elif hours is not None: self.hours = hours = _to_array(_unsexagesimalize(hours)) self.radians = hours / 24.0 * tau self.preference = (preference if preference is not None else 'hours' if hours is not None else 'degrees') self.signed = signed @classmethod def from_degrees(cls, degrees, signed=False): degrees = _to_array(_unsexagesimalize(degrees)) self = cls.__new__(cls) self.degrees = degrees self.radians = degrees / 360.0 * tau self.preference = 'degrees' self.signed = signed return self radians = getset('radians', 'Radians (𝜏 = 2𝜋 in a circle).') # Deprecated names, to support legacy code. @reify def _hours(self): return self.hours @reify def _degrees(self): return self.degrees @reify def hours(self): r"""Hours (24\ |h| in a circle).""" return self.radians * 24.0 / tau @reify def degrees(self): """Degrees (360° in a circle).""" return self.radians * 360.0 / tau def arcminutes(self): """Return the angle in arcminutes.""" return self.degrees * 60.0 def arcseconds(self): """Return the angle in arcseconds.""" return self.degrees * 3600.0 def mas(self): """Return the angle in milliarcseconds.""" return self.degrees * 3600000.0 def __str__(self): size = self.radians.size if size == 0: return 'Angle []' if self.preference == 'degrees': v = self.degrees fmt = _dsgn.format if self.signed else _dfmt.format places = 1 else: v = self.hours fmt = _hfmt.format places = 2 if size >= 2: return '{0} values from {1} to {2}'.format( len(v), _sfmt(fmt, v[0], places), _sfmt(fmt, v[-1], places)) return _sfmt(fmt, v, places) def __repr__(self): if self.radians.size == 0: return '<{0} []>'.format(type(self).__name__) else: return '<{0} {1}>'.format(type(self).__name__, self) def hms(self, warn=True): """Convert to a tuple (hours, minutes, seconds). All three quantities will have the same sign as the angle itself. """ if warn and self.preference != 'hours': raise WrongUnitError('hms') sign, units, minutes, seconds = _sexagesimalize_to_float(self.hours) return sign * units, sign * minutes, sign * seconds def signed_hms(self, warn=True): """Convert to a tuple (sign, hours, minutes, seconds). The ``sign`` will be either +1 or -1, and the other quantities will all be positive. """ if warn and self.preference != 'hours': raise WrongUnitError('signed_hms') return _sexagesimalize_to_float(self.hours) def hstr(self, places=2, warn=True, format=_hfmt): """Return a string like ``12h 07m 30.00s``; see `Formatting angles`. .. versionadded:: 1.39 Added the ``format=`` parameter. """ if warn and self.preference != 'hours': raise WrongUnitError('hstr') hours = self.hours shape = getattr(hours, 'shape', ()) fmt = format.format # `format()` method of `format` string if shape: return [_sfmt(fmt, h, places) for h in hours] return _sfmt(fmt, hours, places) def dms(self, warn=True): """Convert to a tuple (degrees, minutes, seconds). All three quantities will have the same sign as the angle itself. """ if warn and self.preference != 'degrees': raise WrongUnitError('dms') sign, units, minutes, seconds = _sexagesimalize_to_float(self.degrees) return sign * units, sign * minutes, sign * seconds def signed_dms(self, warn=True): """Convert to a tuple (sign, degrees, minutes, seconds). The ``sign`` will be either +1 or -1, and the other quantities will all be positive. """ if warn and self.preference != 'degrees': raise WrongUnitError('signed_dms') return _sexagesimalize_to_float(self.degrees) def dstr(self, places=1, warn=True, format=None): """Return a string like ``181deg 52' 30.0"``; see `Formatting angles`. .. versionadded:: 1.39 Added the ``format=`` parameter. """ if warn and self.preference != 'degrees': raise WrongUnitError('dstr') degrees = self.degrees signed = self.signed if format is None: format = _dsgn if signed else _dfmt fmt = format.format # `format()` method of `format` string shape = getattr(degrees, 'shape', ()) if shape: return [_sfmt(fmt, d, places) for d in degrees] return _sfmt(fmt, degrees, places) def to(self, unit): """Convert this angle to the given AstroPy unit.""" from astropy.units import rad return (self.radians * rad).to(unit) # Or should this do: from astropy.coordinates import Angle from astropy.units import rad return Angle(self.radians, rad).to(unit) class WrongUnitError(ValueError): def __init__(self, name): unit = 'hours' if (name.startswith('h') or '_h' in name) else 'degrees' usual = 'hours' if (unit == 'degrees') else 'degrees' message = ('this angle is usually expressed in {0}, not {1};' ' if you want to use {1} anyway,'.format(usual, unit)) if name == unit: message += ' then please use the attribute _{0}'.format(unit) else: message += ' then call {0}() with warn=False'.format(name) self.args = (message,) def _sexagesimalize_to_float(value): """Decompose `value` into units, minutes, and seconds. Note that this routine is not appropriate for displaying a value, because rounding to the smallest digit of display is necessary before showing a value to the user. Use `_sexagesimalize_to_int()` for data being displayed to the user. This routine simply decomposes the floating point `value` into a sign (+1.0 or -1.0), units, minutes, and seconds, returning the result in a four-element tuple. >>> _sexagesimalize_to_float(12.05125) (1.0, 12.0, 3.0, 4.5) >>> _sexagesimalize_to_float(-12.05125) (-1.0, 12.0, 3.0, 4.5) """ sign = np.sign(value) n = abs(value) minutes, seconds = divmod(n * 3600.0, 60.0) units, minutes = divmod(minutes, 60.0) return sign, units, minutes, seconds def _sexagesimalize_to_int(value, places=0): """Decompose `value` into units, minutes, seconds, and second fractions. This routine prepares a value for sexagesimal display, with its seconds fraction expressed as an integer with `places` digits. The result is a tuple of five integers: ``(sign [either +1 or -1], units, minutes, seconds, second_fractions)`` The integers are properly rounded per astronomical convention so that, for example, given ``places=3`` the result tuple ``(1, 11, 22, 33, 444)`` means that the input was closer to 11u 22' 33.444" than to either 33.443" or 33.445" in its value. """ power = 10 ** places n = int((power * 3600 * value + 0.5) // 1.0) sign = np.sign(n) n, fraction = divmod(abs(n), power) n, seconds = divmod(n, 60) n, minutes = divmod(n, 60) return sign, n, minutes, seconds, fraction def _sfmt(fmt, value, places): """Decompose floating point `value` into sexagesimal, and format.""" if isnan(value): return 'nan' sgn, h, m, s, fraction = _sexagesimalize_to_int(value, places) sign = '-' if sgn < 0.0 else '' return fmt(sign, h, m, s, fraction, places) def wms(whole, minutes=0.0, seconds=0.0): """Return a quantity expressed with 1/60 minutes and 1/3600 seconds.""" return (whole + copysign(minutes, whole) / 60.0 + copysign(seconds, whole) / 3600.0) def _unsexagesimalize(value): """Return `value` after interpreting a (units, minutes, seconds) tuple. When `value` is not a tuple, it is simply returned. >>> _unsexagesimalize(3.25) 3.25 An input tuple is interpreted as units, minutes, and seconds. Note that only the sign of `units` is significant! So all of the following tuples convert into exactly the same value: >>> '%f' % _unsexagesimalize((-1, 2, 3)) '-1.034167' >>> '%f' % _unsexagesimalize((-1, -2, 3)) '-1.034167' >>> '%f' % _unsexagesimalize((-1, -2, -3)) '-1.034167' """ if isinstance(value, tuple): components = iter(value) value = next(components) factor = 1.0 for component in components: factor *= 60.0 value += copysign(component, value) / factor return value def _interpret_angle(name, angle_object, angle_float, unit='degrees'): """Return an angle in radians from one of two arguments. It is common for Skyfield routines to accept both an argument like `alt` that takes an Angle object as well as an `alt_degrees` that can be given a bare float or a sexagesimal tuple. A pair of such arguments can be passed to this routine for interpretation. """ if angle_object is not None: if isinstance(angle_object, Angle): return angle_object.radians elif angle_float is not None: return _unsexagesimalize(angle_float) / 360.0 * tau raise ValueError('you must either provide the {0}= parameter with' ' an Angle argument or supply the {0}_{1}= parameter' ' with a numeric argument'.format(name, unit)) def _ltude(value, name, psuffix, nsuffix): # Support for old deprecated Topos argument interpretation. if not isinstance(value, str): return _unsexagesimalize(value) value = value.strip().upper() if value.endswith(psuffix): sign = +1.0 elif value.endswith(nsuffix): sign = -1.0 else: raise ValueError('your {0} string {1!r} does not end with either {2!r}' ' or {3!r}'.format(name, value, psuffix, nsuffix)) try: value = float(value[:-1]) except ValueError: raise ValueError('your {0} string {1!r} cannot be parsed as a floating' ' point number'.format(name, value)) return sign * value