1# -*- coding: utf-8 -*- 2import datetime as dt_module 3import re 4import sys 5from collections import namedtuple 6from datetime import date, datetime, timedelta 7from numpy import ( 8 array, concatenate, cos, float_, int64, isnan, isinf, linspace, 9 nan, ndarray, nonzero, pi, rollaxis, searchsorted, sin, where, zeros_like, 10) 11from time import strftime, struct_time 12from ._compatibility import interp 13from .constants import ASEC2RAD, B1950, DAY_S, T0, tau 14from .curvelib import Splines, build_spline_given_ends 15from .descriptorlib import reify 16from .earthlib import sidereal_time, earth_rotation_angle 17from .framelib import ICRS_to_J2000 as B 18from .functions import (A, mxm, mxmxm, load_bundled_npy, rot_x, rot_y, rot_z, 19 _to_array, _reconcile) 20from .nutationlib import ( 21 build_nutation_matrix, equation_of_the_equinoxes_complimentary_terms, 22 iau2000a_radians, mean_obliquity, 23) 24from .precessionlib import compute_precession 25 26DAY_US = 86400000000.0 27GREGORIAN_START = 2299161 28GREGORIAN_START_ENGLAND = 2361222 29_OLD_PYTHON = sys.version_info < (2, 7) 30def _cat(*args): return concatenate(args, axis=1) 31 32CalendarTuple = namedtuple('CalendarTuple', 'year month day hour minute second') 33 34class CalendarArray(ndarray): 35 @property 36 def year(self): return self[0] 37 @property 38 def month(self): return self[1] 39 @property 40 def day(self): return self[2] 41 @property 42 def hour(self): return self[3] 43 @property 44 def minute(self): return self[4] 45 @property 46 def second(self): return self[5] 47 48if hasattr(dt_module, 'timezone'): 49 utc = dt_module.timezone.utc 50else: 51 class UTC(dt_module.tzinfo): 52 'UTC' 53 zero = timedelta(0) 54 def utcoffset(self, dt): 55 return self.zero 56 def tzname(self, dt): 57 return 'UTC' 58 def dst(self, dt): 59 return self.zero 60 61 utc = UTC() 62 63# Much of the following code is adapted from the USNO's "novas.c". 64 65_time_zero = dt_module.time(tzinfo=utc) 66MONTH_NAMES = A['0', 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 67 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] 68 69tt_minus_tai = array(32.184 / DAY_S) 70 71class Timescale(object): 72 """The data necessary to express dates in different timescales. 73 74 A `Timescale` provides time objects with the data tables they need 75 to translate between different time scales: the schedule of UTC leap 76 seconds, and the value of ∆T over time. Most programs create a 77 single `Timescale` which they use to build their `Time` objects: 78 79 >>> from skyfield.api import load 80 >>> ts = load.timescale() 81 >>> t = ts.utc(1980, 3, 1, 9, 30) 82 >>> t 83 <Time tt=2444299.896425741> 84 85 See :ref:`downloading-timescale-files` if you are interested in 86 checking how recent the data is in the files loaded by the 87 timescale. 88 89 """ 90 _utcnow = datetime.utcnow 91 polar_motion_table = None 92 93 def __init__(self, delta_t_recent, leap_dates, leap_offsets): 94 if callable(delta_t_recent): 95 # Let a caller completely override our approach to ∆T by 96 # passing a function of their own. 97 self.delta_t_function = delta_t_recent 98 else: 99 self.delta_t_table = delta_t_recent # so users can see it 100 self.delta_t_function = build_delta_t(delta_t_recent) 101 102 self.leap_dates, self.leap_offsets = leap_dates, leap_offsets 103 self.J2000 = Time(self, float_(T0)) 104 self.B1950 = Time(self, float_(B1950)) 105 self.julian_calendar_cutoff = None 106 107 # Our internal leap-second table has three columns: 108 # 109 # _leap_utc Integer UTC seconds since JD 0 w/leap seconds missing. 110 # _leap_offsets Integer offset between UTC and TAI. 111 # _leap_tai Integer TAI seconds since JD 0. 112 113 is_legacy_table = isinf(leap_dates[-1]) 114 if is_legacy_table: 115 leap_dates = leap_dates[2:-1] 116 leap_offsets = leap_offsets[3:] 117 118 one_zero = [[-1,0]] 119 self._leap_utc = (leap_dates[:,None] * DAY_S + one_zero).flatten() 120 self._leap_offsets = (leap_offsets[:,None] + one_zero).flatten() 121 self._leap_tai = self._leap_utc + self._leap_offsets 122 123 def now(self): 124 """Return the current date and time as a `Time` object. 125 126 For the return value to be correct, your operating system time 127 and timezone settings must be set so that the Python Standard 128 Library constructor ``datetime.datetime.utcnow()`` returns a 129 correct UTC date and time. 130 131 """ 132 return self.from_datetime(self._utcnow().replace(tzinfo=utc)) 133 134 def from_datetime(self, datetime): 135 """Return a `Time` for a Python ``datetime``. 136 137 The ``datetime`` must be “timezone-aware”: it must have a time 138 zone object as its ``tzinfo`` attribute instead of ``None``. 139 140 .. versionadded:: 1.24 141 142 """ 143 return self._utc(_datetime_to_utc_tuple(datetime)) 144 145 def from_datetimes(self, datetime_list): 146 """Return a `Time` for a list of Python ``datetime`` objects. 147 148 The ``datetime`` objects must each be “timezone-aware”: they 149 must each have a time zone object as their ``tzinfo`` attribute 150 instead of ``None``. 151 152 .. versionadded:: 1.24 153 154 """ 155 tuples = (_datetime_to_utc_tuple(d) for d in datetime_list) 156 return self._utc(array(value) for value in zip(*tuples)) 157 158 def utc(self, year, month=1, day=1, hour=0, minute=0, second=0.0): 159 """Build a `Time` from a UTC `calendar date`. 160 161 .. versionadded:: 1.24 162 Passing a Python ``datetime`` or a list of datetimes as the 163 first argument has been deprecated (and was never supported 164 for the other time scale methods). Instead, use the methods 165 :meth:`~skyfield.timelib.Timescale.from_datetime()` and 166 :meth:`~skyfield.timelib.Timescale.from_datetimes()`. 167 168 """ 169 # TODO: someday deprecate passing datetime objects here, as 170 # there are now separate constructors for them. 171 if isinstance(year, datetime): 172 return self.from_datetime(year) 173 if isinstance(year, date): 174 return self.from_datetime(datetime.combine(year, _time_zero)) 175 if hasattr(year, '__len__') and isinstance(year[0], datetime): 176 return self.from_datetimes(year) 177 178 tup = year, month, day, hour, minute, second 179 return self._utc(tup) 180 181 def _utc(self, tup): 182 # Build a Time from a UTC tuple, carefully preserving its exact 183 # second number in the Time's hidden TAI seconds field. 184 year, month, day, hour, minute, second = tup 185 a = _to_array 186 cutoff = self.julian_calendar_cutoff 187 188 # Figure out exactly the TAI second number. 189 seconds = (julian_day(a(year), a(month), a(day), cutoff) - 0.5) * DAY_S 190 seconds, sfr = divmod(seconds, 1.0) # in case there were any fractions 191 seconds += interp(seconds, self._leap_utc, self._leap_offsets) 192 more = a(hour) * 3600.0 + a(minute) * 60.0 + a(second) 193 seconds2, sfr = divmod(sfr + more, 1.0) 194 seconds += seconds2 195 196 # For the other timescales, use the usual Julian date + fraction. 197 whole, fraction = divmod(seconds, DAY_S) 198 fraction += sfr 199 fraction /= DAY_S 200 201 t = Time(self, whole, fraction + tt_minus_tai) 202 t.tai_fraction = fraction 203 t._tai_seconds = seconds, sfr 204 return t 205 206 def _jd(self, year, month, day, hour, minute, second): 207 a = _to_array 208 cutoff = self.julian_calendar_cutoff 209 whole = julian_day(a(year), a(month), a(day), cutoff) - 0.5 210 fraction = (a(second) + a(minute) * 60.0 + a(hour) * 3600.0) / DAY_S 211 return _reconcile(whole, fraction) 212 213 def _cal(self, whole, fraction): 214 return calendar_tuple(whole, fraction, self.julian_calendar_cutoff) 215 216 def _strftime(self, format, jd, fraction, seconds_bump=None): 217 # Python forces an unhappy choice upon us: either use the faster 218 # time.strftime() and lose support for '%f', or use the slower 219 # datetime.strftime() and crash if years are negative. We take the 220 # first option, but then patch '%f' support back in by secretly 221 # passing the microseconds string as the time zone name. After all, 222 # the routines supported by this function never use time zones. 223 # What could go wrong? 224 225 offset, ms = _strftime_offset_seconds(format) 226 fraction = fraction + offset / DAY_S 227 year, month, day, hour, minute, second = self._cal(jd, fraction) 228 z = year * 0 229 230 # TODO: will this always produce the same whole number that 231 # calendar_tuple() produces internally? Or should we make a private 232 # version of calendar_tuple() that returns it to us for use here? 233 weekday = (fraction + 0.5 + _to_array(jd)).astype(int) % 7 234 235 if _format_uses_day_of_year(format): 236 start_of_year = julian_day(year, 1, 1, self.julian_calendar_cutoff) 237 yday = (jd + fraction + 1.5 - start_of_year).astype(int) 238 else: 239 yday = z 240 241 if ms: 242 format = format[:ms.start()] + '%Z' + format[ms.end():] 243 second = (second * 1e6).astype(int) 244 second, usec = divmod(second, 1000000) 245 if seconds_bump is not None: 246 second += seconds_bump 247 if getattr(jd, 'ndim', 0): 248 u = ['%06d' % u for u in usec] 249 tup = (year, month, day, hour, minute, second, 250 weekday, yday, z, u) 251 return [strftime(format, struct_time(t)) for t in zip(*tup)] 252 u = '%06d' % usec 253 tup = year, month, day, hour, minute, second, weekday, yday, z, u 254 return strftime(format, struct_time(tup)) 255 else: 256 second = second.astype(int) 257 if seconds_bump is not None: 258 second += seconds_bump 259 tup = year, month, day, hour, minute, second, weekday, yday, z 260 if getattr(jd, 'ndim', 0): 261 return [strftime(format, item) for item in zip(*tup)] 262 return strftime(format, tup) 263 264 def tai(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0, 265 jd=None): 266 """Build a `Time` from an International Atomic Time `calendar date`. 267 268 .. versionadded:: 1.6 269 Passing a Julian date with ``jd=`` has been deprecated; 270 instead, use :meth:`~skyfield.timelib.Timescale.tai_jd()`. 271 272 """ 273 if jd is not None: 274 return self.tai_jd(jd) # deprecate someday 275 whole, fraction = self._jd(year, month, day, hour, minute, second) 276 t = Time(self, whole, fraction + tt_minus_tai) 277 t.tai_fraction = fraction 278 return t 279 280 def tai_jd(self, jd, fraction=None): 281 """Build a `Time` from an International Atomic Time Julian date.""" 282 jd, fraction = _normalize_jd_and_fraction(jd, fraction) 283 t = Time(self, jd, fraction + tt_minus_tai) 284 t.tai_fraction = fraction 285 return t 286 287 def tt(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0, 288 jd=None): 289 """Build a `Time` from a Terrestrial Time `calendar date`. 290 291 .. versionadded:: 1.6 292 Passing a Julian date with ``jd=`` has been deprecated; 293 instead, use :meth:`~skyfield.timelib.Timescale.tt_jd()`. 294 295 """ 296 if jd is not None: 297 return self.tt_jd(jd) # deprecate someday 298 whole, fraction = self._jd(year, month, day, hour, minute, second) 299 return Time(self, whole, fraction) 300 301 def tt_jd(self, jd, fraction=None): 302 """Build a `Time` from a Terrestrial Time Julian date.""" 303 jd, fraction = _normalize_jd_and_fraction(jd, fraction) 304 return Time(self, jd, fraction) 305 306 def J(self, year): 307 """Build a `Time` from a Terrestrial Time Julian year or array. 308 309 Julian years are convenient uniform periods of exactly 365.25 310 days of Terrestrial Time, centered on 2000 January 1 12h TT = 311 Julian year 2000.0. 312 313 """ 314 tt = _to_array(year) * 365.25 + 1721045.0 315 return Time(self, tt, 0.0) 316 317 def tdb(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0, 318 jd=None): 319 """Build a `Time` from a Barycentric Dynamical Time `calendar date`. 320 321 .. versionadded:: 1.6 322 Passing a Julian date with ``jd=`` has been deprecated; 323 instead, use :meth:`~skyfield.timelib.Timescale.tdb_jd()`. 324 325 """ 326 if jd is not None: 327 return self.tdb_jd(jd) # deprecate someday 328 whole, fraction = self._jd(year, month, day, hour, minute, second) 329 jd = whole + fraction # TODO: why do tests break if we pass separately 330 return Time(self, jd, - tdb_minus_tt(jd) / DAY_S) 331 332 def tdb_jd(self, jd, fraction=None): 333 """Build a `Time` from a Barycentric Dynamical Time Julian date.""" 334 jd, fraction = _normalize_jd_and_fraction(jd, fraction) 335 t = Time(self, jd, fraction - tdb_minus_tt(jd, fraction) / DAY_S) 336 t.tdb_fraction = fraction 337 return t 338 339 def ut1(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0, 340 jd=None): 341 """Build a `Time` from a UT1 Universal Time `calendar date`. 342 343 .. versionadded:: 1.6 344 Passing a Julian date with ``jd=`` has been deprecated; 345 instead, use :meth:`~skyfield.timelib.Timescale.ut1_jd()`. 346 347 """ 348 if jd is None: # TODO: deprecate the jd parameter to this method 349 whole, fraction = self._jd(year, month, day, hour, minute, second) 350 jd = whole + fraction # TODO: can we pass high precision on? 351 return self.ut1_jd(jd) 352 353 def ut1_jd(self, jd): 354 """Build a `Time` from a UT1 Universal Time Julian date.""" 355 ut1 = _to_array(jd) 356 357 # Estimate TT = UT1, to get a rough Delta T estimate. 358 tt_approx = ut1 359 delta_t_approx = self.delta_t_function(tt_approx) 360 361 # Use the rough Delta T to make a much better estimate of TT, 362 # then generate an even better Delta T. 363 tt_approx = ut1 + delta_t_approx / DAY_S 364 delta_t_approx = self.delta_t_function(tt_approx) 365 366 # We can now estimate TT with an error of < 1e-9 seconds within 367 # 10 centuries of either side of the present; for details, see: 368 # https://github.com/skyfielders/astronomy-notebooks 369 # and look for the notebook "error-in-timescale-ut1.ipynb". 370 delta_t_approx /= DAY_S 371 t = Time(self, ut1, delta_t_approx) 372 t.ut1_fraction = 0.0 * ut1 373 return t 374 375 def from_astropy(self, t): 376 """Build a Skyfield `Time` from an AstroPy time object.""" 377 return self.tt(jd=t.tt.jd) 378 379 def linspace(self, t0, t1, num=50): 380 """Return ``num`` times spaced uniformly between ``t0`` to ``t1``. 381 382 This routine is named after, and powered by, the NumPy routine 383 `linspace()`_. 384 385 .. _linspace(): https://numpy.org/doc/stable/reference/generated/numpy.linspace.html 386 387 """ 388 whole0 = t0.whole 389 frac0 = t0.tt_fraction 390 whole1 = t1.whole 391 frac1 = t1.tt_fraction 392 return Time( 393 self, 394 linspace(whole0, whole1, num), 395 linspace(frac0, frac1, num), 396 ) 397 398class Time(object): 399 """A single moment in history, or an array of several moments. 400 401 Skyfield programs don’t usually instantiate this class directly, but 402 instead build time objects using one of the timescale methods listed 403 at `timescale-summary`. If you do attempt the low-level operation 404 of building a time object yourself, either leave ``tt_fraction`` at 405 its default value of ``None`` — in which case Skyfield will assume 406 the fraction is zero — or provide a ``tt_fraction`` array that has 407 exactly the same dimensions as your ``tt`` array. 408 409 """ 410 def __init__(self, ts, tt, tt_fraction=None): 411 if tt_fraction is None: 412 tt_fraction = zeros_like(tt) 413 self.ts = ts 414 self.whole = tt 415 self.tt_fraction = tt_fraction 416 self.shape = getattr(tt, 'shape', ()) 417 418 def __len__(self): 419 return self.shape[0] 420 421 def __repr__(self): 422 size = getattr(self.tt, 'size', -1) 423 if size > 3: 424 rstr = '<Time tt=[{0} ... {1}] len={2}>' 425 return rstr.format(self.tt[0], self.tt[-1], size) 426 else: 427 return ('<Time tt={0}>'.format(self.tt) 428 .replace('[ ', '[').replace(' ', ' ')) 429 430 def __getitem__(self, index): 431 # TODO: also copy cached matrices? 432 # TODO: raise non-IndexError exception if this Time is not an array; 433 # otherwise, a `for` loop over it will not raise an error. 434 t = Time(self.ts, self.whole[index], self.tt_fraction[index]) 435 d = self.__dict__ 436 for name in 'tai_fraction', 'tdb_fraction', 'ut1_fraction': 437 value = d.get(name) 438 if value is not None: 439 setattr(t, name, value[index]) 440 return t 441 442 def astimezone(self, tz): 443 """Convert to a Python ``datetime`` in a particular timezone ``tz``. 444 445 If this time is an array, then an array of datetimes is returned 446 instead of a single value. 447 448 """ 449 dt, leap_second = self.astimezone_and_leap_second(tz) 450 return dt 451 452 def astimezone_and_leap_second(self, tz): 453 """Convert to a Python ``datetime`` and leap second in a timezone. 454 455 Convert this time to a Python ``datetime`` and a leap second:: 456 457 dt, leap_second = t.astimezone_and_leap_second(tz) 458 459 The argument ``tz`` should be a ``datetime`` compatible 460 timezone. 461 462 The leap second value is provided because a Python ``datetime`` 463 can only number seconds ``0`` through ``59``, but leap seconds 464 have a designation of at least ``60``. The leap second return 465 value will normally be ``0``, but will instead be ``1`` if the 466 date and time are a UTC leap second. Add the leap second value 467 to the ``second`` field of the ``datetime`` to learn the real 468 name of the second. 469 470 If this time is an array, then an array of ``datetime`` objects 471 and an array of leap second integers is returned, instead of a 472 single value each. 473 474 """ 475 dt, leap_second = self.utc_datetime_and_leap_second() 476 normalize = getattr(tz, 'normalize', None) 477 if self.shape and normalize is not None: 478 dt = array([normalize(d.astimezone(tz)) for d in dt]) 479 elif self.shape: 480 dt = array([d.astimezone(tz) for d in dt]) 481 elif normalize is not None: 482 dt = normalize(dt.astimezone(tz)) 483 else: 484 dt = dt.astimezone(tz) 485 return dt, leap_second 486 487 def toordinal(self): 488 """Return the proleptic Gregorian ordinal of the UTC date. 489 490 This method makes Skyfield `Time` objects compatible with Python 491 `datetime`_ objects, which also provide a ``toordinal()`` 492 method. Thanks to this method, a `Time` can often be used 493 directly as a coordinate for a plot. 494 495 """ 496 whole, fraction, is_leap_second = self._utc_seconds(0.0) 497 return (whole - 148731076800.0 + fraction) / DAY_S 498 499 def utc_datetime(self): 500 """Convert to a Python ``datetime`` in UTC. 501 502 If this time is an array, then a list of datetimes is returned 503 instead of a single value. 504 505 """ 506 dt, leap_second = self.utc_datetime_and_leap_second() 507 return dt 508 509 def utc_datetime_and_leap_second(self): 510 """Convert to a Python ``datetime`` in UTC, plus a leap second value. 511 512 Convert this time to a `datetime`_ object and a leap second:: 513 514 dt, leap_second = t.utc_datetime_and_leap_second() 515 516 The leap second value is provided because a Python ``datetime`` 517 can only number seconds ``0`` through ``59``, but leap seconds 518 have a designation of at least ``60``. The leap second return 519 value will normally be ``0``, but will instead be ``1`` if the 520 date and time are a UTC leap second. Add the leap second value 521 to the ``second`` field of the ``datetime`` to learn the real 522 name of the second. 523 524 If this time is an array, then an array of ``datetime`` objects 525 and an array of leap second integers is returned, instead of a 526 single value each. 527 528 """ 529 year, month, day, hour, minute, second = self._utc_tuple(0.5e-6) 530 micro = (second * 1e6).astype(int) 531 second, micro = divmod(micro, 1000000) 532 leap_second = second // 60 533 second -= leap_second # fit within limited bounds of Python datetime 534 if self.shape: 535 zone = [utc] * self.shape[0] 536 argsets = zip(year, month, day, hour, minute, second, micro, zone) 537 dt = array([datetime(*args) for args in argsets]) 538 else: 539 dt = datetime(year, month, day, hour, minute, second, micro, utc) 540 return dt, leap_second 541 542 def utc_iso(self, delimiter='T', places=0): 543 """Convert to an ISO 8601 string like ``2014-01-18T01:35:38Z`` in UTC. 544 545 If this time is an array of dates, then a sequence of strings is 546 returned instead of a single string. 547 548 """ 549 # "places" used to be the 1st argument, so continue to allow an 550 # integer in that spot. TODO: deprecate this in Skyfield 2.0 551 # and remove it in 3.0. 552 if isinstance(delimiter, int): 553 places = delimiter 554 delimiter = 'T' 555 556 if places: 557 power_of_ten = 10 ** places 558 offset = 0.5 / power_of_ten 559 year, month, day, hour, minute, second = self._utc_tuple(offset) 560 second, fraction = divmod(second, 1.0) 561 fraction *= power_of_ten 562 format = '%04d-%02d-%02d{0}%02d:%02d:%02d.%0{1}dZ'.format( 563 delimiter, places) 564 args = (year, month, day, hour, minute, second, fraction) 565 else: 566 format = '%04d-%02d-%02d{0}%02d:%02d:%02dZ'.format(delimiter) 567 args = self._utc_tuple(0.5) 568 569 if self.shape: 570 return [format % tup for tup in zip(*args)] 571 else: 572 return format % args 573 574 def utc_jpl(self): 575 """Convert to a string like ``A.D. 2014-Jan-18 01:35:37.5000 UTC``. 576 577 Returns a string for this date and time in UTC, in the format 578 used by the JPL HORIZONS system. If this time is an array of 579 dates, then a sequence of strings is returned instead of a 580 single string. 581 582 """ 583 year, month, day, hour, minute, second = self._utc_tuple(0.00005) 584 second, fraction = divmod(second, 1.0) 585 fraction *= 1e4 586 bc = year < 1 587 year = abs(year - bc) 588 era = where(bc, 'B.C.', 'A.D.') 589 format = '%s %04d-%s-%02d %02d:%02d:%02d.%04d UTC' 590 args = (era, year, MONTH_NAMES[month], day, 591 hour, minute, second, fraction) 592 593 if self.shape: 594 return [format % tup for tup in zip(*args)] 595 else: 596 return format % args 597 598 def utc_strftime(self, format='%Y-%m-%d %H:%M:%S UTC'): 599 """Format the UTC time using a Python datetime formatting string. 600 601 This calls Python’s ``time.strftime()`` to format the date and 602 time. A single string is returned or else a whole array of 603 strings, depending on whether this time object is an array. 604 The most commonly used formats are: 605 606 * ``%Y`` four-digit year, ``%y`` two-digit year 607 * ``%m`` month number, ``%B`` name, ``%b`` abbreviation 608 * ``%d`` day of month 609 * ``%H`` hour 610 * ``%M`` minute 611 * ``%S`` second 612 * ``%A`` day of week, ``%a`` its abbreviation 613 614 You can find the full list, along with options that control 615 field widths and leading zeros, at: 616 617 https://docs.python.org/3/library/time.html#time.strftime 618 619 If the smallest time unit in your format is minutes or seconds, 620 then the time is rounded to the nearest minute or second. 621 Otherwise the value is truncated rather than rounded. 622 623 """ 624 offset, uses_ms = _strftime_offset_seconds(format) 625 year, month, day, hour, minute, second, jd = self._utc_tuple(offset, 1) 626 start_of_year = julian_day(year, 1, 1, self.ts.julian_calendar_cutoff) 627 weekday = jd % 7 628 yday = jd + 1 - start_of_year 629 return _strftime(format, year, month, day, hour, minute, second, 630 weekday, yday, uses_ms) 631 632 def _utc_tuple(self, offset, return_jd=False): 633 """Return UTC as (year, month, day, hour, minute, second.fraction). 634 635 The `offset` in seconds is added to the UTC time before it is 636 split into its components. This is useful if the user is going 637 to round the result before displaying it. If the result is 638 going to be displayed as seconds, for example, set `offset` to 639 0.5 and then throw away the fraction; if the result is going to 640 be displayed as minutes, set `offset` to 30.0 and then throw 641 away the seconds; and so forth. 642 643 """ 644 second, sfr, is_leap_second = self._utc_seconds(offset) 645 second = second.astype(int64) 646 second -= is_leap_second 647 jd, second = divmod(second + 43200, 86400) 648 cutoff = self.ts.julian_calendar_cutoff 649 year, month, day = compute_calendar_date(jd, cutoff) 650 minute, second = divmod(second, 60) 651 hour, minute = divmod(minute, 60) 652 second += is_leap_second 653 if not return_jd: 654 return year, month, day, hour, minute, second + sfr 655 return year, month, day, hour, minute, second + sfr, jd 656 657 def _utc_seconds(self, offset): 658 """Return integer seconds since JD 0.0, plus a 0 ≤ fraction < 1.""" 659 seconds, fr = self._tai_seconds 660 seconds2, fr = divmod(fr + offset, 1.0) 661 seconds = seconds + seconds2 # not +=, which would modify cached array 662 ts = self.ts 663 tai_minus_utc = interp(seconds, ts._leap_tai, ts._leap_offsets) 664 tai_minus_utc, is_leap_second = divmod(tai_minus_utc, 1.0) 665 is_leap_second = is_leap_second > 0.0 666 return seconds - tai_minus_utc, fr, is_leap_second 667 668 @reify 669 def _tai_seconds(self): 670 # If UTC was supplied, this will already hold an exact value; otherwise: 671 seconds, fr = divmod(self.whole * DAY_S, 1.0) 672 seconds2, fr = divmod(fr + self.tai_fraction * DAY_S, 1.0) 673 seconds += seconds2 674 return seconds, fr 675 676 def _leap_seconds(self): 677 # TODO: should this be reified? 678 ts = self.ts 679 seconds, fr = self._tai_seconds 680 return interp(seconds, ts._leap_tai, ts._leap_offsets) 681 682 # Calendar tuples. 683 684 def tai_calendar(self): 685 """TAI as a (year, month, day, hour, minute, second) `calendar date`.""" 686 return self.ts._cal(self.whole, self.tai_fraction) 687 688 def tt_calendar(self): 689 """TT as a (year, month, day, hour, minute, second) `calendar date`.""" 690 return self.ts._cal(self.whole, self.tt_fraction) 691 692 def tdb_calendar(self): 693 """TDB as a (year, month, day, hour, minute, second) `calendar date`.""" 694 return self.ts._cal(self.whole, self.tdb_fraction) 695 696 def ut1_calendar(self): 697 """UT1 as a (year, month, day, hour, minute, second) `calendar date`.""" 698 return self.ts._cal(self.whole, self.ut1_fraction) 699 700 # Date formatting. 701 702 def tai_strftime(self, format='%Y-%m-%d %H:%M:%S TAI'): 703 """Format TAI with a datetime strftime() format string.""" 704 return self.ts._strftime(format, self.whole, self.tai_fraction) 705 706 def tt_strftime(self, format='%Y-%m-%d %H:%M:%S TT'): 707 """Format TT with a datetime strftime() format string.""" 708 return self.ts._strftime(format, self.whole, self.tt_fraction) 709 710 def tdb_strftime(self, format='%Y-%m-%d %H:%M:%S TDB'): 711 """Format TDB with a datetime strftime() format string.""" 712 return self.ts._strftime(format, self.whole, self.tdb_fraction) 713 714 def ut1_strftime(self, format='%Y-%m-%d %H:%M:%S UT1'): 715 """Format UT1 with a datetime strftime() format string.""" 716 return self.ts._strftime(format, self.whole, self.ut1_fraction) 717 718 # Convenient caching of several expensive functions of time. 719 720 @reify 721 def M(self): 722 """3×3 rotation matrix: ICRS → equinox of this date.""" 723 724 # Compute N and P instead of asking for self.N and self.P to 725 # avoid keeping copies of them, since Skyfield itself never uses 726 # them again once M has been computed. But we do check in case 727 # a user has forced them to be built already. 728 729 d = self.__dict__ 730 731 P = d.get('P') 732 if P is None: 733 P = self.precession_matrix() 734 735 N = d.get('N') 736 if N is None: 737 N = self.nutation_matrix() 738 739 return mxmxm(N, P, B) 740 741 @reify 742 def MT(self): 743 """3×3 rotation matrix: equinox of this date → ICRS.""" 744 return rollaxis(self.M, 1) 745 746 @reify 747 def C(self): 748 # Calculate the Equation of Origins in cycles 749 eq_origins = (earth_rotation_angle(self.ut1) - self.gast / 24.0) 750 R = rot_z(2 * pi * eq_origins) 751 return mxm(R, self.M) 752 753 @reify 754 def CT(self): 755 return rollaxis(self.C, 1) 756 757 @reify 758 def _nutation_angles_radians(self): 759 # TODO: add psi and eps corrections support back in here, rather 760 # than at points of use. 761 return iau2000a_radians(self) 762 763 def _nutation_angles(self, angles): 764 # Sample code shared with early adopters suggested that setting 765 # this attribute manually could avoid the expense of IAU 2000A, 766 # so this setter continues to support the pattern. 767 768 d_psi, d_eps = angles 769 self._nutation_angles_radians = ( 770 d_psi / 1e7 * ASEC2RAD, 771 d_eps / 1e7 * ASEC2RAD, 772 ) 773 774 _nutation_angles = property(None, _nutation_angles) 775 776 @reify 777 def _mean_obliquity_radians(self): 778 # Cached because it is used to compute both gast and N. 779 return mean_obliquity(self.tdb) * ASEC2RAD 780 781 # Conversion between timescales. 782 783 @reify 784 def J(self): 785 """Return a floating point Julian year or array of years for this date. 786 787 Julian years are convenient uniform periods of exactly 365.25 788 days of Terrestrial Time, centered on 2000 January 1 12h TT = 789 Julian year 2000.0. 790 791 """ 792 return (self.whole - 1721045.0 + self.tt_fraction) / 365.25 793 794 @reify 795 def utc(self): 796 """A tuple ``(year, month, day, hour, minute, seconds)`` in UTC.""" 797 utc = self._utc_tuple(0.0) 798 return (array(utc).view(CalendarArray) if self.shape 799 else CalendarTuple(*utc)) 800 801 @reify 802 def tai_fraction(self): 803 return self.tt_fraction - tt_minus_tai 804 805 @reify 806 def tdb_fraction(self): 807 fr = self.tt_fraction 808 return fr + tdb_minus_tt(self.whole, fr) / DAY_S 809 810 @reify 811 def ut1_fraction(self): 812 return self.tt_fraction - self.delta_t / DAY_S 813 814 @reify 815 def delta_t(self): 816 return self.ts.delta_t_function(self.tt) 817 818 @reify 819 def dut1(self): 820 return 32.184 + self._leap_seconds() - self.delta_t 821 822 @reify 823 def gmst(self): 824 """Greenwich Mean Sidereal Time (GMST) in hours.""" 825 return sidereal_time(self) 826 827 @reify 828 def gast(self): 829 """Greenwich Apparent Sidereal Time (GAST) in hours.""" 830 d_psi, _ = self._nutation_angles_radians 831 tt = self.tt 832 # TODO: move this into an eqeq function? 833 c_terms = equation_of_the_equinoxes_complimentary_terms(tt) 834 eq_eq = d_psi * cos(self._mean_obliquity_radians) + c_terms 835 return (self.gmst + eq_eq / tau * 24.0) % 24.0 836 837 # Low-precision floats generated from internal float pairs. 838 839 @property 840 def tai(self): 841 return self.whole + self.tai_fraction 842 843 @property 844 def tt(self): 845 return self.whole + self.tt_fraction 846 847 @property 848 def tdb(self): 849 return self.whole + self.tdb_fraction 850 851 @property 852 def ut1(self): 853 return self.whole + self.ut1_fraction 854 855 # Earth position as a function of time. 856 857 def polar_motion_angles(self): 858 table = self.ts.polar_motion_table 859 if table is None: 860 return 0.0, 0.0, 0.0 861 sprime = -47.0e-6 * (self.tdb - T0) / 36525.0 862 tt, x, y = table 863 return sprime, interp(self.tt, tt, x), interp(self.tt, tt, y) 864 865 def polar_motion_matrix(self): 866 sprime, x, y = self.polar_motion_angles() 867 return mxmxm( 868 rot_x(y * ASEC2RAD), 869 rot_y(x * ASEC2RAD), 870 rot_z(-sprime * ASEC2RAD), 871 ) 872 873 def nutation_matrix(self): 874 """Compute the 3×3 nutation matrix N for this date.""" 875 d_psi, d_eps = self._nutation_angles_radians 876 mean_obliquity = self._mean_obliquity_radians 877 true_obliquity = mean_obliquity + d_eps 878 return build_nutation_matrix(mean_obliquity, true_obliquity, d_psi) 879 880 def precession_matrix(self): 881 """Compute the 3×3 precession matrix P for this date.""" 882 return compute_precession(self.tdb) 883 884 # Various dunders. 885 886 def __eq__(self, other_time): 887 if isinstance(other_time, Time): 888 return self.__sub__(other_time) == 0.0 889 return False 890 891 def __add__(self, other_time): 892 if isinstance(other_time, timedelta): 893 w = other_time.days 894 f = other_time.seconds / DAY_S + other_time.microseconds / DAY_US 895 elif isinstance(other_time, (int, float)): 896 w, f = divmod(other_time, 1.0) 897 else: 898 return NotImplemented 899 900 return self.ts.tt_jd(self.whole + w, self.tt_fraction + f) 901 902 def __sub__(self, other_time): 903 if isinstance(other_time, Time): 904 return self.whole - other_time.whole + ( 905 self.tt_fraction - other_time.tt_fraction 906 ) 907 elif isinstance(other_time, timedelta): 908 w = other_time.days 909 f = other_time.seconds / DAY_S + other_time.microseconds / DAY_US 910 elif isinstance(other_time, (int, float)): 911 w, f = divmod(other_time, 1.0) 912 else: 913 return NotImplemented 914 915 return self.ts.tt_jd(self.whole - w, self.tt_fraction - f) 916 917 def __hash__(self): 918 # Someone wanted to use Time objects with functools.lru_cache so 919 # we make this attempt to support hashability; beware that it 920 # will return the same hash for very closely spaced times that 921 # all round to the same floating point TT. 922 return hash(self.tt) 923 924 # Deprecated attributes that were once used internally, consuming 925 # memory with matrices that are never used again by Skyfield once 926 # t.M has been computed. 927 928 P = reify(precession_matrix) 929 N = reify(nutation_matrix) 930 P.__doc__ = N.__doc__ = None # omit from Sphinx documentation 931 932 @reify 933 def PT(self): return rollaxis(self.P, 1) 934 @reify 935 def NT(self): return rollaxis(self.N, 1) 936 937def julian_day(year, month=1, day=1, julian_before=None): 938 """Given a calendar date, return a Julian day integer. 939 940 Uses the proleptic Gregorian calendar unless ``julian_before`` is 941 set to a specific Julian day, in which case the Julian calendar is 942 used for dates older than that. 943 944 """ 945 # Support months <1 and >12 by overflowing cleanly into adjacent years. 946 y, month = divmod(month - 1, 12) 947 year = year + y 948 month += 1 949 950 # See the Explanatory Supplement to the Astronomical Almanac 15.11. 951 janfeb = month <= 2 952 g = year + 4716 - janfeb 953 f = (month + 9) % 12 954 e = 1461 * g // 4 + day - 1402 955 J = e + (153 * f + 2) // 5 956 957 mask = 1 if (julian_before is None) else (J >= julian_before) 958 J += (38 - (g + 184) // 100 * 3 // 4) * mask 959 return J 960 961def julian_date(year, month=1, day=1, hour=0, minute=0, second=0.0): 962 """Given a proleptic Gregorian calendar date and time, build a Julian date. 963 964 The difference between a “Julian day” and a “Julian date” is that 965 the “day” is the integer part, while the “date” includes a fraction 966 indicating the time. 967 968 """ 969 return julian_day(year, month, day) - 0.5 + ( 970 second + minute * 60.0 + hour * 3600.0) / DAY_S 971 972def julian_date_of_besselian_epoch(b): 973 return 2415020.31352 + (b - 1900.0) * 365.242198781 974 975def compute_calendar_date(jd_integer, julian_before=None): 976 """Convert Julian day ``jd_integer`` into a calendar (year, month, day). 977 978 Uses the proleptic Gregorian calendar unless ``julian_before`` is 979 set to a specific Julian day, in which case the Julian calendar is 980 used for dates older than that. 981 982 """ 983 use_gregorian = (julian_before is None) or (jd_integer >= julian_before) 984 985 # See the Explanatory Supplement to the Astronomical Almanac 15.11. 986 f = jd_integer + 1401 987 f += use_gregorian * ((4 * jd_integer + 274277) // 146097 * 3 // 4 - 38) 988 e = 4 * f + 3 989 g = e % 1461 // 4 990 h = 5 * g + 2 991 day = h % 153 // 5 + 1 992 month = (h // 153 + 2) % 12 + 1 993 year = e // 1461 - 4716 + (12 + 2 - month) // 12 994 return year, month, day 995 996calendar_date = compute_calendar_date # old name, in case anyone used it 997 998def calendar_tuple(jd_float, fraction=0.0, julian_before=None): 999 """Return a (year, month, day, hour, minute, second.fraction) tuple.""" 1000 jd_float = _to_array(jd_float) 1001 whole1, fraction1 = divmod(jd_float, 1.0) 1002 whole2, fraction = divmod(fraction1 + fraction + 0.5, 1.0) 1003 whole = (whole1 + whole2).astype(int) 1004 year, month, day = compute_calendar_date(whole, julian_before) 1005 second = fraction * 86400.0 1006 minute, second = divmod(second, 60.0) 1007 minute = minute.astype(int) 1008 hour, minute = divmod(minute, 60) 1009 return year, month, day, hour, minute, second 1010 1011def tdb_minus_tt(jd_tdb, fraction_tdb=0.0): 1012 """Computes how far TDB is in advance of TT, given TDB. 1013 1014 Given that the two time scales never diverge by more than 2ms, TT 1015 can also be given as the argument to perform the conversion in the 1016 other direction. 1017 1018 """ 1019 t = (jd_tdb - T0 + fraction_tdb) / 36525.0 1020 1021 # USNO Circular 179, eq. 2.6. 1022 return (0.001657 * sin ( 628.3076 * t + 6.2401) 1023 + 0.000022 * sin ( 575.3385 * t + 4.2970) 1024 + 0.000014 * sin (1256.6152 * t + 6.1969) 1025 + 0.000005 * sin ( 606.9777 * t + 4.0212) 1026 + 0.000005 * sin ( 52.9691 * t + 0.4444) 1027 + 0.000002 * sin ( 21.3299 * t + 5.5431) 1028 + 0.000010 * t * sin ( 628.3076 * t + 4.2490)) 1029 1030class DeltaT(object): 1031 def __init__(self, table_tt, table_delta_t, long_term_function): 1032 self.table_tt = table_tt 1033 self.table_delta_t = table_delta_t 1034 self.long_term_function = long_term_function 1035 1036 def __call__(self, tt): 1037 delta_t = interp(tt, self.table_tt, self.table_delta_t, nan, nan) 1038 if getattr(delta_t, 'shape', None): 1039 nan_indexes = nonzero(isnan(delta_t)) 1040 if nan_indexes: 1041 J = (tt[nan_indexes] - 1721045.0) / 365.25 1042 delta_t[nan_indexes] = self.long_term_function(J) 1043 else: 1044 if isnan(delta_t): 1045 J = (tt - 1721045.0) / 365.25 1046 delta_t = self.long_term_function(J) 1047 return delta_t 1048 1049delta_t_parabola_stephenson_morrison_hohenkerk_2016 = Splines( 1050 [1825.0, 1925.0, 0.0, 32.5, 0.0, -320.0]) 1051 1052delta_t_parabola_morrison_stephenson_2004 = Splines( 1053 [1820.0, 1920.0, 0.0, 32.0, 0.0, -20.0]) 1054 1055def build_delta_t(delta_t_recent): 1056 parabola = delta_t_parabola_stephenson_morrison_hohenkerk_2016 1057 s15_table = load_bundled_npy('delta_t.npz')['Table-S15.2020.txt'] 1058 table_tt, table_delta_t = delta_t_recent 1059 1060 p = parabola 1061 pd = parabola.derivative 1062 s = Splines(s15_table) 1063 sd = s.derivative 1064 1065 long_term_parabola_width = p.upper[0] - p.lower[0] 1066 1067 # How many years wide we make the splines that connect the tables to 1068 # the long-term parabola; tuned by hand until the derivative graphed 1069 # by `work_on_delta_t_discontinuities()` in `design/delta_t.py` 1070 # doesn't look too terrible. 1071 patch_width = 800.0 1072 1073 # To the left of the Table-S15 splines, design a spline connecting 1074 # them to the long-term parabola. 1075 1076 x1 = s.lower[0] # For the current table, this = -720.0. 1077 x0 = x1 - patch_width 1078 left = build_spline_given_ends(x0, p(x0), pd(x0), x1, s(x1), sd(x1)) 1079 1080 # And to the left of that, put the pure long-term parabola. 1081 1082 x1 = x0 1083 x0 = x1 - long_term_parabola_width 1084 far_left = build_spline_given_ends(x0, p(x0), pd(x0), x1, p(x1), pd(x1)) 1085 1086 # Truncate the splines table where the daily table starts, and 1087 # adjust the final spline to remove any discontinuity. 1088 1089 x = (table_tt[0] - 1721045.0) / 365.25 # TT to J centuries 1090 1091 i = searchsorted(s15_table[0], x) 1092 s15_table = s15_table[:,:i] 1093 1094 desired_y = table_delta_t[0] 1095 current_y = s(x) 1096 x0, x1, a3, a2, a1, a0 = s15_table[:,-1] 1097 t = (x - x0) / (x1 - x0) 1098 a1 = a1 + (desired_y - current_y) / t # adjust linear term 1099 s15_table[:,-1] = x0, x1, a3, a2, a1, a0 1100 1101 # To the right of the recent ∆T table, design a spline connecting 1102 # smoothly to the long-term parabola. 1103 1104 x0 = (table_tt[-1] - 1721045.0) / 365.25 # TT to J centuries 1105 x1 = (x0 + patch_width) // 100.0 * 100.0 # Choose multiple of 100 years 1106 y0 = table_delta_t[-1] 1107 1108 lookback = min(366, len(table_delta_t)) # Slope of last year of ∆T. 1109 slope = (table_delta_t[-1] - table_delta_t[-lookback]) * lookback / 365.0 1110 right = build_spline_given_ends(x0, y0, slope, x1, p(x1), pd(x1)) 1111 1112 # At the far right, finish with the pure long-term parabola. 1113 1114 x0 = x1 1115 x1 = x0 + long_term_parabola_width 1116 far_right = build_spline_given_ends(x0, p(x0), pd(x0), x1, p(x1), pd(x1)) 1117 1118 curve = Splines(_cat( 1119 array([far_left]).T, 1120 array([left]).T, 1121 s15_table, 1122 array([right]).T, 1123 array([far_right]).T, 1124 )) 1125 1126 return DeltaT(table_tt, table_delta_t, curve) 1127 1128def build_delta_t_table(delta_t_recent): 1129 """Deprecated: the Delta T interpolation table used in Skyfield <= 1.37.""" 1130 bundled = _cat( 1131 load_bundled_npy('morrison_stephenson_deltat.npy')[:,:22], 1132 load_bundled_npy('historic_deltat.npy'), 1133 ) 1134 recent_start_time = delta_t_recent[0,0] 1135 i = searchsorted(bundled[0], recent_start_time) 1136 century = 36524.0 1137 start_tt = bundled[0,0] - century 1138 start_J = Time(None, start_tt).J 1139 end_tt = delta_t_recent[0,-1] + century 1140 end_J = Time(None, end_tt).J 1141 return _cat( 1142 [[start_tt], [delta_t_parabola_morrison_stephenson_2004(start_J)]], 1143 bundled[:,:i], 1144 delta_t_recent, 1145 [[end_tt], [delta_t_parabola_morrison_stephenson_2004(end_J)]], 1146 ) 1147 1148_format_uses_milliseconds = re.compile(r'%[-_0^#EO]*f').search 1149_format_uses_seconds = re.compile(r'%[-_0^#EO]*[STXc]').search 1150_format_uses_minutes = re.compile(r'%[-_0^#EO]*[MR]').search 1151_format_uses_day_of_year = re.compile(r'%[-_0^#EO]*j').search 1152 1153def _datetime_to_utc_tuple(dt): 1154 z = dt.tzinfo 1155 if z is None: 1156 raise ValueError(_naive_complaint) 1157 if z is not utc: 1158 dt = dt.astimezone(utc) 1159 return (dt.year, dt.month, dt.day, 1160 dt.hour, dt.minute, dt.second + dt.microsecond / 1e6) 1161 1162def _normalize_jd_and_fraction(jd, fraction): 1163 jd = _to_array(jd) 1164 if fraction is None: 1165 jd, fraction = divmod(jd, 1.0) 1166 else: 1167 jd, fraction = _reconcile(jd, _to_array(fraction)) 1168 return jd, fraction 1169 1170def _strftime_offset_seconds(format): 1171 uses_ms = _format_uses_milliseconds(format) 1172 if uses_ms: 1173 if _OLD_PYTHON: 1174 raise ValueError('strftime() "%f" not supported under Python 2') 1175 offset = 1e-16 # encourage .0 to not turn into .999999 1176 elif _format_uses_seconds(format): 1177 offset = 0.5 1178 elif _format_uses_minutes(format): 1179 offset = 30.0 1180 else: 1181 offset = 0.0 1182 return offset, uses_ms 1183 1184def _strftime(format, year, month, day, hour, minute, second, 1185 weekday, yday, uses_ms): 1186 zero = year * 0 1187 1188 # TODO: avoid computing yday if this is false: 1189 #_format_uses_day_of_year(format) 1190 1191 if uses_ms: 1192 format = format[:uses_ms.start()] + '%Z' + format[uses_ms.end():] 1193 second = (second * 1e6).astype(int) 1194 second, usec = divmod(second, 1000000) 1195 if getattr(year, 'ndim', 0): 1196 u = ['%06d' % u for u in usec] 1197 tup = (year, month, day, hour, minute, second, 1198 weekday, yday, zero, u) 1199 return [strftime(format, struct_time(t)) for t in zip(*tup)] 1200 u = '%06d' % usec 1201 tup = year, month, day, hour, minute, second, weekday, yday, zero, u 1202 return strftime(format, struct_time(tup)) 1203 else: 1204 second = second.astype(int) 1205 tup = year, month, day, hour, minute, second, weekday, yday, zero 1206 if getattr(year, 'ndim', 0): 1207 return [strftime(format, item) for item in zip(*tup)] 1208 return strftime(format, tup) 1209 1210_naive_complaint = """cannot interpret a datetime that lacks a timezone 1211 1212You must either specify that your datetime is in UTC: 1213 1214 from skyfield.api import utc 1215 d = datetime(..., tzinfo=utc) # to build a new datetime 1216 d = d.replace(tzinfo=utc) # to fix an existing datetime 1217 1218Or use a timezone object like those provided by the third-party `pytz` library: 1219 1220 from pytz import timezone 1221 eastern = timezone('US/Eastern') 1222 d = eastern.localize(datetime(2014, 1, 16, 1, 32, 9))""" 1223