Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.
 
 
 
 

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