1# -*- coding: utf-8 -*- 2 3 4# PyMeeus: Python module implementing astronomical algorithms. 5# Copyright (C) 2018 Dagoberto Salazar 6# 7# This program is free software: you can redistribute it and/or modify 8# it under the terms of the GNU Lesser General Public License as published by 9# the Free Software Foundation, either version 3 of the License, or 10# (at your option) any later version. 11# 12# This program is distributed in the hope that it will be useful, 13# but WITHOUT ANY WARRANTY; without even the implied warranty of 14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15# GNU Lesser General Public License for more details. 16# 17# You should have received a copy of the GNU Lesser General Public License 18# along with this program. If not, see <https://www.gnu.org/licenses/>. 19 20 21import calendar 22import datetime 23from math import radians, cos, sin, asin, sqrt, acos, degrees 24 25from pymeeus.base import TOL, get_ordinal_suffix, iint 26from pymeeus.Angle import Angle 27 28 29""" 30.. module:: Epoch 31 :synopsis: Class to handle time 32 :license: GNU Lesser General Public License v3 (LGPLv3) 33 34.. moduleauthor:: Dagoberto Salazar 35""" 36 37 38DAY2SEC = 86400.0 39"""Number of seconds per day""" 40 41DAY2MIN = 1440.0 42"""Number of minutes per day""" 43 44DAY2HOURS = 24.0 45"""Number of hours per day""" 46 47LEAP_TABLE = { 48 1972.5: 1, 49 1973.0: 2, 50 1974.0: 3, 51 1975.0: 4, 52 1976.0: 5, 53 1977.0: 6, 54 1978.0: 7, 55 1979.0: 8, 56 1980.0: 9, 57 1981.5: 10, 58 1982.5: 11, 59 1983.5: 12, 60 1985.5: 13, 61 1988.0: 14, 62 1990.0: 15, 63 1991.0: 16, 64 1992.5: 17, 65 1993.5: 18, 66 1994.5: 19, 67 1996.0: 20, 68 1997.5: 21, 69 1999.0: 22, 70 2006.0: 23, 71 2009.0: 24, 72 2012.5: 25, 73 2015.5: 26, 74 2017.0: 27, 75} 76"""This table represents the point in time FROM WHERE the given number of leap 77seconds is valid. Given that leap seconds are (so far) always added at 78June 30th or December 31st, a leap second added in 1997/06/30 is represented 79here as '1997.5', while a leap second added in 2005/12/31 appears here as 80'2006.0'.""" 81 82 83class Epoch(object): 84 """ 85 Class Epoch deals with the tasks related to time handling. 86 87 The constructor takes either a single JDE value, another Epoch object, or a 88 series of values representing year, month, day, hours, minutes, seconds. 89 This series of values are by default supposed to be in **Terrestial Time** 90 (TT). 91 92 This is not necesarily the truth, though. For instance, the time of a 93 current observation is tipically in UTC time (civil time), not in TT, and 94 there is some offset between those two time references. 95 96 When a UTC time is provided, the parameter **utc=True** must be given. 97 Then, the input is converted to International Atomic Time (TAI) using an 98 internal table of leap seconds, and from there, it is converted to (and 99 stored as) Terrestrial Time (TT). 100 101 Given that leap seconds are added or subtracted in a rather irregular 102 basis, it is not possible to predict them in advance, and the internal leap 103 seconds table will become outdated at some point in time. To counter this, 104 you have two options: 105 106 - Download an updated version of this Pymeeus package. 107 - Use the argument **leap_seconds** in the constructor or :meth:`set` 108 method to provide the correct number of leap seconds (w.r.t. TAI) to be 109 applied. 110 111 .. note:: Providing the **leap_seconds** argument will automatically set 112 the argument **utc** to True. 113 114 For instance, if at some time in the future the TAI-UTC difference is 43 115 seconds, you should set **leap_seconds=43** if you don't have an updated 116 version of this class. 117 118 In order to know which is the most updated leap second value stored in this 119 class, you may use the :meth:`get_last_leap_second()` method. 120 121 .. note:: The current version of UTC was implemented in January 1st, 1972. 122 Therefore, for dates before the correction in **NOT** carried out, even 123 if the **utc** argument is set to True, and it is supposed that the 124 input data is already in TT scale. 125 126 .. note:: For conversions between TT and Universal Time (UT), please use 127 the method :meth:`tt2ut`. 128 129 .. note:: Internally, time values are stored as a Julian Ephemeris Day 130 (JDE), based on the uniform scale of Dynamical Time, or more 131 specifically, Terrestial Time (TT) (itself the redefinition of 132 Terrestrial Dynamical Time, TDT). 133 134 .. note:: The UTC-TT conversion is composed of three corrections: 135 136 a. TT-TAI, comprising 32.184 s, 137 b. TAI-UTC(1972), 10 s, and 138 c. UTC(1972)-UTC(target) 139 140 item c. is the corresponding amount of leap seconds to the target Epoch. 141 When you do, for instance, **leap_seconds=43**, you modify the c. part. 142 143 .. note:: Given that this class stores the epoch as JDE, if the JDE value 144 is in the order of millions of days then, for a computer with 15-digit 145 accuracy, the final time resolution is about 10 milliseconds. That is 146 considered enough for most applications of this class. 147 """ 148 149 def __init__(self, *args, **kwargs): 150 """Epoch constructor. 151 152 This constructor takes either a single JDE value, another Epoch object, 153 or a series of values representing year, month, day, hours, minutes, 154 seconds. This series of values are by default supposed to be in 155 **Terrestial Time** (TT). 156 157 It is also possible that the year, month, etc. arguments be provided in 158 a tuple or list. Moreover, it is also possible provide :class:`date` or 159 :class:`datetime` objects for initialization. 160 161 The **month** value can be provided as an integer (1 = January, 2 = 162 February, etc), or it can be provided with short (Jan, Feb,...) or long 163 (January, February,...) names. Also, hours, minutes, seconds can be 164 provided separately, or as decimals of the day value. 165 166 When a UTC time is provided, the parameter **utc=True** must be given. 167 Then, the input is converted to International Atomic Time (TAI) using 168 an internal table of leap seconds, and from there, it is converted to 169 (and stored as) Terrestrial Time (TT). If **utc** is not provided, it 170 is supposed that the input data is already in TT scale. 171 172 If a value is provided with the **leap_seconds** argument, then that 173 value will be used for the UTC->TAI conversion, and the internal leap 174 seconds table will be bypassed. 175 176 :param args: Either JDE, Epoch, date, datetime or year, month, day, 177 hours, minutes, seconds values, by themselves or inside a tuple or 178 list 179 :type args: int, float, :py:class:`Epoch`, tuple, list, date, 180 datetime 181 :param utc: Whether the provided epoch is a civil time (UTC) 182 :type utc: bool 183 :param leap_seconds: This is the value to be used in the UTC->TAI 184 conversion, instead of taking it from internal leap seconds table. 185 :type leap_seconds: int, float 186 187 :returns: Epoch object. 188 :rtype: :py:class:`Epoch` 189 :raises: ValueError if input values are in the wrong range. 190 :raises: TypeError if input values are of wrong type. 191 192 >>> e = Epoch(1987, 6, 19.5) 193 >>> print(e) 194 2446966.0 195 """ 196 197 # Initialize field 198 self._jde = 0.0 199 self.set(*args, **kwargs) # Use 'set()' method to handle the setup 200 201 def set(self, *args, **kwargs): 202 """Method used to set the value of this object. 203 204 This method takes either a single JDE value, or a series of values 205 representing year, month, day, hours, minutes, seconds. This series of 206 values are by default supposed to be in **Terrestial Time** (TT). 207 208 It is also possible to provide another Epoch object as input for the 209 :meth:`set` method, or the year, month, etc arguments can be provided 210 in a tuple or list. Moreover, it is also possible provide :class:`date` 211 or :class:`datetime` objects for initialization. 212 213 The **month** value can be provided as an integer (1 = January, 2 = 214 February, etc), or it can be provided as short (Jan, Feb, ...) or long 215 (January, February, ...) names. Also, hours, minutes, seconds can be 216 provided separately, or as decimals of the day value. 217 218 When a UTC time is provided, the parameter **utc=True** must be given. 219 Then, the input is converted to International Atomic Time (TAI) using 220 an internal table of leap seconds, and from there, it is converted to 221 (and stored as) Terrestrial Time (TT). If **utc** is not provided, it 222 is supposed that the input data is already in TT scale. 223 224 If a value is provided with the **leap_seconds** argument, then that 225 value will be used for the UTC->TAI conversion, and the internal leap 226 seconds table will be bypassed. 227 228 .. note:: The UTC to TT correction is only carried out for dates after 229 January 1st, 1972. 230 231 :param args: Either JDE, Epoch, date, datetime or year, month, day, 232 hours, minutes, seconds values, by themselves or inside a tuple or 233 list 234 :type args: int, float, :py:class:`Epoch`, tuple, list, date, 235 datetime 236 :param utc: Whether the provided epoch is a civil time (UTC) 237 :type utc: bool 238 :param leap_seconds: This is the value to be used in the UTC->TAI 239 conversion, instead of taking it from internal leap seconds table. 240 :type leap_seconds: int, float 241 242 :returns: None. 243 :rtype: None 244 :raises: ValueError if input values are in the wrong range. 245 :raises: TypeError if input values are of wrong type. 246 247 >>> e = Epoch() 248 >>> e.set(1987, 6, 19.5) 249 >>> print(e) 250 2446966.0 251 >>> e.set(1977, 'Apr', 26.4) 252 >>> print(e) 253 2443259.9 254 >>> e.set(1957, 'October', 4.81) 255 >>> print(e) 256 2436116.31 257 >>> e.set(333, 'Jan', 27, 12) 258 >>> print(e) 259 1842713.0 260 >>> e.set(1900, 'Jan', 1) 261 >>> print(e) 262 2415020.5 263 >>> e.set(-1001, 'august', 17.9) 264 >>> print(e) 265 1355671.4 266 >>> e.set(-4712, 1, 1.5) 267 >>> print(e) 268 0.0 269 >>> e.set((1600, 12, 31)) 270 >>> print(e) 271 2305812.5 272 >>> e.set([1988, 'JUN', 19, 12]) 273 >>> print(e) 274 2447332.0 275 >>> d = datetime.date(2000, 1, 1) 276 >>> e.set(d) 277 >>> print(e) 278 2451544.5 279 >>> e.set(837, 'Apr', 10, 7, 12) 280 >>> print(e) 281 2026871.8 282 >>> d = datetime.datetime(837, 4, 10, 7, 12, 0, 0) 283 >>> e.set(d) 284 >>> print(e) 285 2026871.8 286 """ 287 288 # Clean up the internal parameters 289 self._jde = 0.0 290 # If no arguments are given, return. Internal values are 0.0 291 if len(args) == 0: 292 return 293 # If we have only one argument, it can be a JDE or another Epoch object 294 elif len(args) == 1: 295 if isinstance(args[0], Epoch): 296 self._jde = args[0]._jde 297 return 298 elif isinstance(args[0], (int, float)): 299 self._jde = args[0] 300 return 301 elif isinstance(args[0], (tuple, list)): 302 year, month, day, hours, minutes, sec = \ 303 self._check_values(*args[0]) 304 elif isinstance(args[0], datetime.datetime): 305 d = args[0] 306 year, month, day, hours, minutes, sec = self._check_values( 307 d.year, 308 d.month, 309 d.day, 310 d.hour, 311 d.minute, 312 d.second + d.microsecond / 1e6, 313 ) 314 elif isinstance(args[0], datetime.date): 315 d = args[0] 316 year, month, day, hours, minutes, sec = self._check_values( 317 d.year, d.month, d.day 318 ) 319 else: 320 raise TypeError("Invalid input type") 321 elif len(args) == 2: 322 # Insuficient data to set the Epoch 323 raise ValueError("Invalid number of input values") 324 elif len(args) >= 3: # Year, month, day 325 year, month, day, hours, minutes, sec = self._check_values(*args) 326 day += hours / DAY2HOURS + minutes / DAY2MIN + sec / DAY2SEC 327 # Handle the 'leap_seconds' argument, if pressent 328 if "leap_seconds" in kwargs: 329 # Compute JDE 330 self._jde = self._compute_jde(year, month, day, utc2tt=False, 331 leap_seconds=kwargs["leap_seconds"]) 332 elif "utc" in kwargs: 333 self._jde = self._compute_jde(year, month, day, 334 utc2tt=kwargs["utc"]) 335 else: 336 self._jde = self._compute_jde(year, month, day, utc2tt=False) 337 338 def _compute_jde(self, y, m, d, utc2tt=True, leap_seconds=0.0): 339 """Method to compute the Julian Ephemeris Day (JDE). 340 341 .. note:: The UTC to TT correction is only carried out for dates after 342 January 1st, 1972. 343 344 :param y: Year 345 :type y: int 346 :param m: Month 347 :type m: int 348 :param d: Day 349 :type d: float 350 :param utc2tt: Whether correction UTC to TT is done automatically. 351 :type utc2tt: bool 352 :param leap_seconds: Number of leap seconds to apply 353 :type leap_seconds: float 354 355 :returns: Julian Ephemeris Day (JDE) 356 :rtype: float 357 """ 358 359 # The best approach here is first convert to JDE, and then adjust secs 360 if m <= 2: 361 y -= 1 362 m += 12 363 a = iint(y / 100.0) 364 b = 0.0 365 if not Epoch.is_julian(y, m, iint(d)): 366 b = 2.0 - a + iint(a / 4.0) 367 jde = (iint(365.25 * (y + 4716.0)) + 368 iint(30.6001 * (m + 1.0)) + d + b - 1524.5) 369 # If enabled, let's convert from UTC to TT, adding the needed seconds 370 deltasec = 0.0 371 # In this case, UTC to TT correction is applied automatically 372 if utc2tt: 373 if y >= 1972: 374 deltasec = 32.184 # Difference between TT and TAI 375 deltasec += 10.0 # Difference between UTC and TAI in 1972 376 deltasec += Epoch.leap_seconds(y, m) 377 else: # Correction is NOT automatic 378 if leap_seconds != 0.0: # We apply provided leap seconds 379 if y >= 1972: 380 deltasec = 32.184 # Difference between TT and TAI 381 deltasec += 10.0 # Difference between UTC-TAI in 1972 382 deltasec += leap_seconds 383 return jde + deltasec / DAY2SEC 384 385 def _check_values(self, *args): 386 """This method takes the input arguments to 'set()' method (year, 387 month, day, etc) and carries out some sanity checks on them. 388 389 It returns a tuple containing those values separately, assigning zeros 390 to those arguments which were not provided. 391 392 :param args: Year, month, day, hours, minutes, seconds values. 393 :type args: int, float 394 395 :returns: Tuple with year, month, day, hours, minutes, seconds values. 396 :rtype: tuple 397 :raises: ValueError if input values are in the wrong range, or too few 398 arguments given as input. 399 """ 400 401 # This list holds the maximum amount of days a given month can have 402 maxdays = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] 403 404 # Initialize some variables 405 year = -9999 406 month = -9999 407 day = -9999 408 hours = 0.0 409 minutes = 0.0 410 sec = 0.0 411 412 # Carry out some basic checks 413 if len(args) < 3: 414 raise ValueError("Invalid number of input values") 415 elif len(args) >= 3: # Year, month, day 416 year = args[0] 417 month = args[1] 418 day = args[2] 419 if len(args) >= 4: # Year, month, day, hour 420 hours = args[3] 421 if len(args) >= 5: # Year, month, day, hour, minutes 422 minutes = args[4] 423 if len(args) >= 6: # Year, month, day, hour, minutes, seconds 424 sec = args[5] 425 if year < -4712: # No negative JDE will be allowed 426 raise ValueError("Invalid value for the input year") 427 if day < 1 or day > 31: 428 raise ValueError("Invalid value for the input day") 429 if hours < 0 or hours > 23: 430 raise ValueError("Invalid value for the input hours") 431 if minutes < 0 or minutes > 59: 432 raise ValueError("Invalid value for the input minutes") 433 if sec < 0 or sec > 59: 434 raise ValueError("Invalid value for the input seconds") 435 436 # Test the days according to the month 437 month = Epoch.get_month(month) 438 limit_day = maxdays[month - 1] 439 # We need extra tests if month is '2' (February) 440 if month == 2: 441 if Epoch.is_leap(year): 442 limit_day = 29 443 if day > limit_day: 444 raise ValueError("Invalid value for the input day") 445 446 # We are ready to return the parameters 447 return year, month, day, hours, minutes, sec 448 449 @staticmethod 450 def check_input_date(*args, **kwargs): 451 """Method to check that the input is a proper date. 452 453 This method returns an Epoch object, and the **leap_seconds** argument 454 then controls the way the UTC->TT conversion is handled for that new 455 object. If **leap_seconds** argument is set to a value different than 456 zero, then that value will be used for the UTC->TAI conversion, and the 457 internal leap seconds table will be bypassed. On the other hand, if it 458 is set to zero, then the UTC to TT correction is disabled, and it is 459 supposed that the input data is already in TT scale. 460 461 :param args: Either Epoch, date, datetime or year, month, day values, 462 by themselves or inside a tuple or list 463 :type args: int, float, :py:class:`Epoch`, datetime, date, tuple, 464 list 465 :param leap_seconds: If different from zero, this is the value to be 466 used in the UTC->TAI conversion. If equals to zero, conversion is 467 disabled. If not given, UTC to TT conversion is carried out 468 (default). 469 :type leap_seconds: int, float 470 471 :returns: Epoch object corresponding to the input date 472 :rtype: :py:class:`Epoch` 473 :raises: ValueError if input values are in the wrong range. 474 :raises: TypeError if input values are of wrong type. 475 """ 476 477 t = Epoch() 478 if len(args) == 0: 479 raise ValueError("Invalid input: No date given") 480 # If we have only one argument, it can be an Epoch, a date, a datetime 481 # or a tuple/list 482 elif len(args) == 1: 483 if isinstance(args[0], Epoch): 484 t = args[0] 485 elif isinstance(args[0], (tuple, list)): 486 if len(args[0]) >= 3: 487 t = Epoch(args[0][0], args[0][1], args[0][2], **kwargs) 488 else: 489 raise ValueError("Invalid input") 490 elif isinstance(args[0], datetime.datetime) or isinstance( 491 args[0], datetime.date 492 ): 493 t = Epoch(args[0].year, args[0].month, args[0].day, **kwargs) 494 else: 495 raise TypeError("Invalid input type") 496 elif len(args) == 2: 497 raise ValueError("Invalid input: Date given is not valid") 498 elif len(args) >= 3: 499 # We will rely on Epoch capacity to handle improper input 500 t = Epoch(args[0], args[1], args[2], **kwargs) 501 return t 502 503 @staticmethod 504 def is_julian(year, month, day): 505 """This method returns True if given date is in the Julian calendar. 506 507 :param year: Year 508 :type y: int 509 :param month: Month 510 :type m: int 511 :param day: Day 512 :type day: int 513 514 :returns: Whether the provided date belongs to Julian calendar or not. 515 :rtype: bool 516 517 >>> Epoch.is_julian(1997, 5, 27.1) 518 False 519 >>> Epoch.is_julian(1397, 7, 7.0) 520 True 521 """ 522 523 if ( 524 (year < 1582) 525 or (year == 1582 and month < 10) 526 or (year == 1582 and month == 10 and day < 5.0) 527 ): 528 return True 529 else: 530 return False 531 532 def julian(self): 533 """This method returns True if this Epoch object holds a date in the 534 Julian calendar. 535 536 :returns: Whether this Epoch object holds a date belonging to Julian 537 calendar or not. 538 :rtype: bool 539 540 >>> e = Epoch(1997, 5, 27.1) 541 >>> e.julian() 542 False 543 >>> e = Epoch(1397, 7, 7.0) 544 >>> e.julian() 545 True 546 """ 547 548 y, m, d = self.get_date() 549 return Epoch.is_julian(y, m, d) 550 551 @staticmethod 552 def get_month(month, as_string=False): 553 """Method to get the month as a integer in the [1, 12] range, or as a 554 full name. 555 556 :param month: Month, in numeric, short name or long name format 557 :type month: int, float, str 558 :param as_string: Whether the output will be numeric, or a long name. 559 :type as_string: bool 560 561 :returns: Month as integer in the [1, 12] range, or as a long name. 562 :rtype: int, str 563 :raises: ValueError if input month value is invalid. 564 565 >>> Epoch.get_month(4.0) 566 4 567 >>> Epoch.get_month('Oct') 568 10 569 >>> Epoch.get_month('FEB') 570 2 571 >>> Epoch.get_month('August') 572 8 573 >>> Epoch.get_month('august') 574 8 575 >>> Epoch.get_month('NOVEMBER') 576 11 577 >>> Epoch.get_month(9.0, as_string=True) 578 'September' 579 >>> Epoch.get_month('Feb', as_string=True) 580 'February' 581 >>> Epoch.get_month('March', as_string=True) 582 'March' 583 """ 584 585 months_mmm = [ 586 "Jan", 587 "Feb", 588 "Mar", 589 "Apr", 590 "May", 591 "Jun", 592 "Jul", 593 "Aug", 594 "Sep", 595 "Oct", 596 "Nov", 597 "Dec", 598 ] 599 600 months_full = [ 601 "January", 602 "February", 603 "March", 604 "April", 605 "May", 606 "June", 607 "July", 608 "August", 609 "September", 610 "October", 611 "November", 612 "December", 613 ] 614 615 if isinstance(month, (int, float)): 616 month = int(month) # Truncate if it has decimals 617 if month >= 1 and month <= 12: 618 if not as_string: 619 return month 620 else: 621 return months_full[month - 1] 622 else: 623 raise ValueError("Invalid value for the input month") 624 elif isinstance(month, str): 625 month = month.strip().capitalize() 626 if len(month) == 3: 627 if month in months_mmm: 628 if not as_string: 629 return months_mmm.index(month) + 1 630 else: 631 return months_full[months_mmm.index(month)] 632 else: 633 raise ValueError("Invalid value for the input month") 634 else: 635 if month in months_full: 636 if not as_string: 637 return months_full.index(month) + 1 638 else: 639 return month 640 else: 641 raise ValueError("Invalid value for the input month") 642 643 @staticmethod 644 def is_leap(year): 645 """Method to check if a given year is a leap year. 646 647 :param year: Year to be checked. 648 :type year: int, float 649 650 :returns: Whether or not year is a leap year. 651 :rtype: bool 652 :raises: ValueError if input year value is invalid. 653 654 >>> Epoch.is_leap(2003) 655 False 656 >>> Epoch.is_leap(2012) 657 True 658 >>> Epoch.is_leap(1900) 659 False 660 >>> Epoch.is_leap(-1000) 661 True 662 >>> Epoch.is_leap(1000) 663 True 664 """ 665 666 if isinstance(year, (int, float)): 667 # Mind the difference between Julian and Gregorian calendars 668 if year >= 1582: 669 year = iint(year) 670 return calendar.isleap(year) 671 else: 672 return (abs(year) % 4) == 0 673 else: 674 raise ValueError("Invalid value for the input year") 675 676 def leap(self): 677 """This method checks if the current Epoch object holds a leap year. 678 679 :returns: Whether or the year in this Epoch object is a leap year. 680 :rtype: bool 681 682 >>> e = Epoch(2003, 1, 1) 683 >>> e.leap() 684 False 685 >>> e = Epoch(2012, 1, 1) 686 >>> e.leap() 687 True 688 >>> e = Epoch(1900, 1, 1) 689 >>> e.leap() 690 False 691 >>> e = Epoch(-1000, 1, 1) 692 >>> e.leap() 693 True 694 >>> e = Epoch(1000, 1, 1) 695 >>> e.leap() 696 True 697 """ 698 699 y, m, d = self.get_date() 700 return Epoch.is_leap(y) 701 702 @staticmethod 703 def get_doy(yyyy, mm, dd): 704 """This method returns the Day Of Year (DOY) for the given date. 705 706 :param yyyy: Year, in four digits format 707 :type yyyy: int, float 708 :param mm: Month, in numeric format (1 = January, 2 = February, etc) 709 :type mm: int, float 710 :param dd: Day, in numeric format 711 :type dd: int, float 712 713 :returns: Day Of Year (DOY). 714 :rtype: float 715 :raises: ValueError if input values correspond to a wrong date. 716 717 >>> Epoch.get_doy(1999, 1, 29) 718 29.0 719 >>> Epoch.get_doy(1978, 11, 14) 720 318.0 721 >>> Epoch.get_doy(2017, 12, 31.7) 722 365.7 723 >>> Epoch.get_doy(2012, 3, 3.1) 724 63.1 725 >>> Epoch.get_doy(-400, 2, 29.9) 726 60.9 727 """ 728 729 # Let's carry out first some basic checks 730 if dd < 1 or dd >= 32 or mm < 1 or mm > 12: 731 raise ValueError("Invalid input data") 732 day = int(dd) 733 frac = dd % 1 734 if yyyy >= 1: # datetime's minimum year is 1 735 try: 736 d = datetime.date(yyyy, mm, day) 737 except ValueError: 738 raise ValueError("Invalid input date") 739 doy = d.timetuple().tm_yday 740 else: 741 k = 2 if Epoch.is_leap(yyyy) else 1 742 doy = (iint((275.0 * mm) / 9.0) - 743 k * iint((mm + 9.0) / 12.0) + day - 30.0) 744 return float(doy + frac) 745 746 @staticmethod 747 def doy2date(year, doy): 748 """This method takes a year and a Day Of Year values, and returns the 749 corresponding date. 750 751 :param year: Year, in four digits format 752 :type year: int, float 753 :param doy: Day of Year number 754 :type doy: int, float 755 756 :returns: Year, month, day. 757 :rtype: tuple 758 :raises: ValueError if either input year or doy values are invalid. 759 760 >>> t = Epoch.doy2date(1999, 29) 761 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 762 1999/1/29.0 763 >>> t = Epoch.doy2date(2017, 365.7) 764 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 765 2017/12/31.7 766 >>> t = Epoch.doy2date(2012, 63.1) 767 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 768 2012/3/3.1 769 >>> t = Epoch.doy2date(-1004, 60) 770 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 771 -1004/2/29.0 772 >>> t = Epoch.doy2date(0, 60) 773 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 774 0/2/29.0 775 >>> t = Epoch.doy2date(1, 60) 776 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 777 1/3/1.0 778 >>> t = Epoch.doy2date(-1, 60) 779 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 780 -1/3/1.0 781 >>> t = Epoch.doy2date(-2, 60) 782 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 783 -2/3/1.0 784 >>> t = Epoch.doy2date(-3, 60) 785 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 786 -3/3/1.0 787 >>> t = Epoch.doy2date(-4, 60) 788 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 789 -4/2/29.0 790 >>> t = Epoch.doy2date(-5, 60) 791 >>> print("{}/{}/{}".format(t[0], t[1], round(t[2], 1))) 792 -5/3/1.0 793 """ 794 795 if isinstance(year, (int, float)) and isinstance(doy, (int, float)): 796 frac = float(doy % 1) 797 doy = int(doy) 798 if year >= 1: # datetime's minimum year is 1 799 ref = datetime.date(year, 1, 1) 800 mydate = datetime.date.fromordinal(ref.toordinal() + doy - 1) 801 return year, mydate.month, mydate.day + frac 802 else: 803 # The algorithm provided by Meeus doesn't work for years below 804 # +1. This little hack solves that problem (the 'if' result is 805 # inverted here). 806 k = 1 if Epoch.is_leap(year) else 2 807 if doy < 32: 808 m = 1 809 else: 810 m = iint((9.0 * (k + doy)) / 275.0 + 0.98) 811 d = (doy - iint((275.0 * m) / 9.0) + 812 k * iint((m + 9.0) / 12.0) + 30) 813 return year, int(m), d + frac 814 else: 815 raise ValueError("Invalid input values") 816 817 @staticmethod 818 def leap_seconds(year, month): 819 """Returns the leap seconds accumulated for the given year and month. 820 821 :param year: Year 822 :type year: int 823 :param month: Month, in numeric format ([1:12] range) 824 :type month: int 825 826 :returns: Leap seconds accumulated for given year and month. 827 :rtype: int 828 829 >>> Epoch.leap_seconds(1972, 4) 830 0 831 >>> Epoch.leap_seconds(1972, 6) 832 0 833 >>> Epoch.leap_seconds(1972, 7) 834 1 835 >>> Epoch.leap_seconds(1983, 6) 836 11 837 >>> Epoch.leap_seconds(1983, 7) 838 12 839 >>> Epoch.leap_seconds(1985, 8) 840 13 841 >>> Epoch.leap_seconds(2016, 11) 842 26 843 >>> Epoch.leap_seconds(2017, 1) 844 27 845 >>> Epoch.leap_seconds(2018, 7) 846 27 847 """ 848 849 list_years = sorted(LEAP_TABLE.keys()) 850 # First test the extremes of the table 851 if (year + month / 12.0) <= list_years[0]: 852 return 0 853 if (year + month / 12.0) >= list_years[-1]: 854 return LEAP_TABLE[list_years[-1]] 855 lyear = (year + 0.25) if month <= 6 else (year + 0.75) 856 idx = 0 857 while lyear > list_years[idx]: 858 idx += 1 859 return LEAP_TABLE[list_years[idx - 1]] 860 861 @staticmethod 862 def get_last_leap_second(): 863 """Method to get the date and value of the last leap second added to 864 the table 865 866 :returns: Tuple with year, month, day, leap second value. 867 :rtype: tuple 868 """ 869 870 list_years = sorted(LEAP_TABLE.keys()) 871 lyear = list_years[-1] 872 lseconds = LEAP_TABLE[lyear] 873 year = iint(lyear) 874 # So far, leap seconds are added either on June 30th or December 31th 875 if lyear % 1 == 0.0: 876 year -= 1 877 month = 12 878 day = 31.0 879 else: 880 month = 6 881 day = 30.0 882 return year, month, day, lseconds 883 884 @staticmethod 885 def utc2local(): 886 """Method to return the difference between UTC and local time. 887 888 By default, dates in this Epoch class are handled in Coordinated 889 Universal Time (UTC). This method provides you the seconds that you 890 have to add or subtract to UTC time to convert to your local time. 891 892 Please bear in mind that, in order for this method to work, you 893 operative system must be correctly configured, with the right time and 894 corresponding time zone. 895 896 :returns: Difference in seconds between local and UTC time. 897 :rtype: float 898 """ 899 900 localhour = datetime.datetime.now().hour 901 utchour = datetime.datetime.utcnow().hour 902 localminute = datetime.datetime.now().minute 903 utcminute = datetime.datetime.utcnow().minute 904 return ((localhour - utchour) * 3600.0 + 905 (localminute - utcminute) * 60.0) 906 907 @staticmethod 908 def easter(year): 909 """Method to return the Easter day for given year. 910 911 .. note:: This method is valid for both Gregorian and Julian years. 912 913 :param year: Year 914 :type year: int 915 916 :returns: Easter month and day, as a tuple 917 :rtype: tuple 918 :raises: TypeError if input values are of wrong type. 919 920 >>> Epoch.easter(1991) 921 (3, 31) 922 >>> Epoch.easter(1818) 923 (3, 22) 924 >>> Epoch.easter(1943) 925 (4, 25) 926 >>> Epoch.easter(2000) 927 (4, 23) 928 >>> Epoch.easter(1954) 929 (4, 18) 930 >>> Epoch.easter(179) 931 (4, 12) 932 >>> Epoch.easter(1243) 933 (4, 12) 934 """ 935 936 # This algorithm is describes in pages 67-69 of Meeus book 937 if not isinstance(year, (int, float)): 938 raise TypeError("Invalid input type") 939 year = int(year) 940 if year >= 1583: 941 # In this case, we are using the Gregorian calendar 942 a = year % 19 943 b = iint(year / 100.0) 944 c = year % 100 945 d = iint(b / 4.0) 946 e = b % 4 947 f = iint((b + 8.0) / 25.0) 948 g = iint((b - f + 1.0) / 3.0) 949 h = (19 * a + b - d - g + 15) % 30 950 i = iint(c / 4.0) 951 k = c % 4 952 ll = (32 + 2 * (e + i) - h - k) % 7 953 m = iint((a + 11 * h + 22 * ll) / 451.0) 954 n = iint((h + ll - 7 * m + 114) / 31.0) 955 p = (h + ll - 7 * m + 114) % 31 956 return (n, p + 1) 957 else: 958 # The Julian calendar is used here 959 a = year % 4 960 b = year % 7 961 c = year % 19 962 d = (19 * c + 15) % 30 963 e = (2 * a + 4 * b - d + 34) % 7 964 f = iint((d + e + 114) / 31.0) 965 g = (d + e + 114) % 31 966 return (f, g + 1) 967 968 @staticmethod 969 def jewish_pesach(year): 970 """Method to return the Jewish Easter (Pesach) day for given year. 971 972 .. note:: This method is valid for both Gregorian and Julian years. 973 974 :param year: Year 975 :type year: int 976 977 :returns: Jewish Easter (Pesach) month and day, as a tuple 978 :rtype: tuple 979 :raises: TypeError if input values are of wrong type. 980 981 >>> Epoch.jewish_pesach(1990) 982 (4, 10) 983 """ 984 985 # This algorithm is described in pages 71-73 of Meeus book 986 if not isinstance(year, (int, float)): 987 raise TypeError("Invalid input type") 988 year = iint(year) 989 c = iint(year / 100.0) 990 s = 0 if year < 1583 else iint((3.0 * c - 5.0) / 4.0) 991 a = (12 * (year + 1)) % 19 992 b = year % 4 993 q = (-1.904412361576 + 1.554241796621 * a + 994 0.25 * b - 0.003177794022 * year + s) 995 j = (iint(q) + 3 * year + 5 * b + 2 + s) % 7 996 r = q - iint(q) 997 if j == 2 or j == 4 or j == 6: 998 d = iint(q) + 23 999 elif j == 1 and a > 6 and r > 0.632870370: 1000 d = iint(q) + 24 1001 elif j == 0 and a > 11 and r > 0.897723765: 1002 d = iint(q) + 23 1003 else: 1004 d = iint(q) + 22 1005 if d > 31: 1006 return (4, d - 31) 1007 else: 1008 return (3, d) 1009 1010 @staticmethod 1011 def moslem2gregorian(year, month, day): 1012 """Method to convert a date in the Moslen calendar to the Gregorian 1013 (or Julian) calendar. 1014 1015 .. note:: This method is valid for both Gregorian and Julian years. 1016 1017 :param year: Year 1018 :type year: int 1019 :param month: Month 1020 :type month: int 1021 :param day: Day 1022 :type day: int 1023 1024 :returns: Date in Gregorian (Julian) calendar: year, month and day, as 1025 a tuple 1026 :rtype: tuple 1027 :raises: TypeError if input values are of wrong type. 1028 1029 >>> Epoch.moslem2gregorian(1421, 1, 1) 1030 (2000, 4, 6) 1031 """ 1032 1033 # First, check that input types are correct 1034 if ( 1035 not isinstance(year, (int, float)) 1036 or not isinstance(month, (int, float)) 1037 or not isinstance(day, (int, float)) 1038 ): 1039 raise TypeError("Invalid input type") 1040 if day < 1 or day > 30 or month < 1 or month > 12 or year < 1: 1041 raise ValueError("Invalid input data") 1042 # This algorithm is described in pages 73-75 of Meeus book 1043 # Note: Ramadan is month Nr. 9 1044 h = iint(year) 1045 m = iint(month) 1046 d = iint(day) 1047 n = d + iint(29.5001 * (m - 1) + 0.99) 1048 q = iint(h / 30.0) 1049 r = h % 30 1050 a = iint((11.0 * r + 3.0) / 30.0) 1051 w = 404 * q + 354 * r + 208 + a 1052 q1 = iint(w / 1461.0) 1053 q2 = w % 1461 1054 g = 621 + 4 * iint(7.0 * q + q1) 1055 k = iint(q2 / 365.2422) 1056 e = iint(365.2422 * k) 1057 j = q2 - e + n - 1 1058 x = g + k 1059 if j > 366 and x % 4 == 0: 1060 j -= 366 1061 x += 1 1062 elif j > 365 and x % 4 > 0: 1063 j -= 365 1064 x += 1 1065 1066 # Check if date is in Gregorian calendar. '277' is DOY of October 4th 1067 if (x > 1583) or (x == 1582 and j > 277): 1068 jd = iint(365.25 * (x - 1.0)) + 1721423 + j 1069 alpha = iint((jd - 1867216.25) / 36524.25) 1070 beta = jd if jd < 2299161 else (jd + 1 + alpha - iint(alpha / 4.0)) 1071 b = beta + 1524 1072 c = iint((b - 122.1) / 365.25) 1073 d = iint(365.25 * c) 1074 e = iint((b - d) / 30.6001) 1075 day = b - d - iint(30.6001 * e) 1076 month = (e - 1) if e < 14 else (e - 13) 1077 year = (c - 4716) if month > 2 else (c - 4715) 1078 return year, month, day 1079 else: 1080 # It is a Julian date. We have year and DOY 1081 return Epoch.doy2date(x, j) 1082 1083 @staticmethod 1084 def gregorian2moslem(year, month, day): 1085 """Method to convert a date in the Gregorian (or Julian) calendar to 1086 the Moslen calendar. 1087 1088 :param year: Year 1089 :type year: int 1090 :param month: Month 1091 :type month: int 1092 :param day: Day 1093 :type day: int 1094 1095 :returns: Date in Moslem calendar: year, month and day, as a tuple 1096 :rtype: tuple 1097 :raises: TypeError if input values are of wrong type. 1098 1099 >>> Epoch.gregorian2moslem(1991, 8, 13) 1100 (1412, 2, 2) 1101 """ 1102 1103 # First, check that input types are correct 1104 if ( 1105 not isinstance(year, (int, float)) 1106 or not isinstance(month, (int, float)) 1107 or not isinstance(day, (int, float)) 1108 ): 1109 raise TypeError("Invalid input type") 1110 if day < 1 or day > 31 or month < 1 or month > 12 or year < -4712: 1111 raise ValueError("Invalid input data") 1112 # This algorithm is described in pages 75-76 of Meeus book 1113 x = iint(year) 1114 m = iint(month) 1115 d = iint(day) 1116 if m < 3: 1117 x -= 1 1118 m += 12 1119 alpha = iint(x / 100.0) 1120 beta = 2 - alpha + iint(alpha / 4.0) 1121 b = iint(365.25 * x) + iint(30.6001 * (m + 1.0)) + d + 1722519 + beta 1122 c = iint((b - 122.1) / 365.25) 1123 d = iint(365.25 * c) 1124 e = iint((b - d) / 30.6001) 1125 d = b - d - iint(30.6001 * e) 1126 m = (e - 1) if e < 14 else (e - 13) 1127 x = (c - 4716) if month > 2 else (c - 4715) 1128 w = 1 if x % 4 == 0 else 2 1129 n = iint((275.0 * m) / 9.0) - w * iint((m + 9.0) / 12.0) + d - 30 1130 a = x - 623 1131 b = iint(a / 4.0) 1132 c = a % 4 1133 c1 = 365.2501 * c 1134 c2 = iint(c1) 1135 if c1 - c2 > 0.5: 1136 c2 += 1 1137 dp = 1461 * b + 170 + c2 1138 q = iint(dp / 10631.0) 1139 r = dp % 10631 1140 j = iint(r / 354.0) 1141 k = r % 354 1142 o = iint((11.0 * j + 14.0) / 30.0) 1143 h = 30 * q + j + 1 1144 jj = k - o + n - 1 1145 # jj is the number of the day in the moslem year h. If jj > 354 we need 1146 # to know if h is a leap year 1147 if jj > 354: 1148 cl = h % 30 1149 dl = (11 * cl + 3) % 30 1150 if dl < 19: 1151 jj -= 354 1152 h += 1 1153 else: 1154 jj -= 355 1155 h += 1 1156 if jj == 0: 1157 jj = 355 1158 h -= 1 1159 # Now, let's convert DOY jj to month and day 1160 if jj == 355: 1161 m = 12 1162 d = 30 1163 else: 1164 s = iint((jj - 1.0) / 29.5) 1165 m = 1 + s 1166 d = iint(jj - 29.5 * s) 1167 return h, m, d 1168 1169 def __str__(self): 1170 """Method used when trying to print the object. 1171 1172 :returns: Internal JDE value as a string. 1173 :rtype: string 1174 1175 >>> e = Epoch(1987, 6, 19.5) 1176 >>> print(e) 1177 2446966.0 1178 """ 1179 1180 return str(self._jde) 1181 1182 def __repr__(self): 1183 """Method providing the 'official' string representation of the object. 1184 It provides a valid expression that could be used to recreate the 1185 object. 1186 1187 :returns: As string with a valid expression to recreate the object 1188 :rtype: string 1189 1190 >>> e = Epoch(1987, 6, 19.5) 1191 >>> repr(e) 1192 'Epoch(2446966.0)' 1193 """ 1194 1195 return "{}({})".format(self.__class__.__name__, self._jde) 1196 1197 def get_date(self, **kwargs): 1198 """This method converts the internal JDE value back to a date. 1199 1200 Use **utc=True** to enable the TT to UTC conversion mechanism, or 1201 provide a non zero value to **leap_seconds** to apply a specific leap 1202 seconds value. 1203 1204 :param utc: Whether the TT to UTC conversion mechanism will be enabled 1205 :type utc: bool 1206 :param leap_seconds: Optional value for leap seconds. 1207 :type leap_seconds: int, float 1208 1209 :returns: Year, month, day in a tuple 1210 :rtype: tuple 1211 1212 >>> e = Epoch(2436116.31) 1213 >>> y, m, d = e.get_date() 1214 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1215 1957/10/4.81 1216 >>> e = Epoch(1988, 1, 27) 1217 >>> y, m, d = e.get_date() 1218 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1219 1988/1/27.0 1220 >>> e = Epoch(1842713.0) 1221 >>> y, m, d = e.get_date() 1222 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1223 333/1/27.5 1224 >>> e = Epoch(1507900.13) 1225 >>> y, m, d = e.get_date() 1226 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1227 -584/5/28.63 1228 """ 1229 1230 jd = self._jde + 0.5 1231 z = iint(jd) 1232 f = jd % 1 1233 if z < 2299161: 1234 a = z 1235 else: 1236 alpha = iint((z - 1867216.25) / 36524.25) 1237 a = z + 1 + alpha - iint(alpha / 4.0) 1238 b = a + 1524 1239 c = iint((b - 122.1) / 365.25) 1240 d = iint(365.25 * c) 1241 e = iint((b - d) / 30.6001) 1242 day = b - d - iint(30.6001 * e) + f 1243 if e < 14: 1244 month = e - 1 1245 elif e == 14 or e == 15: 1246 month = e - 13 1247 if month > 2: 1248 year = c - 4716 1249 elif month == 1 or month == 2: 1250 year = c - 4715 1251 year = int(year) 1252 month = int(month) 1253 1254 tt2utc = False 1255 if "utc" in kwargs: 1256 tt2utc = kwargs["utc"] 1257 if "leap_seconds" in kwargs: 1258 tt2utc = False 1259 leap_seconds = kwargs["leap_seconds"] 1260 else: 1261 leap_seconds = 0.0 1262 # If enabled, let's convert from TT to UTC, subtracting needed seconds 1263 deltasec = 0.0 1264 # In this case, TT to UTC correction is applied automatically, but only 1265 # for dates after July 1st, 1972 1266 if tt2utc: 1267 if year > 1972 or (year == 1972 and month >= 7): 1268 deltasec = 32.184 # Difference between TT and TAI 1269 deltasec += 10.0 # Difference between UTC and TAI in 1972 1270 deltasec += Epoch.leap_seconds(year, month) 1271 else: # Correction is NOT automatic 1272 if leap_seconds != 0.0: # We apply provided leap seconds 1273 if year > 1972 or (year == 1972 and month >= 7): 1274 deltasec = 32.184 # Difference between TT and TAI 1275 deltasec += 10.0 # Difference between UTC-TAI in 1972 1276 deltasec += leap_seconds 1277 1278 if deltasec != 0.0: 1279 doy = Epoch.get_doy(year, month, day) 1280 doy -= deltasec / DAY2SEC 1281 # Check that we didn't change year 1282 if doy < 1.0: 1283 year -= 1 1284 doy = 366.0 + doy if Epoch.is_leap(year) else 365.0 + doy 1285 year, month, day = Epoch.doy2date(year, doy) 1286 return year, month, day 1287 1288 def get_full_date(self, **kwargs): 1289 """This method converts the internal JDE value back to a full date. 1290 1291 Use **utc=True** to enable the TT to UTC conversion mechanism, or 1292 provide a non zero value to **leap_seconds** to apply a specific leap 1293 seconds value. 1294 1295 :param utc: Whether the TT to UTC conversion mechanism will be enabled 1296 :type utc: bool 1297 :param leap_seconds: Optional value for leap seconds. 1298 :type leap_seconds: int, float 1299 1300 :returns: Year, month, day, hours, minutes, seconds in a tuple 1301 :rtype: tuple 1302 1303 >>> e = Epoch(2436116.31) 1304 >>> y, m, d, h, mi, s = e.get_full_date() 1305 >>> print("{}/{}/{} {}:{}:{}".format(y, m, d, h, mi, round(s, 1))) 1306 1957/10/4 19:26:24.0 1307 >>> e = Epoch(1988, 1, 27) 1308 >>> y, m, d, h, mi, s = e.get_full_date() 1309 >>> print("{}/{}/{} {}:{}:{}".format(y, m, d, h, mi, round(s, 1))) 1310 1988/1/27 0:0:0.0 1311 >>> e = Epoch(1842713.0) 1312 >>> y, m, d, h, mi, s = e.get_full_date() 1313 >>> print("{}/{}/{} {}:{}:{}".format(y, m, d, h, mi, round(s, 1))) 1314 333/1/27 12:0:0.0 1315 >>> e = Epoch(1507900.13) 1316 >>> y, m, d, h, mi, s = e.get_full_date() 1317 >>> print("{}/{}/{} {}:{}:{}".format(y, m, d, h, mi, round(s, 1))) 1318 -584/5/28 15:7:12.0 1319 """ 1320 1321 y, m, d = self.get_date(**kwargs) 1322 r = d % 1 1323 d = int(d) 1324 h = int(r * 24.0) 1325 r = r * 24 - h 1326 mi = int(r * 60.0) 1327 s = 60.0 * (r * 60.0 - mi) 1328 return y, m, d, h, mi, s 1329 1330 @staticmethod 1331 def tt2ut(year, month): 1332 """This method provides an approximation of the difference, in seconds, 1333 between Terrestrial Time and Universal Time, denoted **DeltaT**, where: 1334 DeltaT = TT - UT. 1335 1336 Here we depart from Meeus book and use the polynomial expressions from: 1337 1338 https://eclipse.gsfc.nasa.gov/LEcat5/deltatpoly.html 1339 1340 Which are regarded as more elaborate and precise than Meeus'. 1341 1342 Please note that, by definition, the UTC time used internally in this 1343 Epoch class by default is kept within 0.9 seconds from UT. Therefore, 1344 UTC is in itself a quite good approximation to UT, arguably better than 1345 some of the results provided by this method. 1346 1347 :param year: Year we want to compute DeltaT for. 1348 :type year: int, float 1349 :param month: Month we want to compute DeltaT for. 1350 :type month: int, float 1351 1352 :returns: DeltaT, in seconds 1353 :rtype: float 1354 1355 >>> round(Epoch.tt2ut(1642, 1), 1) 1356 62.1 1357 >>> round(Epoch.tt2ut(1680, 1), 1) 1358 15.3 1359 >>> round(Epoch.tt2ut(1700, 1), 1) 1360 8.8 1361 >>> round(Epoch.tt2ut(1726, 1), 1) 1362 10.9 1363 >>> round(Epoch.tt2ut(1750, 1), 1) 1364 13.4 1365 >>> round(Epoch.tt2ut(1774, 1), 1) 1366 16.7 1367 >>> round(Epoch.tt2ut(1800, 1), 1) 1368 13.7 1369 >>> round(Epoch.tt2ut(1820, 1), 1) 1370 11.9 1371 >>> round(Epoch.tt2ut(1890, 1), 1) 1372 -6.1 1373 >>> round(Epoch.tt2ut(1928, 2), 1) 1374 24.2 1375 >>> round(Epoch.tt2ut(1977, 2), 1) 1376 47.7 1377 >>> round(Epoch.tt2ut(1998, 1), 1) 1378 63.0 1379 >>> round(Epoch.tt2ut(2015, 7), 1) 1380 69.3 1381 """ 1382 1383 y = year + (month - 0.5) / 12.0 1384 if year < -500: 1385 u = (year - 1820.0) / 100.0 1386 dt = -20.0 + 32.0 * u * u 1387 elif year >= -500 and year < 500: 1388 u = y / 100.0 1389 dt = 10583.6 + u * ( 1390 -1014.41 1391 + u 1392 * ( 1393 33.78311 1394 + u 1395 * ( 1396 -5.952053 1397 + (u * (-0.1798452 + 1398 u * (0.022174192 + 0.0090316521 * u))) 1399 ) 1400 ) 1401 ) 1402 elif year >= 500 and year < 1600: 1403 dt = 1574.2 + u * ( 1404 -556.01 1405 + u 1406 * ( 1407 71.23472 1408 + u 1409 * ( 1410 0.319781 1411 + (u * (-0.8503463 + 1412 u * (-0.005050998 + 0.0083572073 * u))) 1413 ) 1414 ) 1415 ) 1416 elif year >= 1600 and year < 1700: 1417 t = y - 1600.0 1418 dt = 120.0 + t * (-0.9808 + t * (-0.01532 + t / 7129.0)) 1419 elif year >= 1700 and year < 1800: 1420 t = y - 1700.0 1421 dt = 8.83 + t * ( 1422 0.1603 + t * (-0.0059285 + t * (0.00013336 - t / 1174000.0)) 1423 ) 1424 elif year >= 1800 and year < 1860: 1425 t = y - 1800.0 1426 dt = 13.72 + t * ( 1427 -0.332447 1428 + t 1429 * ( 1430 0.0068612 1431 + t 1432 * ( 1433 0.0041116 1434 + t 1435 * ( 1436 -0.00037436 1437 + t 1438 * (0.0000121272 + t * (-0.0000001699 + 1439 0.000000000875 * t)) 1440 ) 1441 ) 1442 ) 1443 ) 1444 elif year >= 1860 and year < 1900: 1445 t = y - 1860.0 1446 dt = 7.62 + t * ( 1447 0.5737 1448 + t 1449 * (-0.251754 + t * (0.01680668 + 1450 t * (-0.0004473624 + t / 233174.0))) 1451 ) 1452 elif year >= 1900 and year < 1920: 1453 t = y - 1900.0 1454 dt = -2.79 + t * ( 1455 1.494119 + t * (-0.0598939 + t * (0.0061966 - 0.000197 * t)) 1456 ) 1457 elif year >= 1920 and year < 1941: 1458 t = y - 1920.0 1459 dt = 21.20 + t * (0.84493 + t * (-0.076100 + 0.0020936 * t)) 1460 elif year >= 1941 and year < 1961: 1461 t = y - 1950.0 1462 dt = 29.07 + t * (0.407 + t * (-1.0 / 233.0 + t / 2547.0)) 1463 elif year >= 1961 and year < 1986: 1464 t = y - 1975.0 1465 dt = 45.45 + t * (1.067 + t * (-1.0 / 260.0 - t / 718.0)) 1466 elif year >= 1986 and year < 2005: 1467 t = y - 2000.0 1468 dt = 63.86 + t * ( 1469 0.3345 1470 + t 1471 * (-0.060374 + t * (0.0017275 + 1472 t * (0.000651814 + 0.00002373599 * t))) 1473 ) 1474 elif year >= 2005 and year < 2050: 1475 t = y - 2000.0 1476 dt = 62.92 + t * (0.32217 + 0.005589 * t) 1477 elif year >= 2050 and year < 2150: 1478 dt = (-20.0 + 32.0 * ((y - 1820.0) / 100.0) ** 2 - 1479 0.5628 * (2150.0 - y)) 1480 else: 1481 u = (year - 1820.0) / 100.0 1482 dt = -20.0 + 32.0 * u * u 1483 return dt 1484 1485 def dow(self, as_string=False): 1486 """Method to return the day of week corresponding to this Epoch. 1487 1488 By default, this method returns an integer value: 0 for Sunday, 1 for 1489 Monday, etc. However, when **as_string=True** is passed, the names of 1490 the days are returned. 1491 1492 :param as_string: Whether result will be given as a integer or as a 1493 string. False by default. 1494 :type as_string: bool 1495 1496 :returns: Day of the week, as a integer or as a string. 1497 :rtype: int, str 1498 1499 >>> e = Epoch(1954, 'June', 30) 1500 >>> e.dow() 1501 3 1502 >>> e = Epoch(2018, 'Feb', 14.9) 1503 >>> e.dow(as_string=True) 1504 'Wednesday' 1505 >>> e = Epoch(2018, 'Feb', 15) 1506 >>> e.dow(as_string=True) 1507 'Thursday' 1508 >>> e = Epoch(2018, 'Feb', 15.99) 1509 >>> e.dow(as_string=True) 1510 'Thursday' 1511 >>> e.set(2018, 'Jul', 15.4) 1512 >>> e.dow(as_string=True) 1513 'Sunday' 1514 >>> e.set(2018, 'Jul', 15.9) 1515 >>> e.dow(as_string=True) 1516 'Sunday' 1517 """ 1518 1519 jd = iint(self._jde - 0.5) + 2.0 1520 doy = iint(jd % 7) 1521 if not as_string: 1522 return doy 1523 else: 1524 day_names = [ 1525 "Sunday", 1526 "Monday", 1527 "Tuesday", 1528 "Wednesday", 1529 "Thursday", 1530 "Friday", 1531 "Saturday", 1532 ] 1533 return day_names[doy] 1534 1535 def mean_sidereal_time(self): 1536 """Method to compute the _mean_ sidereal time at Greenwich for the 1537 epoch stored in this object. It represents the Greenwich hour angle of 1538 the mean vernal equinox. 1539 1540 .. note:: If you require the result as an angle, you should convert the 1541 result from this method to hours with decimals (with 1542 :const:`DAY2HOURS`), and then multiply by 15 deg/hr. Alternatively, 1543 you can convert the result to hours with decimals, and feed this 1544 value to an :class:`Angle` object, setting **ra=True**, and making 1545 use of :class:`Angle` facilities for further handling. 1546 1547 :returns: Mean sidereal time, in days 1548 :rtype: float 1549 1550 >>> e = Epoch(1987, 4, 10) 1551 >>> round(e.mean_sidereal_time(), 9) 1552 0.549147764 1553 >>> e = Epoch(1987, 4, 10, 19, 21, 0.0) 1554 >>> round(e.mean_sidereal_time(), 9) 1555 0.357605204 1556 """ 1557 1558 jd0 = iint(self()) + 0.5 if self() % 1 >= 0.5 else iint(self()) - 0.5 1559 t = (jd0 - 2451545.0) / 36525.0 1560 theta0 = 6.0 / DAY2HOURS + 41.0 / DAY2MIN + 50.54841 / DAY2SEC 1561 s = t * (8640184.812866 + t * (0.093104 - 0.0000062 * t)) 1562 theta0 += (s % DAY2SEC) / DAY2SEC 1563 deltajd = self() - jd0 1564 if abs(deltajd) < TOL: # In this case, we are done 1565 return theta0 % 1 1566 else: 1567 deltajd *= 1.00273790935 1568 return (theta0 + deltajd) % 1 1569 1570 def apparent_sidereal_time(self, true_obliquity, nutation_longitude): 1571 """Method to compute the _apparent_ sidereal time at Greenwich for the 1572 epoch stored in this object. It represents the Greenwich hour angle of 1573 the true vernal equinox. 1574 1575 .. note:: If you require the result as an angle, you should convert the 1576 result from this method to hours with decimals (with 1577 :const:`DAY2HOURS`), and then multiply by 15 deg/hr. Alternatively, 1578 you can convert the result to hours with decimals, and feed this 1579 value to an :class:`Angle` object, setting **ra=True**, and making 1580 use of :class:`Angle` facilities for further handling. 1581 1582 :param true_obliquity: The true obliquity of the ecliptic as an int, 1583 float or :class:`Angle`, in degrees. You can use the method 1584 `Earth.true_obliquity()` to find it. 1585 :type true_obliquity: int, float, :class:`Angle` 1586 :param nutation_longitude: The nutation in longitude as an int, float 1587 or :class:`Angle`, in degrees. You can use method 1588 `Earth.nutation_longitude()` to find it. 1589 :type nutation_longitude: int, float, :class:`Angle` 1590 1591 :returns: Apparent sidereal time, in days 1592 :rtype: float 1593 :raises: TypeError if input value is of wrong type. 1594 1595 >>> e = Epoch(1987, 4, 10) 1596 >>> round(e.apparent_sidereal_time(23.44357, (-3.788)/3600.0), 8) 1597 0.54914508 1598 """ 1599 1600 if not ( 1601 isinstance(true_obliquity, (int, float, Angle)) 1602 and isinstance(nutation_longitude, (int, float, Angle)) 1603 ): 1604 raise TypeError("Invalid input value") 1605 if isinstance(true_obliquity, Angle): 1606 true_obliquity = float(true_obliquity) # Convert to a float 1607 if isinstance(nutation_longitude, Angle): 1608 nutation_longitude = float(nutation_longitude) 1609 mean_stime = self.mean_sidereal_time() 1610 epsilon = radians(true_obliquity) # Convert to radians 1611 delta_psi = nutation_longitude * 3600.0 # From degrees to seconds 1612 # Correction is in seconds of arc: It must be converted to seconds of 1613 # time, and then to days (sidereal time is given here in days) 1614 return mean_stime + ((delta_psi * cos(epsilon)) / 15.0) / DAY2SEC 1615 1616 def mjd(self): 1617 """This method returns the Modified Julian Day (MJD). 1618 1619 :returns: Modified Julian Day (MJD). 1620 :rtype: float 1621 1622 >>> e = Epoch(1858, 'NOVEMBER', 17) 1623 >>> e.mjd() 1624 0.0 1625 """ 1626 1627 return self._jde - 2400000.5 1628 1629 def jde(self): 1630 """Method to return the internal value of the Julian Ephemeris Day. 1631 1632 :returns: The internal value of the Julian Ephemeris Day. 1633 :rtype: float 1634 1635 >>> a = Epoch(-1000, 2, 29.0) 1636 >>> print(a.jde()) 1637 1355866.5 1638 """ 1639 1640 return self._jde 1641 1642 def year(self): 1643 """This method returns the contents of this object as a year with 1644 decimals. 1645 1646 :returns: Year with decimals. 1647 :rtype: float 1648 1649 >>> e = Epoch(1993, 'October', 1) 1650 >>> print(round(e.year(), 4)) 1651 1993.7479 1652 """ 1653 1654 y, m, d = self.get_date() 1655 doy = Epoch.get_doy(y, m, d) 1656 # We must substract 1 from doy in order to compute correctly 1657 doy -= 1 1658 days_of_year = 365.0 1659 if self.leap(): 1660 days_of_year = 366.0 1661 return y + doy / days_of_year 1662 1663 def rise_set(self, latitude, longitude, altitude=0.0): 1664 """This method computes the times of rising and setting of the Sun. 1665 1666 .. note:: The algorithm used is the one explained in the article 1667 "Sunrise equation" of the Wikipedia at: 1668 https://en.wikipedia.org/wiki/Sunrise_equation 1669 1670 .. note:: This algorithm is only valid within the artic and antartic 1671 circles (+/- 66d 33'). Outside that range this method returns 1672 invalid values (JDE < 0) 1673 1674 .. note:: The results are given in UTC time. 1675 1676 :param latitude: Latitude of the observer, as an Angle object. Positive 1677 to the North 1678 :type latitude: :py:class:`Angle` 1679 :param longitude: Longitude of the observer, as an Angle object. 1680 Positive to the East 1681 :type longitude: :py:class:`Angle` 1682 :param altitude: Altitude of the observer, as meters above sea level 1683 :type altitude: int, float 1684 1685 :returns: Two :py:class:`Epoch` objects representing rising time and 1686 setting time, in a tuple 1687 :rtype: tuple 1688 :raises: TypeError if input values are of wrong type. 1689 1690 >>> e = Epoch(2019, 4, 2) 1691 >>> latitude = Angle(48, 8, 0) 1692 >>> longitude = Angle(11, 34, 0) 1693 >>> altitude = 520.0 1694 >>> rising, setting = e.rise_set(latitude, longitude, altitude) 1695 >>> y, m, d, h, mi, s = rising.get_full_date() 1696 >>> print("{}:{}".format(h, mi)) 1697 4:48 1698 >>> y, m, d, h, mi, s = setting.get_full_date() 1699 >>> print("{}:{}".format(h, mi)) 1700 17:48 1701 """ 1702 1703 if not (isinstance(latitude, Angle) and isinstance(longitude, Angle) 1704 and isinstance(altitude, (int, float))): 1705 raise TypeError("Invalid input types") 1706 # Check that latitude is within valid range 1707 limit = Angle(66, 33, 0) 1708 if latitude > limit or latitude < -limit: 1709 error = Epoch(-1) 1710 return error, error 1711 # Let's start computing the number of days since 2000/1/1 12:00 (cjd) 1712 # Compute fractional Julian Day for leap seconds and terrestrial time 1713 # We need current epoch without hours, minutes and seconds 1714 year, month, day = self.get_date() 1715 e = Epoch(year, month, day) 1716 frac = (10.0 + 32.184 + Epoch.leap_seconds(year, month)) / 86400.0 1717 cjd = e.jde() - 2451545.0 + frac 1718 # Compute mean solar noon 1719 jstar = cjd - (float(longitude) / 360.0) 1720 # Solar mean anomaly 1721 m = (357.5291 + 0.98560028 * jstar) % 360 1722 mr = radians(m) 1723 # Equation of the center 1724 c = 1.9148 * sin(mr) + 0.02 * sin(2.0 * mr) + 0.0003 * sin(3.0 * mr) 1725 # Ecliptic longitude 1726 lambd = (m + c + 180.0 + 102.9372) % 360 1727 lr = radians(lambd) 1728 # Solar transit 1729 jtran = 2451545.5 + jstar + 0.0053 * sin(mr) - 0.0069 * sin(2.0 * lr) 1730 # NOTE: The original algorithm indicates a value of 2451545.0, but that 1731 # leads to transit times around midnight, which is an error 1732 # Declination of the Sun 1733 sin_delta = sin(lr) * sin(radians(23.44)) 1734 delta = asin(sin_delta) 1735 cos_delta = cos(delta) 1736 # Hour angle 1737 # First, correct by elevation 1738 corr = -0.83 - 2.076 * sqrt(altitude) / 60.0 1739 cos_om = ((sin(radians(corr)) - sin(latitude.rad()) * sin_delta) / 1740 (cos(latitude.rad()) * cos_delta)) 1741 # Finally, compute rising and setting times 1742 omega = degrees(acos(cos_om)) 1743 jrise = Epoch(jtran - (omega / 360.0)) 1744 jsett = Epoch(jtran + (omega / 360.0)) 1745 return jrise, jsett 1746 1747 def __call__(self): 1748 """Method used when Epoch is called only with parenthesis. 1749 1750 :returns: The internal value of the Julian Ephemeris Day. 1751 :rtype: float 1752 1753 >>> a = Epoch(-122, 1, 1.0) 1754 >>> print(a()) 1755 1676497.5 1756 """ 1757 1758 return self._jde 1759 1760 def __add__(self, b): 1761 """This method defines the addition between an Epoch and some days. 1762 1763 :param b: Value to be added, in days. 1764 :type b: int, float 1765 1766 :returns: A new Epoch object. 1767 :rtype: :py:class:`Epoch` 1768 :raises: TypeError if operand is of wrong type. 1769 1770 >>> a = Epoch(1991, 7, 11) 1771 >>> b = a + 10000 1772 >>> y, m, d = b.get_date() 1773 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1774 2018/11/26.0 1775 """ 1776 1777 if isinstance(b, (int, float)): 1778 return Epoch(self._jde + float(b)) 1779 else: 1780 raise TypeError("Wrong operand type") 1781 1782 def __sub__(self, b): 1783 """This method defines the subtraction between Epochs or between an 1784 Epoch and a given number of days. 1785 1786 :param b: Value to be subtracted, either an Epoch or days. 1787 :type b: py:class:`Epoch`, int, float 1788 1789 :returns: A new Epoch object if parameter 'b' is in days, or the 1790 difference between provided Epochs, in days. 1791 :rtype: :py:class:`Epoch`, float 1792 :raises: TypeError if operand is of wrong type. 1793 1794 >>> a = Epoch(1986, 2, 9.0) 1795 >>> print(round(a(), 2)) 1796 2446470.5 1797 >>> b = Epoch(1910, 4, 20.0) 1798 >>> print(round(b(), 2)) 1799 2418781.5 1800 >>> c = a - b 1801 >>> print(round(c, 2)) 1802 27689.0 1803 >>> a = Epoch(2003, 12, 31.0) 1804 >>> b = a - 365.5 1805 >>> y, m, d = b.get_date() 1806 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1807 2002/12/30.5 1808 """ 1809 1810 if isinstance(b, (int, float)): 1811 return Epoch(self._jde - b) 1812 elif isinstance(b, Epoch): 1813 return float(self._jde - b._jde) 1814 else: 1815 raise TypeError("Invalid operand type") 1816 1817 def __iadd__(self, b): 1818 """This method defines the accumulative addition to this Epoch. 1819 1820 :param b: Value to be added, in days. 1821 :type b: int, float 1822 1823 :returns: This Epoch. 1824 :rtype: :py:class:`Epoch` 1825 :raises: TypeError if operand is of wrong type. 1826 1827 >>> a = Epoch(2003, 12, 31.0) 1828 >>> a += 32.5 1829 >>> y, m, d = a.get_date() 1830 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1831 2004/2/1.5 1832 """ 1833 1834 if isinstance(b, (int, float)): 1835 self = self + b 1836 return self 1837 else: 1838 raise TypeError("Wrong operand type") 1839 1840 def __isub__(self, b): 1841 """This method defines the accumulative subtraction to this Epoch. 1842 1843 :param b: Value to be subtracted, in days. 1844 :type b: int, float 1845 1846 :returns: This Epoch. 1847 :rtype: :py:class:`Epoch` 1848 :raises: TypeError if operand is of wrong type. 1849 1850 >>> a = Epoch(2001, 12, 31.0) 1851 >>> a -= 2*365 1852 >>> y, m, d = a.get_date() 1853 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1854 2000/1/1.0 1855 """ 1856 1857 if isinstance(b, (int, float)): 1858 self = self - b 1859 return self 1860 else: 1861 raise TypeError("Wrong operand type") 1862 1863 def __radd__(self, b): 1864 """This method defines the addition to a Epoch by the right 1865 1866 :param b: Value to be added, in days. 1867 :type b: int, float 1868 1869 :returns: A new Epoch object. 1870 :rtype: :py:class:`Epoch` 1871 :raises: TypeError if operand is of wrong type. 1872 1873 >>> a = Epoch(2004, 2, 27.8) 1874 >>> b = 2.2 + a 1875 >>> y, m, d = b.get_date() 1876 >>> print("{}/{}/{}".format(y, m, round(d, 2))) 1877 2004/3/1.0 1878 """ 1879 1880 if isinstance(b, (int, float)): 1881 return self.__add__(b) # It is the same as by the left 1882 else: 1883 raise TypeError("Wrong operand type") 1884 1885 def __int__(self): 1886 """This method returns the internal JDE value as an int. 1887 1888 :returns: Internal JDE value as an int. 1889 :rtype: int 1890 1891 >>> a = Epoch(2434923.85) 1892 >>> int(a) 1893 2434923 1894 """ 1895 1896 return int(self._jde) 1897 1898 def __float__(self): 1899 """This method returns the internal JDE value as a float. 1900 1901 :returns: Internal JDE value as a float. 1902 :rtype: float 1903 1904 >>> a = Epoch(2434923.85) 1905 >>> float(a) 1906 2434923.85 1907 """ 1908 1909 return float(self._jde) 1910 1911 def __eq__(self, b): 1912 """This method defines the 'is equal' operator between Epochs. 1913 1914 .. note:: For the comparison, the **base.TOL** value is used. 1915 1916 :returns: A boolean. 1917 :rtype: bool 1918 :raises: TypeError if input values are of wrong type. 1919 1920 >>> a = Epoch(2007, 5, 20.0) 1921 >>> b = Epoch(2007, 5, 20.000001) 1922 >>> a == b 1923 False 1924 >>> a = Epoch(2004, 10, 15.7) 1925 >>> b = Epoch(a) 1926 >>> a == b 1927 True 1928 >>> a = Epoch(2434923.85) 1929 >>> a == 2434923.85 1930 True 1931 """ 1932 1933 if isinstance(b, (int, float)): 1934 return abs(self._jde - float(b)) < TOL 1935 elif isinstance(b, Epoch): 1936 return abs(self._jde - b._jde) < TOL 1937 else: 1938 raise TypeError("Wrong operand type") 1939 1940 def __ne__(self, b): 1941 """This method defines the 'is not equal' operator between Epochs. 1942 1943 .. note:: For the comparison, the **base.TOL** value is used. 1944 1945 :returns: A boolean. 1946 :rtype: bool 1947 1948 >>> a = Epoch(2007, 5, 20.0) 1949 >>> b = Epoch(2007, 5, 20.000001) 1950 >>> a != b 1951 True 1952 >>> a = Epoch(2004, 10, 15.7) 1953 >>> b = Epoch(a) 1954 >>> a != b 1955 False 1956 >>> a = Epoch(2434923.85) 1957 >>> a != 2434923.85 1958 False 1959 """ 1960 1961 return not self.__eq__(b) # '!=' == 'not(==)' 1962 1963 def __lt__(self, b): 1964 """This method defines the 'is less than' operator between Epochs. 1965 1966 :returns: A boolean. 1967 :rtype: bool 1968 :raises: TypeError if input values are of wrong type. 1969 1970 >>> a = Epoch(2004, 10, 15.7) 1971 >>> b = Epoch(2004, 10, 15.7) 1972 >>> a < b 1973 False 1974 """ 1975 1976 if isinstance(b, (int, float)): 1977 return self._jde < float(b) 1978 elif isinstance(b, Epoch): 1979 return self._jde < b._jde 1980 else: 1981 raise TypeError("Wrong operand type") 1982 1983 def __ge__(self, b): 1984 """This method defines 'is equal or greater' operator between Epochs. 1985 1986 :returns: A boolean. 1987 :rtype: bool 1988 :raises: TypeError if input values are of wrong type. 1989 1990 >>> a = Epoch(2004, 10, 15.71) 1991 >>> b = Epoch(2004, 10, 15.7) 1992 >>> a >= b 1993 True 1994 """ 1995 1996 return not self.__lt__(b) # '>=' == 'not(<)' 1997 1998 def __gt__(self, b): 1999 """This method defines the 'is greater than' operator between Epochs. 2000 2001 :returns: A boolean. 2002 :rtype: bool 2003 :raises: TypeError if input values are of wrong type. 2004 2005 >>> a = Epoch(2004, 10, 15.71) 2006 >>> b = Epoch(2004, 10, 15.7) 2007 >>> a > b 2008 True 2009 >>> a = Epoch(-207, 11, 5.2) 2010 >>> b = Epoch(-207, 11, 5.2) 2011 >>> a > b 2012 False 2013 """ 2014 2015 if isinstance(b, (int, float)): 2016 return self._jde > float(b) 2017 elif isinstance(b, Epoch): 2018 return self._jde > b._jde 2019 else: 2020 raise TypeError("Wrong operand type") 2021 2022 def __le__(self, b): 2023 """This method defines 'is equal or less' operator between Epochs. 2024 2025 :returns: A boolean. 2026 :rtype: bool 2027 :raises: TypeError if input values are of wrong type. 2028 2029 >>> a = Epoch(2004, 10, 15.71) 2030 >>> b = Epoch(2004, 10, 15.7) 2031 >>> a <= b 2032 False 2033 >>> a = Epoch(-207, 11, 5.2) 2034 >>> b = Epoch(-207, 11, 5.2) 2035 >>> a <= b 2036 True 2037 """ 2038 2039 return not self.__gt__(b) # '<=' == 'not(>)' 2040 2041 2042JDE2000 = Epoch(2000, 1, 1.5) 2043"""Standard epoch for January 1st, 2000 at 12h corresponding to JDE2451545.0""" 2044 2045 2046def main(): 2047 2048 # Let's define a small helper function 2049 def print_me(msg, val): 2050 print("{}: {}".format(msg, val)) 2051 2052 # Let's do some work with the Epoch class 2053 print("\n" + 35 * "*") 2054 print("*** Use of Epoch class") 2055 print(35 * "*" + "\n") 2056 2057 e = Epoch(1987, 6, 19.5) 2058 print_me("JDE for 1987/6/19.5", e) 2059 2060 # Redefine the Epoch object 2061 e.set(333, "Jan", 27, 12) 2062 print_me("JDE for 333/1/27.5", e) 2063 2064 # We can create an Epoch from a 'date' or 'datetime' object 2065 d = datetime.datetime(837, 4, 10, 7, 12, 0, 0) 2066 f = Epoch(d) 2067 print_me("JDE for 837/4/10.3", f) 2068 2069 print("") 2070 2071 # Check if a given date belong to the Julian or Gregorian calendar 2072 print_me("Is 1590/4/21.4 a Julian date?", Epoch.is_julian(1590, 4, 21.4)) 2073 2074 print("") 2075 2076 # We can also check if a given year is leap or not 2077 print_me("Is -1000 a leap year?", Epoch.is_leap(-1000)) 2078 print_me("Is 1800 a leap year?", Epoch.is_leap(1800)) 2079 print_me("Is 2012 a leap year?", Epoch.is_leap(2012)) 2080 2081 print("") 2082 2083 # Get the Day Of Year corresponding to a given date 2084 print_me("Day Of Year (DOY) of 1978/11/14", Epoch.get_doy(1978, 11, 14)) 2085 print_me("Day Of Year (DOY) of -400/2/29.9", Epoch.get_doy(-400, 2, 29.9)) 2086 2087 print("") 2088 2089 # Now the opposite: Get a date from a DOY 2090 t = Epoch.doy2date(2017, 365.7) 2091 s = str(t[0]) + "/" + str(t[1]) + "/" + str(round(t[2], 2)) 2092 print_me("Date from DOY 2017:365.7", s) 2093 2094 t = Epoch.doy2date(-4, 60) 2095 s = str(t[0]) + "/" + str(t[1]) + "/" + str(round(t[2], 2)) 2096 print_me("Date from DOY -4:60", s) 2097 2098 print("") 2099 2100 # There is an internal table which we can use to get the leap seconds 2101 print_me("Number of leap seconds applied up to July 1983", 2102 Epoch.leap_seconds(1983, 7)) 2103 2104 print("") 2105 2106 # We can convert the internal JDE value back to a date 2107 e = Epoch(2436116.31) 2108 y, m, d = e.get_date() 2109 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2110 print_me("Date from JDE 2436116.31", s) 2111 2112 print("") 2113 2114 # It is possible to get the day of the week corresponding to a given date 2115 e = Epoch(2018, "Feb", 15) 2116 print_me("The day of week of 2018/2/15 is", e.dow(as_string=True)) 2117 2118 print("") 2119 2120 # In some cases it is useful to get the Modified Julian Day (MJD) 2121 e = Epoch(1923, "August", 23) 2122 print_me("Modified Julian Day for 1923/8/23", round(e.mjd(), 2)) 2123 2124 print("") 2125 2126 # If your system is appropriately configured, you can get the difference in 2127 # seconds between your local time and UTC 2128 print_me( 2129 "To convert from local system time to UTC you must add/subtract" 2130 + " this amount of seconds", 2131 Epoch.utc2local(), 2132 ) 2133 2134 print("") 2135 2136 # Compute DeltaT = TT - UT differences for various dates 2137 print_me("DeltaT (TT - UT) for Feb/333", round(Epoch.tt2ut(333, 2), 1)) 2138 print_me("DeltaT (TT - UT) for Jan/1642", round(Epoch.tt2ut(1642, 1), 1)) 2139 print_me("DeltaT (TT - UT) for Feb/1928", round(Epoch.tt2ut(1928, 1), 1)) 2140 print_me("DeltaT (TT - UT) for Feb/1977", round(Epoch.tt2ut(1977, 2), 1)) 2141 print_me("DeltaT (TT - UT) for Jan/1998", round(Epoch.tt2ut(1998, 1), 1)) 2142 2143 print("") 2144 2145 # The difference between civil day and sidereal day is almost 4 minutes 2146 e = Epoch(1987, 4, 10) 2147 st1 = round(e.mean_sidereal_time(), 9) 2148 e = Epoch(1987, 4, 11) 2149 st2 = round(e.mean_sidereal_time(), 9) 2150 ds = (st2 - st1) * DAY2MIN 2151 msg = "{}m {}s".format(iint(ds), (ds % 1) * 60.0) 2152 print_me("Difference between sidereal time 1987/4/11 and 1987/4/10", msg) 2153 2154 print("") 2155 2156 print( 2157 "When correcting for nutation-related effects, we get the " 2158 + "'apparent' sidereal time:" 2159 ) 2160 e = Epoch(1987, 4, 10) 2161 print("e = Epoch(1987, 4, 10)") 2162 print_me( 2163 "e.apparent_sidereal_time(23.44357, (-3.788)/3600.0)", 2164 e.apparent_sidereal_time(23.44357, (-3.788) / 3600.0), 2165 ) 2166 # 0.549145082637 2167 2168 print("") 2169 2170 # Epoch class can also provide the date of Easter for a given year 2171 # Let's spice up the output a little bit, calling dow() and get_month() 2172 month, day = Epoch.easter(2019) 2173 e = Epoch(2019, month, day) 2174 s = ( 2175 e.dow(as_string=True) 2176 + ", " 2177 + str(day) 2178 + get_ordinal_suffix(day) 2179 + " of " 2180 + Epoch.get_month(month, as_string=True) 2181 ) 2182 print_me("Easter day for 2019", s) 2183 # I know Easter is always on Sunday, by the way... ;-) 2184 2185 print("") 2186 2187 # Compute the date of the Jewish Easter (Pesach) for a given year 2188 month, day = Epoch.jewish_pesach(1990) 2189 s = ( 2190 str(day) 2191 + get_ordinal_suffix(day) 2192 + " of " 2193 + Epoch.get_month(month, as_string=True) 2194 ) 2195 print_me("Jewish Pesach day for 1990", s) 2196 2197 print("") 2198 2199 # Now, convert a date in the Moslem calendar to the Gregorian calendar 2200 y, m, d = Epoch.moslem2gregorian(1421, 1, 1) 2201 print_me("The date 1421/1/1 in the Moslem calendar is, in Gregorian " + 2202 "calendar", "{}/{}/{}".format(y, m, d)) 2203 y, m, d = Epoch.moslem2gregorian(1439, 9, 1) 2204 print_me( 2205 "The start of Ramadan month (9/1) for Gregorian year 2018 is", 2206 "{}/{}/{}".format(y, m, d), 2207 ) 2208 # We can go from the Gregorian calendar back to the Moslem calendar too 2209 print_me( 2210 "Date 1991/8/13 in Gregorian calendar is, in Moslem calendar", 2211 "{}/{}/{}".format(*Epoch.gregorian2moslem(1991, 8, 13)), 2212 ) 2213 # Note: The '*' before 'Epoch' will _unpack_ the tuple into components 2214 2215 print("") 2216 2217 # It is possible to carry out some algebraic operations with Epochs 2218 2219 # Add 10000 days to a given date 2220 a = Epoch(1991, 7, 11) 2221 b = a + 10000 2222 y, m, d = b.get_date() 2223 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2224 print_me("1991/7/11 plus 10000 days is", s) 2225 2226 # Subtract two Epochs to find the number of days between them 2227 a = Epoch(1986, 2, 9.0) 2228 b = Epoch(1910, 4, 20.0) 2229 print_me("The number of days between 1986/2/9 and 1910/4/20 is", 2230 round(a - b, 2)) 2231 2232 # We can also subtract a given amount of days from an Epoch 2233 a = Epoch(2003, 12, 31.0) 2234 b = a - 365.5 2235 y, m, d = b.get_date() 2236 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2237 print_me("2003/12/31 minus 365.5 days is", s) 2238 2239 # Accumulative addition and subtraction of days is also allowed 2240 a = Epoch(2003, 12, 31.0) 2241 a += 32.5 2242 y, m, d = a.get_date() 2243 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2244 print_me("2003/12/31 plus 32.5 days is", s) 2245 2246 a = Epoch(2001, 12, 31.0) 2247 a -= 2 * 365 2248 y, m, d = a.get_date() 2249 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2250 print_me("2001/12/31 minus 2*365 days is", s) 2251 2252 # It is also possible to add days from the right 2253 a = Epoch(2004, 2, 27.8) 2254 b = 2.2 + a 2255 y, m, d = b.get_date() 2256 s = str(y) + "/" + str(m) + "/" + str(round(d, 2)) 2257 print_me("2004/2/27.8 plus 2.2 days is", s) 2258 2259 print("") 2260 2261 # Comparison operadors between epochs are also defined 2262 a = Epoch(2007, 5, 20.0) 2263 b = Epoch(2007, 5, 20.000001) 2264 print_me("2007/5/20.0 == 2007/5/20.000001", a == b) 2265 print_me("2007/5/20.0 != 2007/5/20.000001", a != b) 2266 print_me("2007/5/20.0 > 2007/5/20.000001", a > b) 2267 print_me("2007/5/20.0 <= 2007/5/20.000001", a <= b) 2268 2269 print("") 2270 2271 # Compute the time of rise and setting of the Sun in a given day 2272 e = Epoch(2018, 5, 2) 2273 print("On May 2nd, 2018, Sun rising/setting times in Munich were (UTC):") 2274 latitude = Angle(48, 8, 0) 2275 longitude = Angle(11, 34, 0) 2276 altitude = 520.0 2277 rising, setting = e.rise_set(latitude, longitude, altitude) 2278 y, m, d, h, mi, s = rising.get_full_date() 2279 print("Rising time: {}:{}".format(h, mi)) # 3:50 2280 y, m, d, h, mi, s = setting.get_full_date() 2281 print("Setting time: {}:{}".format(h, mi)) # 18:33 2282 2283 2284if __name__ == "__main__": 2285 2286 main() 2287