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.
 
 

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