A library that computes the ephemerides.
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.
 
 

263 lines
8.6 KiB

  1. #!/usr/bin/env python3
  2. from datetime import date as date_type
  3. from skyfield.errors import EphemerisRangeError
  4. from skyfield.timelib import Time
  5. from skyfield.searchlib import find_discrete, find_maxima, find_minima
  6. from numpy import pi
  7. from .model import Event, Star, Planet, ASTERS
  8. from .dateutil import translate_to_timezone
  9. from .enum import EventType
  10. from .exceptions import OutOfRangeDateError
  11. from .core import get_timescale, get_skf_objects, flatten_list
  12. def _search_conjunction(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  13. earth = get_skf_objects()["earth"]
  14. aster1 = None
  15. aster2 = None
  16. def is_in_conjunction(time: Time):
  17. earth_pos = earth.at(time)
  18. _, aster1_lon, _ = (
  19. earth_pos.observe(aster1.get_skyfield_object()).apparent().ecliptic_latlon()
  20. )
  21. _, aster2_lon, _ = (
  22. earth_pos.observe(aster2.get_skyfield_object()).apparent().ecliptic_latlon()
  23. )
  24. return ((aster1_lon.radians - aster2_lon.radians) / pi % 2.0).astype(
  25. "int8"
  26. ) == 0
  27. is_in_conjunction.rough_period = 60.0
  28. computed = []
  29. conjunctions = []
  30. for aster1 in ASTERS:
  31. # Ignore the Sun
  32. if isinstance(aster1, Star):
  33. continue
  34. for aster2 in ASTERS:
  35. if isinstance(aster2, Star) or aster2 == aster1 or aster2 in computed:
  36. continue
  37. times, is_conjs = find_discrete(start_time, end_time, is_in_conjunction)
  38. for i, time in enumerate(times):
  39. if is_conjs[i]:
  40. aster1_pos = (aster1.get_skyfield_object() - earth).at(time)
  41. aster2_pos = (aster2.get_skyfield_object() - earth).at(time)
  42. distance = aster1_pos.separation_from(aster2_pos).degrees
  43. if distance - aster2.get_apparent_radius(
  44. time, earth
  45. ) < aster1.get_apparent_radius(time, earth):
  46. occulting_aster = (
  47. [aster1, aster2]
  48. if aster1_pos.distance().km < aster2_pos.distance().km
  49. else [aster2, aster1]
  50. )
  51. conjunctions.append(
  52. Event(
  53. EventType.OCCULTATION,
  54. occulting_aster,
  55. translate_to_timezone(time.utc_datetime(), timezone),
  56. )
  57. )
  58. else:
  59. conjunctions.append(
  60. Event(
  61. EventType.CONJUNCTION,
  62. [aster1, aster2],
  63. translate_to_timezone(time.utc_datetime(), timezone),
  64. )
  65. )
  66. computed.append(aster1)
  67. return conjunctions
  68. def _search_oppositions(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  69. earth = get_skf_objects()["earth"]
  70. sun = get_skf_objects()["sun"]
  71. aster = None
  72. def is_oppositing(time: Time) -> [bool]:
  73. earth_pos = earth.at(time)
  74. sun_pos = earth_pos.observe(
  75. sun
  76. ).apparent() # Never do this without eyes protection!
  77. aster_pos = earth_pos.observe(get_skf_objects()[aster.skyfield_name]).apparent()
  78. _, lon1, _ = sun_pos.ecliptic_latlon()
  79. _, lon2, _ = aster_pos.ecliptic_latlon()
  80. return (lon1.degrees - lon2.degrees) > 180
  81. is_oppositing.rough_period = 1.0
  82. events = []
  83. for aster in ASTERS:
  84. if not isinstance(aster, Planet) or aster.skyfield_name in ["MERCURY", "VENUS"]:
  85. continue
  86. times, _ = find_discrete(start_time, end_time, is_oppositing)
  87. for time in times:
  88. events.append(
  89. Event(
  90. EventType.OPPOSITION,
  91. [aster],
  92. translate_to_timezone(time.utc_datetime(), timezone),
  93. )
  94. )
  95. return events
  96. def _search_maximal_elongations(
  97. start_time: Time, end_time: Time, timezone: int
  98. ) -> [Event]:
  99. earth = get_skf_objects()["earth"]
  100. sun = get_skf_objects()["sun"]
  101. aster = None
  102. def get_elongation(time: Time):
  103. sun_pos = (sun - earth).at(time)
  104. aster_pos = (aster.get_skyfield_object() - earth).at(time)
  105. separation = sun_pos.separation_from(aster_pos)
  106. return separation.degrees
  107. get_elongation.rough_period = 1.0
  108. events = []
  109. for aster in ASTERS:
  110. if aster.skyfield_name not in ["MERCURY", "VENUS"]:
  111. continue
  112. times, elongations = find_maxima(
  113. start_time, end_time, f=get_elongation, epsilon=1.0 / 24 / 3600, num=12
  114. )
  115. for i, time in enumerate(times):
  116. elongation = elongations[i]
  117. events.append(
  118. Event(
  119. EventType.MAXIMAL_ELONGATION,
  120. [aster],
  121. translate_to_timezone(time.utc_datetime(), timezone),
  122. details="{:.3n}°".format(elongation),
  123. )
  124. )
  125. return events
  126. def _get_moon_distance():
  127. earth = get_skf_objects()["earth"]
  128. moon = get_skf_objects()["moon"]
  129. def get_distance(time: Time):
  130. earth_pos = earth.at(time)
  131. moon_pos = earth_pos.observe(moon).apparent()
  132. return moon_pos.distance().au
  133. get_distance.rough_period = 1.0
  134. return get_distance
  135. def _search_moon_apogee(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  136. moon = ASTERS[1]
  137. events = []
  138. times, _ = find_maxima(
  139. start_time, end_time, f=_get_moon_distance(), epsilon=1.0 / 24 / 60
  140. )
  141. for time in times:
  142. events.append(
  143. Event(
  144. EventType.MOON_APOGEE,
  145. [moon],
  146. translate_to_timezone(time.utc_datetime(), timezone),
  147. )
  148. )
  149. return events
  150. def _search_moon_perigee(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  151. moon = ASTERS[1]
  152. events = []
  153. times, _ = find_minima(
  154. start_time, end_time, f=_get_moon_distance(), epsilon=1.0 / 24 / 60
  155. )
  156. for time in times:
  157. events.append(
  158. Event(
  159. EventType.MOON_PERIGEE,
  160. [moon],
  161. translate_to_timezone(time.utc_datetime(), timezone),
  162. )
  163. )
  164. return events
  165. def get_events(date: date_type = date_type.today(), timezone: int = 0) -> [Event]:
  166. """Calculate and return a list of events for the given date, adjusted to the given timezone if any.
  167. Find events that happen on April 4th, 2020 (show hours in UTC):
  168. >>> get_events(date_type(2020, 4, 4))
  169. [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-04 01:14:39.063308+00:00 end=None details=None />]
  170. Find events that happen on April 4th, 2020 (show timezones in UTC+2):
  171. >>> get_events(date_type(2020, 4, 4), 2)
  172. [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-04 03:14:39.063267+02:00 end=None details=None />]
  173. Find events that happen on April 3rd, 2020 (show timezones in UTC-2):
  174. >>> get_events(date_type(2020, 4, 3), -2)
  175. [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-03 23:14:39.063388-02:00 end=None details=None />]
  176. :param date: the date for which the events must be calculated
  177. :param timezone: the timezone to adapt the results to. If not given, defaults to 0.
  178. :return: a list of events found for the given date.
  179. """
  180. start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)
  181. end_time = get_timescale().utc(date.year, date.month, date.day + 1, -timezone)
  182. try:
  183. found_events = []
  184. for fun in [
  185. _search_oppositions,
  186. _search_conjunction,
  187. _search_maximal_elongations,
  188. _search_moon_apogee,
  189. _search_moon_perigee,
  190. ]:
  191. found_events.append(fun(start_time, end_time, timezone))
  192. return sorted(flatten_list(found_events), key=lambda event: event.start_time)
  193. except EphemerisRangeError as error:
  194. start_date = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  195. end_date = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  196. start_date = date_type(start_date.year, start_date.month, start_date.day)
  197. end_date = date_type(end_date.year, end_date.month, end_date.day)
  198. raise OutOfRangeDateError(start_date, end_date) from error