You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

ephemerides.py 6.4 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167
  1. #!/usr/bin/env python3
  2. # Kosmorro - Compute The Next Ephemerides
  3. # Copyright (C) 2019 Jérôme Deuchnord <jerome@deuchnord.fr>
  4. #
  5. # This program is free software: you can redistribute it and/or modify
  6. # it under the terms of the GNU Affero General Public License as
  7. # published by the Free Software Foundation, either version 3 of the
  8. # License, or (at your option) any later version.
  9. #
  10. # This program is distributed in the hope that it will be useful,
  11. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13. # GNU Affero General Public License for more details.
  14. #
  15. # You should have received a copy of the GNU Affero General Public License
  16. # along with this program. If not, see <https://www.gnu.org/licenses/>.
  17. import datetime
  18. from typing import Union
  19. from skyfield import almanac
  20. from skyfield.timelib import Time
  21. from skyfield.constants import tau
  22. from .data import Object, Position, AsterEphemerides, MoonPhase
  23. from .core import get_skf_objects, get_timescale, get_iau2000b, ASTERS, MONTHS, skyfield_to_moon_phase
  24. RISEN_ANGLE = -0.8333
  25. class EphemeridesComputer:
  26. def __init__(self, position: Union[Position, None]):
  27. if position is not None:
  28. position.observation_planet = get_skf_objects()['earth']
  29. self.position = position
  30. def get_sun(self, start_time, end_time) -> dict:
  31. times, is_risen = almanac.find_discrete(start_time,
  32. end_time,
  33. almanac.sunrise_sunset(get_skf_objects(), self.position))
  34. sunrise = times[0] if is_risen[0] else times[1]
  35. sunset = times[1] if not is_risen[1] else times[0]
  36. return {'rise': sunrise, 'set': sunset}
  37. @staticmethod
  38. def get_moon_phase(year, month, day) -> MoonPhase:
  39. earth = get_skf_objects()['earth']
  40. moon = get_skf_objects()['moon']
  41. sun = get_skf_objects()['sun']
  42. def moon_phase_at(time: Time):
  43. time._nutation_angles = get_iau2000b(time)
  44. current_earth = earth.at(time)
  45. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon('date')
  46. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon('date')
  47. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  48. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  49. today = get_timescale().utc(year, month, day)
  50. time1 = get_timescale().utc(year, month, day - 10)
  51. time2 = get_timescale().utc(year, month, day + 10)
  52. times, phase = almanac.find_discrete(time1, time2, moon_phase_at)
  53. return skyfield_to_moon_phase(times, phase, today)
  54. @staticmethod
  55. def get_asters_ephemerides_for_aster(aster, date: datetime.date, position: Position) -> Object:
  56. skyfield_aster = get_skf_objects()[aster.skyfield_name]
  57. def get_angle(time: Time) -> float:
  58. return position.get_planet_topos().at(time).observe(skyfield_aster).apparent().altaz()[0].degrees
  59. def is_risen(time: Time) -> bool:
  60. return get_angle(time) > RISEN_ANGLE
  61. get_angle.rough_period = 1.0
  62. is_risen.rough_period = 0.5
  63. start_time = get_timescale().utc(date.year, date.month, date.day)
  64. end_time = get_timescale().utc(date.year, date.month, date.day, 23, 59, 59)
  65. rise_times, arr = almanac.find_discrete(start_time, end_time, is_risen)
  66. try:
  67. culmination_time, _ = almanac._find_maxima(start_time, end_time, get_angle, epsilon=1./3600/24)
  68. except ValueError:
  69. culmination_time = None
  70. if len(rise_times) == 2:
  71. rise_time = rise_times[0 if arr[0] else 1]
  72. set_time = rise_times[1 if not arr[1] else 0]
  73. else:
  74. rise_time = rise_times[0] if arr[0] else None
  75. set_time = rise_times[0] if not arr[0] else None
  76. culmination_time = culmination_time[0] if culmination_time is not None else None
  77. aster.ephemerides = AsterEphemerides(rise_time, culmination_time, set_time)
  78. return aster
  79. @staticmethod
  80. def is_leap_year(year: int) -> bool:
  81. return (year % 4 == 0 and year % 100 > 0) or (year % 400 == 0)
  82. def compute_ephemerides_for_day(self, year: int, month: int, day: int) -> dict:
  83. return {'moon_phase': self.get_moon_phase(year, month, day),
  84. 'details': [self.get_asters_ephemerides_for_aster(aster, datetime.date(year, month, day), self.position)
  85. for aster in ASTERS] if self.position is not None else []}
  86. def compute_ephemerides_for_month(self, year: int, month: int) -> [dict]:
  87. if month == 2:
  88. max_day = 29 if self.is_leap_year(year) else 28
  89. elif month < 8:
  90. max_day = 30 if month % 2 == 0 else 31
  91. else:
  92. max_day = 31 if month % 2 == 0 else 30
  93. ephemerides = []
  94. for day in range(1, max_day + 1):
  95. ephemerides.append(self.compute_ephemerides_for_day(year, month, day))
  96. return ephemerides
  97. def compute_ephemerides_for_year(self, year: int) -> [dict]:
  98. ephemerides = {'seasons': self.get_seasons(year)}
  99. for month in range(0, 12):
  100. ephemerides[MONTHS[month]] = self.compute_ephemerides_for_month(year, month + 1)
  101. return ephemerides
  102. @staticmethod
  103. def get_seasons(year: int) -> dict:
  104. start_time = get_timescale().utc(year, 1, 1)
  105. end_time = get_timescale().utc(year, 12, 31)
  106. times, almanac_seasons = almanac.find_discrete(start_time, end_time, almanac.seasons(get_skf_objects()))
  107. seasons = {}
  108. for time, almanac_season in zip(times, almanac_seasons):
  109. if almanac_season == 0:
  110. season = 'MARCH'
  111. elif almanac_season == 1:
  112. season = 'JUNE'
  113. elif almanac_season == 2:
  114. season = 'SEPTEMBER'
  115. elif almanac_season == 3:
  116. season = 'DECEMBER'
  117. else:
  118. raise AssertionError
  119. seasons[season] = time.utc_iso()
  120. return seasons
  121. def compute_ephemerides(self, year: int, month: int, day: int):
  122. if day is not None:
  123. return self.compute_ephemerides_for_day(year, month, day)
  124. if month is not None:
  125. return self.compute_ephemerides_for_month(year, month)
  126. return self.compute_ephemerides_for_year(year)