1# -*- coding: utf-8 -*- 2"""Simple distance, velocity, and angle support for Skyfield. 3 4""" 5import numpy as np 6from numpy import abs, copysign, isnan 7from .constants import AU_KM, AU_M, C, DAY_S, tau 8from .descriptorlib import reify 9from .functions import _to_array, length_of 10 11_dfmt = '{0}{1:02}deg {2:02}\' {3:02}.{4:0{5}}"' 12_dsgn = '{0:+>1}{1:02}deg {2:02}\' {3:02}.{4:0{5}}"' 13_hfmt = '{0}{1:02}h {2:02}m {3:02}.{4:0{5}}s' 14 15class UnpackingError(Exception): 16 """You cannot iterate directly over a Skyfield measurement object.""" 17 18class Unit(object): 19 """A measurement that can be expressed in several choices of unit.""" 20 21 def __getitem__(self, *args): 22 cls = self.__class__ 23 name = cls.__name__ 24 s = 'to use this {0}, ask for its value in a particular unit:\n\n{1}' 25 attrs = sorted(k for k, v in cls.__dict__.items() 26 if k[0].islower() and isinstance(v, reify)) 27 examples = [' {0}.{1}'.format(name.lower(), attr) for attr in attrs] 28 raise UnpackingError(s.format(name, '\n'.join(examples))) 29 30 __iter__ = __getitem__ # give advice about both foo[i] and "x,y,z = foo" 31 32class Distance(Unit): 33 """A distance, stored internally as au and available in other units. 34 35 You can initialize a ``Distance`` by providing a single float or a 36 float array as either an ``au=``, ``km=``, or ``m=`` parameter. 37 38 You can access the magnitude of the distance with its three 39 attributes ``.au``, ``.km``, and ``.m``. By default a distance 40 prints itself in astronomical units (au), but you can take control 41 of the formatting and choice of units yourself using standard Python 42 numeric formatting: 43 44 >>> d = Distance(au=1) 45 >>> print(d) 46 1.0 au 47 >>> print('{:.2f} km'.format(d.km)) 48 149597870.70 km 49 50 """ 51 _warned = False 52 53 def __init__(self, au=None, km=None, m=None): 54 if au is not None: 55 self.au = _to_array(au) 56 """Astronomical units.""" 57 elif km is not None: 58 self.km = km = _to_array(km) 59 self.au = km / AU_KM 60 elif m is not None: 61 self.m = m = _to_array(m) 62 self.au = m / AU_M 63 else: 64 raise ValueError('to construct a Distance provide au, km, or m') 65 66 @classmethod 67 def from_au(cls, au): 68 self = cls.__new__(cls) 69 self.au = _to_array(au) 70 return self 71 72 @reify 73 def au(self): # Empty property to provide Sphinx docstring. 74 """Astronomical units (the Earth-Sun distance of 149,597,870,700 m).""" 75 76 @reify 77 def km(self): 78 """Kilometers (1,000 meters).""" 79 return self.au * AU_KM 80 81 @reify 82 def m(self): 83 """Meters.""" 84 return self.au * AU_M 85 86 def __str__(self): 87 n = self.au 88 return ('{0} au' if getattr(n, 'shape', 0) else '{0:.6} au').format(n) 89 90 def __repr__(self): 91 return '<{0} {1}>'.format(type(self).__name__, self) 92 93 def length(self): 94 """Compute the length when this is an |xyz| vector. 95 96 The Euclidean vector length of this vector is returned as a new 97 :class:`~skyfield.units.Distance` object. 98 99 >>> from skyfield.api import Distance 100 >>> d = Distance(au=[1, 1, 0]) 101 >>> d.length() 102 <Distance 1.41421 au> 103 104 """ 105 return Distance(au=length_of(self.au)) 106 107 def light_seconds(self): 108 """Return the length of this vector in light seconds.""" 109 return self.m / C 110 111 def to(self, unit): 112 """Convert this distance to the given AstroPy unit.""" 113 from astropy.units import au 114 return (self.au * au).to(unit) 115 116class Velocity(Unit): 117 """A velocity, stored internally as au/day and available in other units. 118 119 You can initialize a ``Velocity`` by providing a float or float 120 array to its ``au_per_d=`` parameter. 121 122 """ 123 _warned = False 124 125 # TODO: consider reworking this class to return a Rate object. 126 127 def __init__(self, au_per_d=None, km_per_s=None): 128 if km_per_s is not None: 129 self.km_per_s = km_per_s = _to_array(km_per_s) 130 self.au_per_d = km_per_s * DAY_S / AU_KM 131 elif au_per_d is not None: 132 self.au_per_d = _to_array(au_per_d) 133 else: 134 raise ValueError('to construct a Velocity provide' 135 ' au_per_d or km_per_s') 136 137 @reify 138 def au_per_d(self): # Empty property to provide Sphinx docstring. 139 """Astronomical units per day.""" 140 141 @reify 142 def km_per_s(self): 143 """Kilometers per second.""" 144 return self.au_per_d * AU_KM / DAY_S 145 146 @reify 147 def m_per_s(self): 148 """Meters per second.""" 149 return self.au_per_d * AU_M / DAY_S 150 151 def __str__(self): 152 n = self.au_per_d 153 fmt = '{0} au/day' if getattr(n, 'shape', 0) else '{0:.6} au/day' 154 return fmt.format(n) 155 156 def __repr__(self): 157 return '<{0} {1}>'.format(type(self).__name__, self) 158 159 def to(self, unit): 160 """Convert this velocity to the given AstroPy unit.""" 161 from astropy.units import au, d 162 return (self.au_per_d * au / d).to(unit) 163 164class AngleRate(object): 165 """The rate at which an angle is changing.""" 166 167 # TODO: design and implement public constructor. 168 169 @classmethod 170 def _from_radians_per_day(cls, radians_per_day): 171 ar = cls() 172 ar._radians_per_day = radians_per_day 173 return ar 174 175 @reify 176 def radians(self): 177 """:class:`Rate` of change in radians.""" 178 return Rate._from_per_day(self._radians_per_day) 179 180 @reify 181 def degrees(self): 182 """:class:`Rate` of change in degrees.""" 183 return Rate._from_per_day(self._radians_per_day / tau * 360.0) 184 185 @reify 186 def arcminutes(self): 187 """:class:`Rate` of change in arcminutes.""" 188 return Rate._from_per_day(self._radians_per_day / tau * 21600.0) 189 190 @reify 191 def arcseconds(self): 192 """:class:`Rate` of change in arcseconds.""" 193 return Rate._from_per_day(self._radians_per_day / tau * 1296000.0) 194 195 @reify 196 def mas(self): 197 """:class:`Rate` of change in milliarcseconds.""" 198 return Rate._from_per_day(self._radians_per_day / tau * 1.296e9) 199 200 # TODO: str; repr; conversion to AstroPy units 201 202class Rate(object): 203 """Measurement whose denominator is time.""" 204 205 # TODO: design and implement public constructor. 206 207 @classmethod 208 def _from_per_day(cls, per_day): 209 r = cls() 210 r._per_day = per_day 211 return r 212 213 @reify 214 def per_day(self): 215 """Units per day of Terrestrial Time.""" 216 return self._per_day 217 218 @reify 219 def per_hour(self): 220 """Units per hour of Terrestrial Time.""" 221 return self._per_day / 24.0 222 223 @reify 224 def per_minute(self): 225 """Units per minute of Terrestrial Time.""" 226 return self._per_day / 1440.0 227 228 @reify 229 def per_second(self): 230 """Units per second of Terrestrial Time.""" 231 return self._per_day / 86400.0 232 233# Angle units. 234 235_instantiation_instructions = """to instantiate an Angle, try one of: 236 237Angle(angle=another_angle) 238Angle(radians=value) 239Angle(degrees=value) 240Angle(hours=value) 241 242where `value` can be either a Python float, a list of Python floats, 243or a NumPy array of floats""" 244 245class Angle(Unit): 246 247 def __init__(self, angle=None, radians=None, degrees=None, hours=None, 248 preference=None, signed=False): 249 250 if angle is not None: 251 if not isinstance(angle, Angle): 252 raise ValueError(_instantiation_instructions) 253 self.radians = angle.radians 254 elif radians is not None: 255 self.radians = _to_array(radians) 256 elif degrees is not None: 257 self._degrees = degrees = _to_array(_unsexagesimalize(degrees)) 258 self.radians = degrees / 360.0 * tau 259 elif hours is not None: 260 self._hours = hours = _to_array(_unsexagesimalize(hours)) 261 self.radians = hours / 24.0 * tau 262 263 self.preference = (preference if preference is not None 264 else 'hours' if hours is not None 265 else 'degrees') 266 self.signed = signed 267 268 @classmethod 269 def from_degrees(cls, degrees, signed=False): 270 degrees = _to_array(_unsexagesimalize(degrees)) 271 self = cls.__new__(cls) 272 self.degrees = degrees 273 self.radians = degrees / 360.0 * tau 274 self.preference = 'degrees' 275 self.signed = signed 276 return self 277 278 @reify 279 def radians(self): # Empty property to provide Sphinx docstring. 280 """Radians (�� = 2�� in a circle).""" 281 282 @reify 283 def _hours(self): 284 return self.radians * 24.0 / tau 285 286 @reify 287 def _degrees(self): 288 return self.radians * 360.0 / tau 289 290 @reify 291 def hours(self): 292 """Hours (24\ |h| in a circle).""" 293 if self.preference != 'hours': 294 raise WrongUnitError('hours') 295 return self._hours 296 297 @reify 298 def degrees(self): 299 """Degrees (360° in a circle).""" 300 if self.preference != 'degrees': 301 raise WrongUnitError('degrees') 302 return self._degrees 303 304 def arcminutes(self): 305 """Return the angle in arcminutes.""" 306 return self._degrees * 60.0 307 308 def arcseconds(self): 309 """Return the angle in arcseconds.""" 310 return self._degrees * 3600.0 311 312 def mas(self): 313 """Return the angle in milliarcseconds.""" 314 return self._degrees * 3600000.0 315 316 def __str__(self): 317 size = self.radians.size 318 if size == 0: 319 return 'Angle []' 320 if self.preference == 'degrees': 321 v = self._degrees 322 fmt = _dsgn.format if self.signed else _dfmt.format 323 places = 1 324 else: 325 v = self._hours 326 fmt = _hfmt.format 327 places = 2 328 if size >= 2: 329 return '{0} values from {1} to {2}'.format( 330 len(v), _sfmt(fmt, v[0], places), _sfmt(fmt, v[-1], places)) 331 return _sfmt(fmt, v, places) 332 333 def __repr__(self): 334 if self.radians.size == 0: 335 return '<{0} []>'.format(type(self).__name__) 336 else: 337 return '<{0} {1}>'.format(type(self).__name__, self) 338 339 def hms(self, warn=True): 340 """Convert to a tuple (hours, minutes, seconds). 341 342 All three quantities will have the same sign as the angle itself. 343 344 """ 345 if warn and self.preference != 'hours': 346 raise WrongUnitError('hms') 347 sign, units, minutes, seconds = _sexagesimalize_to_float(self._hours) 348 return sign * units, sign * minutes, sign * seconds 349 350 def signed_hms(self, warn=True): 351 """Convert to a tuple (sign, hours, minutes, seconds). 352 353 The ``sign`` will be either +1 or -1, and the other quantities 354 will all be positive. 355 356 """ 357 if warn and self.preference != 'hours': 358 raise WrongUnitError('signed_hms') 359 return _sexagesimalize_to_float(self._hours) 360 361 def hstr(self, places=2, warn=True, format=_hfmt): 362 """Return a string like ``12h 07m 30.00s``; see `Formatting angles`. 363 364 .. versionadded:: 1.39 365 366 Added the ``format=`` parameter. 367 368 """ 369 if warn and self.preference != 'hours': 370 raise WrongUnitError('hstr') 371 hours = self._hours 372 shape = getattr(hours, 'shape', ()) 373 fmt = format.format # `format()` method of `format` string 374 if shape: 375 return [_sfmt(fmt, h, places) for h in hours] 376 return _sfmt(fmt, hours, places) 377 378 def dms(self, warn=True): 379 """Convert to a tuple (degrees, minutes, seconds). 380 381 All three quantities will have the same sign as the angle itself. 382 383 """ 384 if warn and self.preference != 'degrees': 385 raise WrongUnitError('dms') 386 sign, units, minutes, seconds = _sexagesimalize_to_float(self._degrees) 387 return sign * units, sign * minutes, sign * seconds 388 389 def signed_dms(self, warn=True): 390 """Convert to a tuple (sign, degrees, minutes, seconds). 391 392 The ``sign`` will be either +1 or -1, and the other quantities 393 will all be positive. 394 395 """ 396 if warn and self.preference != 'degrees': 397 raise WrongUnitError('signed_dms') 398 return _sexagesimalize_to_float(self._degrees) 399 400 def dstr(self, places=1, warn=True, format=None): 401 """Return a string like ``181deg 52' 30.0"``; see `Formatting angles`. 402 403 .. versionadded:: 1.39 404 405 Added the ``format=`` parameter. 406 407 """ 408 if warn and self.preference != 'degrees': 409 raise WrongUnitError('dstr') 410 degrees = self._degrees 411 signed = self.signed 412 if format is None: 413 format = _dsgn if signed else _dfmt 414 fmt = format.format # `format()` method of `format` string 415 shape = getattr(degrees, 'shape', ()) 416 if shape: 417 return [_sfmt(fmt, d, places) for d in degrees] 418 return _sfmt(fmt, degrees, places) 419 420 def to(self, unit): 421 """Convert this angle to the given AstroPy unit.""" 422 from astropy.units import rad 423 return (self.radians * rad).to(unit) 424 425 # Or should this do: 426 from astropy.coordinates import Angle 427 from astropy.units import rad 428 return Angle(self.radians, rad).to(unit) 429 430class WrongUnitError(ValueError): 431 432 def __init__(self, name): 433 unit = 'hours' if (name.startswith('h') or '_h' in name) else 'degrees' 434 usual = 'hours' if (unit == 'degrees') else 'degrees' 435 message = ('this angle is usually expressed in {0}, not {1};' 436 ' if you want to use {1} anyway,'.format(usual, unit)) 437 if name == unit: 438 message += ' then please use the attribute _{0}'.format(unit) 439 else: 440 message += ' then call {0}() with warn=False'.format(name) 441 self.args = (message,) 442 443def _sexagesimalize_to_float(value): 444 """Decompose `value` into units, minutes, and seconds. 445 446 Note that this routine is not appropriate for displaying a value, 447 because rounding to the smallest digit of display is necessary 448 before showing a value to the user. Use `_sexagesimalize_to_int()` 449 for data being displayed to the user. 450 451 This routine simply decomposes the floating point `value` into a 452 sign (+1.0 or -1.0), units, minutes, and seconds, returning the 453 result in a four-element tuple. 454 455 >>> _sexagesimalize_to_float(12.05125) 456 (1.0, 12.0, 3.0, 4.5) 457 >>> _sexagesimalize_to_float(-12.05125) 458 (-1.0, 12.0, 3.0, 4.5) 459 460 """ 461 sign = np.sign(value) 462 n = abs(value) 463 minutes, seconds = divmod(n * 3600.0, 60.0) 464 units, minutes = divmod(minutes, 60.0) 465 return sign, units, minutes, seconds 466 467def _sexagesimalize_to_int(value, places=0): 468 """Decompose `value` into units, minutes, seconds, and second fractions. 469 470 This routine prepares a value for sexagesimal display, with its 471 seconds fraction expressed as an integer with `places` digits. The 472 result is a tuple of five integers: 473 474 ``(sign [either +1 or -1], units, minutes, seconds, second_fractions)`` 475 476 The integers are properly rounded per astronomical convention so 477 that, for example, given ``places=3`` the result tuple ``(1, 11, 22, 478 33, 444)`` means that the input was closer to 11u 22' 33.444" than 479 to either 33.443" or 33.445" in its value. 480 481 """ 482 power = 10 ** places 483 n = int((power * 3600 * value + 0.5) // 1.0) 484 sign = np.sign(n) 485 n, fraction = divmod(abs(n), power) 486 n, seconds = divmod(n, 60) 487 n, minutes = divmod(n, 60) 488 return sign, n, minutes, seconds, fraction 489 490def _sfmt(fmt, value, places): 491 """Decompose floating point `value` into sexagesimal, and format.""" 492 if isnan(value): 493 return 'nan' 494 sgn, h, m, s, fraction = _sexagesimalize_to_int(value, places) 495 sign = '-' if sgn < 0.0 else '' 496 return fmt(sign, h, m, s, fraction, places) 497 498def wms(whole, minutes=0.0, seconds=0.0): 499 """Return a quantity expressed with 1/60 minutes and 1/3600 seconds.""" 500 return (whole 501 + copysign(minutes, whole) / 60.0 502 + copysign(seconds, whole) / 3600.0) 503 504def _unsexagesimalize(value): 505 """Return `value` after interpreting a (units, minutes, seconds) tuple. 506 507 When `value` is not a tuple, it is simply returned. 508 509 >>> _unsexagesimalize(3.25) 510 3.25 511 512 An input tuple is interpreted as units, minutes, and seconds. Note 513 that only the sign of `units` is significant! So all of the 514 following tuples convert into exactly the same value: 515 516 >>> '%f' % _unsexagesimalize((-1, 2, 3)) 517 '-1.034167' 518 >>> '%f' % _unsexagesimalize((-1, -2, 3)) 519 '-1.034167' 520 >>> '%f' % _unsexagesimalize((-1, -2, -3)) 521 '-1.034167' 522 523 """ 524 if isinstance(value, tuple): 525 components = iter(value) 526 value = next(components) 527 factor = 1.0 528 for component in components: 529 factor *= 60.0 530 value += copysign(component, value) / factor 531 return value 532 533def _interpret_angle(name, angle_object, angle_float, unit='degrees'): 534 """Return an angle in radians from one of two arguments. 535 536 It is common for Skyfield routines to accept both an argument like 537 `alt` that takes an Angle object as well as an `alt_degrees` that 538 can be given a bare float or a sexagesimal tuple. A pair of such 539 arguments can be passed to this routine for interpretation. 540 541 """ 542 if angle_object is not None: 543 if isinstance(angle_object, Angle): 544 return angle_object.radians 545 elif angle_float is not None: 546 return _unsexagesimalize(angle_float) / 360.0 * tau 547 raise ValueError('you must either provide the {0}= parameter with' 548 ' an Angle argument or supply the {0}_{1}= parameter' 549 ' with a numeric argument'.format(name, unit)) 550 551def _ltude(value, name, psuffix, nsuffix): 552 # Support for old deprecated Topos argument interpretation. 553 if not isinstance(value, str): 554 return _unsexagesimalize(value) 555 value = value.strip().upper() 556 if value.endswith(psuffix): 557 sign = +1.0 558 elif value.endswith(nsuffix): 559 sign = -1.0 560 else: 561 raise ValueError('your {0} string {1!r} does not end with either {2!r}' 562 ' or {3!r}'.format(name, value, psuffix, nsuffix)) 563 try: 564 value = float(value[:-1]) 565 except ValueError: 566 raise ValueError('your {0} string {1!r} cannot be parsed as a floating' 567 ' point number'.format(name, value)) 568 return sign * value 569