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.
 
 
 
 

166 lines
6.3 KiB

  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 skyfield import almanac
  19. from skyfield.timelib import Time
  20. from skyfield.constants import tau
  21. from .data import Object, Position, AsterEphemerides, MoonPhase
  22. from .core import get_skf_objects, get_timescale, get_iau2000b, ASTERS, MONTHS, skyfield_to_moon_phase
  23. RISEN_ANGLE = -0.8333
  24. class EphemeridesComputer:
  25. def __init__(self, position: Position):
  26. position.observation_planet = get_skf_objects()['earth']
  27. self.position = position
  28. def get_sun(self, start_time, end_time) -> dict:
  29. times, is_risen = almanac.find_discrete(start_time,
  30. end_time,
  31. almanac.sunrise_sunset(get_skf_objects(), self.position))
  32. sunrise = times[0] if is_risen[0] else times[1]
  33. sunset = times[1] if not is_risen[1] else times[0]
  34. return {'rise': sunrise, 'set': sunset}
  35. @staticmethod
  36. def get_moon_phase(year, month, day) -> MoonPhase:
  37. earth = get_skf_objects()['earth']
  38. moon = get_skf_objects()['moon']
  39. sun = get_skf_objects()['sun']
  40. def moon_phase_at(time: Time):
  41. time._nutation_angles = get_iau2000b(time)
  42. current_earth = earth.at(time)
  43. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon('date')
  44. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon('date')
  45. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  46. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  47. today = get_timescale().utc(year, month, day)
  48. time1 = get_timescale().utc(year, month, day - 10)
  49. time2 = get_timescale().utc(year, month, day + 10)
  50. times, phase = almanac.find_discrete(time1, time2, moon_phase_at)
  51. return skyfield_to_moon_phase(times, phase, today)
  52. @staticmethod
  53. def get_asters_ephemerides_for_aster(aster, date: datetime.date, position: Position) -> Object:
  54. skyfield_aster = get_skf_objects()[aster.skyfield_name]
  55. def get_angle(time: Time) -> float:
  56. return position.get_planet_topos().at(time).observe(skyfield_aster).apparent().altaz()[0].degrees
  57. def is_risen(time: Time) -> bool:
  58. return get_angle(time) > RISEN_ANGLE
  59. get_angle.rough_period = 1.0
  60. is_risen.rough_period = 0.5
  61. start_time = get_timescale().utc(date.year, date.month, date.day)
  62. end_time = get_timescale().utc(date.year, date.month, date.day, 23, 59, 59)
  63. rise_times, arr = almanac.find_discrete(start_time, end_time, is_risen)
  64. try:
  65. culmination_time, _ = almanac._find_maxima(start_time, end_time, get_angle, epsilon=1./3600/24)
  66. except ValueError:
  67. culmination_time = None
  68. if len(rise_times) == 2:
  69. rise_time = rise_times[0 if arr[0] else 1]
  70. set_time = rise_times[1 if not arr[1] else 0]
  71. else:
  72. rise_time = rise_times[0] if arr[0] else None
  73. set_time = rise_times[0] if not arr[0] else None
  74. culmination_time = culmination_time[0] if culmination_time is not None else None
  75. aster.ephemerides = AsterEphemerides(rise_time, culmination_time, set_time)
  76. return aster
  77. @staticmethod
  78. def is_leap_year(year: int) -> bool:
  79. return (year % 4 == 0 and year % 100 > 0) or (year % 400 == 0)
  80. def compute_ephemerides_for_day(self, year: int, month: int, day: int) -> dict:
  81. return {'moon_phase': self.get_moon_phase(year, month, day),
  82. 'details': [self.get_asters_ephemerides_for_aster(aster, datetime.date(year, month, day), self.position)
  83. for aster in ASTERS]}
  84. def compute_ephemerides_for_month(self, year: int, month: int) -> [dict]:
  85. if month == 2:
  86. max_day = 29 if self.is_leap_year(year) else 28
  87. elif month < 8:
  88. max_day = 30 if month % 2 == 0 else 31
  89. else:
  90. max_day = 31 if month % 2 == 0 else 30
  91. ephemerides = []
  92. for day in range(1, max_day + 1):
  93. ephemerides.append(self.compute_ephemerides_for_day(year, month, day))
  94. return ephemerides
  95. def compute_ephemerides_for_year(self, year: int) -> [dict]:
  96. ephemerides = {'seasons': self.get_seasons(year)}
  97. for month in range(0, 12):
  98. ephemerides[MONTHS[month]] = self.compute_ephemerides_for_month(year, month + 1)
  99. return ephemerides
  100. @staticmethod
  101. def get_seasons(year: int) -> dict:
  102. start_time = get_timescale().utc(year, 1, 1)
  103. end_time = get_timescale().utc(year, 12, 31)
  104. times, almanac_seasons = almanac.find_discrete(start_time, end_time, almanac.seasons(get_skf_objects()))
  105. seasons = {}
  106. for time, almanac_season in zip(times, almanac_seasons):
  107. if almanac_season == 0:
  108. season = 'MARCH'
  109. elif almanac_season == 1:
  110. season = 'JUNE'
  111. elif almanac_season == 2:
  112. season = 'SEPTEMBER'
  113. elif almanac_season == 3:
  114. season = 'DECEMBER'
  115. else:
  116. raise AssertionError
  117. seasons[season] = time.utc_iso()
  118. return seasons
  119. def compute_ephemerides(self, year: int, month: int, day: int):
  120. if day is not None:
  121. return self.compute_ephemerides_for_day(year, month, day)
  122. if month is not None:
  123. return self.compute_ephemerides_for_month(year, month)
  124. return self.compute_ephemerides_for_year(year)