A library that computes the ephemerides.
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

531 lignes
19 KiB

  1. #!/usr/bin/env python3
  2. # Kosmorrolib - The Library To Compute Your Ephemerides
  3. # Copyright (C) 2021 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, timedelta
  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 skyfield.units import Angle
  22. from skyfield import almanac, eclipselib
  23. from math import pi
  24. from kosmorrolib.model import (
  25. Event,
  26. Object,
  27. Star,
  28. Planet,
  29. get_aster,
  30. ASTERS,
  31. EARTH,
  32. )
  33. from kosmorrolib.dateutil import translate_to_timezone
  34. from kosmorrolib.enum import EventType, ObjectIdentifier, SeasonType, LunarEclipseType
  35. from kosmorrolib.exceptions import OutOfRangeDateError
  36. from kosmorrolib.core import get_timescale, get_skf_objects, flatten_list
  37. def _search_conjunctions_occultations(
  38. start_time: Time, end_time: Time, timezone: int
  39. ) -> [Event]:
  40. """Function to search conjunction.
  41. **Warning:** this is an internal function, not intended for use by end-developers.
  42. Will return MOON and VENUS opposition on 2021-06-12:
  43. >>> conjunction = _search_conjunctions_occultations(get_timescale().utc(2021, 6, 12), get_timescale().utc(2021, 6, 13), 0)
  44. >>> len(conjunction)
  45. 1
  46. >>> conjunction[0].objects
  47. [<Object type=SATELLITE name=MOON />, <Object type=PLANET name=VENUS />]
  48. Will return nothing if no conjunction happens:
  49. >>> _search_conjunctions_occultations(get_timescale().utc(2021, 6, 17),get_timescale().utc(2021, 6, 18), 0)
  50. []
  51. This function detects occultations too:
  52. >>> _search_conjunctions_occultations(get_timescale().utc(2021, 4, 17),get_timescale().utc(2021, 4, 18), 0)
  53. [<Event type=OCCULTATION objects=[<Object type=SATELLITE name=MOON />, <Object type=PLANET name=MARS />] start=2021-04-17 12:08:16.115650+00:00 end=None details=None />]
  54. """
  55. earth = get_skf_objects()["earth"]
  56. aster1 = None
  57. aster2 = None
  58. def is_in_conjunction(time: Time):
  59. earth_pos = earth.at(time)
  60. _, aster1_lon, _ = (
  61. earth_pos.observe(aster1.skyfield_object).apparent().ecliptic_latlon()
  62. )
  63. _, aster2_lon, _ = (
  64. earth_pos.observe(aster2.skyfield_object).apparent().ecliptic_latlon()
  65. )
  66. return ((aster1_lon.radians - aster2_lon.radians) / pi % 2.0).astype(
  67. "int8"
  68. ) == 0
  69. is_in_conjunction.rough_period = 60.0
  70. computed = []
  71. events = []
  72. for aster1 in ASTERS:
  73. # Ignore the Sun
  74. if isinstance(aster1, Star):
  75. continue
  76. for aster2 in ASTERS:
  77. if isinstance(aster2, Star) or aster2 == aster1 or aster2 in computed:
  78. continue
  79. times, is_conjs = find_discrete(start_time, end_time, is_in_conjunction)
  80. for i, time in enumerate(times):
  81. if is_conjs[i]:
  82. aster1_pos = (aster1.skyfield_object - earth).at(time)
  83. aster2_pos = (aster2.skyfield_object - earth).at(time)
  84. distance = aster1_pos.separation_from(aster2_pos).degrees
  85. if (
  86. distance - aster2.get_apparent_radius(time).degrees
  87. < aster1.get_apparent_radius(time).degrees
  88. ):
  89. occulting_aster = (
  90. [aster1, aster2]
  91. if aster1_pos.distance().km < aster2_pos.distance().km
  92. else [aster2, aster1]
  93. )
  94. events.append(
  95. Event(
  96. EventType.OCCULTATION,
  97. occulting_aster,
  98. translate_to_timezone(time.utc_datetime(), timezone),
  99. )
  100. )
  101. else:
  102. events.append(
  103. Event(
  104. EventType.CONJUNCTION,
  105. [aster1, aster2],
  106. translate_to_timezone(time.utc_datetime(), timezone),
  107. )
  108. )
  109. computed.append(aster1)
  110. return events
  111. def _search_oppositions(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  112. """Function to search oppositions.
  113. **Warning:** this is an internal function, not intended for use by end-developers.
  114. Will return Mars opposition on 2020-10-13:
  115. >>> oppositions = _search_oppositions(get_timescale().utc(2020, 10, 13), get_timescale().utc(2020, 10, 14), 0)
  116. >>> len(oppositions)
  117. 1
  118. >>> oppositions[0].objects[0]
  119. <Object type=PLANET name=MARS />
  120. Will return nothing if no opposition happens:
  121. >>> _search_oppositions(get_timescale().utc(2021, 3, 20), get_timescale().utc(2021, 3, 21), 0)
  122. []
  123. >>> _search_oppositions(get_timescale().utc(2022, 12, 24), get_timescale().utc(2022, 12, 25), 0)
  124. []
  125. """
  126. earth = get_skf_objects()["earth"]
  127. sun = get_skf_objects()["sun"]
  128. aster = None
  129. def is_oppositing(time: Time) -> [bool]:
  130. diff = get_angle(time)
  131. return diff > 180
  132. def get_angle(time: Time):
  133. earth_pos = earth.at(time)
  134. sun_pos = earth_pos.observe(
  135. sun
  136. ).apparent() # Never do this without eyes protection!
  137. aster_pos = earth_pos.observe(aster.skyfield_object).apparent()
  138. _, lon1, _ = sun_pos.ecliptic_latlon()
  139. _, lon2, _ = aster_pos.ecliptic_latlon()
  140. return lon1.degrees - lon2.degrees
  141. is_oppositing.rough_period = 1.0
  142. events = []
  143. for aster in ASTERS:
  144. if not isinstance(aster, Planet) or aster.identifier in [
  145. ObjectIdentifier.MERCURY,
  146. ObjectIdentifier.VENUS,
  147. ]:
  148. continue
  149. times, _ = find_discrete(start_time, end_time, is_oppositing)
  150. for time in times:
  151. if int(get_angle(time)) != 180:
  152. # If the angle is negative, then it is actually a false positive.
  153. # Just ignoring it.
  154. continue
  155. events.append(
  156. Event(
  157. EventType.OPPOSITION,
  158. [aster],
  159. translate_to_timezone(time.utc_datetime(), timezone),
  160. )
  161. )
  162. return events
  163. def _search_maximal_elongations(
  164. start_time: Time, end_time: Time, timezone: int
  165. ) -> [Event]:
  166. """Function to search oppositions.
  167. **Warning:** this is an internal function, not intended for use by end-developers.
  168. Will return Mercury maimum elogation for September 14, 2021:
  169. >>> get_events(date(2021, 9, 14))
  170. [<Event type=MAXIMAL_ELONGATION objects=[<Object type=PLANET name=MERCURY />] start=2021-09-14 04:13:46.664879+00:00 end=None details={'deg': 26.8} />]
  171. """
  172. earth = get_skf_objects()["earth"]
  173. sun = get_skf_objects()["sun"]
  174. def get_elongation(planet: Object):
  175. def f(time: Time):
  176. sun_pos = (sun - earth).at(time)
  177. aster_pos = (planet.skyfield_object - earth).at(time)
  178. separation = sun_pos.separation_from(aster_pos)
  179. return separation.degrees
  180. f.rough_period = 1.0
  181. return f
  182. events = []
  183. for identifier in [ObjectIdentifier.MERCURY, ObjectIdentifier.VENUS]:
  184. planet = get_aster(identifier)
  185. times, elongations = find_maxima(
  186. start_time,
  187. end_time,
  188. f=get_elongation(planet),
  189. epsilon=1.0 / 24 / 3600,
  190. num=12,
  191. )
  192. for i, time in enumerate(times):
  193. elongation = round(elongations[i], 1)
  194. events.append(
  195. Event(
  196. EventType.MAXIMAL_ELONGATION,
  197. [planet],
  198. translate_to_timezone(time.utc_datetime(), timezone),
  199. details={"deg": elongation},
  200. )
  201. )
  202. return events
  203. def _get_distance(to_aster: Object, from_aster: Object):
  204. def get_distance(time: Time):
  205. from_pos = from_aster.skyfield_object.at(time)
  206. to_pos = from_pos.observe(to_aster.skyfield_object)
  207. return to_pos.distance().km
  208. get_distance.rough_period = 1.0
  209. return get_distance
  210. def _search_apogee(to_aster: Object, from_aster: Object = EARTH) -> callable:
  211. """Search for moon apogee
  212. **Warning:** this is an internal function, not intended for use by end-developers.
  213. Get the moon apogee:
  214. >>> _search_apogee(ASTERS[1])(get_timescale().utc(2021, 6, 8), get_timescale().utc(2021, 6, 9), 0)
  215. [<Event type=APOGEE objects=[<Object type=SATELLITE name=MOON />] start=2021-06-08 02:39:40.165271+00:00 end=None details={'distance_km': 406211.04850197025} />]
  216. Get the Earth's apogee:
  217. >>> _search_apogee(EARTH, from_aster=ASTERS[0])(get_timescale().utc(2021, 7, 5), get_timescale().utc(2021, 7, 6), 0)
  218. [<Event type=APOGEE objects=[<Object type=PLANET name=EARTH />] start=2021-07-05 22:35:42.148792+00:00 end=None details={'distance_km': 152100521.91712126} />]
  219. """
  220. def f(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  221. events = []
  222. times, distances = find_maxima(
  223. start_time,
  224. end_time,
  225. f=_get_distance(to_aster, from_aster),
  226. epsilon=1.0 / 24 / 60,
  227. )
  228. for i, time in enumerate(times):
  229. events.append(
  230. Event(
  231. EventType.APOGEE,
  232. [to_aster],
  233. translate_to_timezone(time.utc_datetime(), timezone),
  234. details={"distance_km": distances[i]},
  235. )
  236. )
  237. return events
  238. return f
  239. def _search_perigee(aster: Object, from_aster: Object = EARTH) -> callable:
  240. """Search for moon perigee
  241. **Warning:** this is an internal function, not intended for use by end-developers.
  242. Get the moon perigee:
  243. >>> _search_perigee(ASTERS[1])(get_timescale().utc(2021, 5, 26), get_timescale().utc(2021, 5, 27), 0)
  244. [<Event type=PERIGEE objects=[<Object type=SATELLITE name=MOON />] start=2021-05-26 01:56:01.983455+00:00 end=None details={'distance_km': 357313.9680798693} />]
  245. Get the Earth's perigee:
  246. >>> _search_perigee(EARTH, from_aster=ASTERS[0])(get_timescale().utc(2021, 1, 2), get_timescale().utc(2021, 1, 3), 0)
  247. [<Event type=PERIGEE objects=[<Object type=PLANET name=EARTH />] start=2021-01-02 13:59:00.495905+00:00 end=None details={'distance_km': 147093166.1686309} />]
  248. """
  249. def f(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  250. events = []
  251. times, distances = find_minima(
  252. start_time,
  253. end_time,
  254. f=_get_distance(aster, from_aster),
  255. epsilon=1.0 / 24 / 60,
  256. )
  257. for i, time in enumerate(times):
  258. events.append(
  259. Event(
  260. EventType.PERIGEE,
  261. [aster],
  262. translate_to_timezone(time.utc_datetime(), timezone),
  263. details={"distance_km": distances[i]},
  264. )
  265. )
  266. return events
  267. return f
  268. def _search_earth_season_change(
  269. start_time: Time, end_time: Time, timezone: int
  270. ) -> [Event]:
  271. """Function to find earth season change event.
  272. **Warning:** this is an internal function, not intended for use by end-developers.
  273. Will return JUNE SOLSTICE on 2020/06/20:
  274. >>> season_change = _search_earth_season_change(get_timescale().utc(2020, 6, 20), get_timescale().utc(2020, 6, 21), 0)
  275. >>> len(season_change)
  276. 1
  277. >>> season_change[0].event_type
  278. <EventType.SEASON_CHANGE: 7>
  279. >>> season_change[0].details
  280. {'season': <SeasonType.JUNE_SOLSTICE: 1>}
  281. Will return nothing if there is no season change event in the period of time being calculated:
  282. >>> _search_earth_season_change(get_timescale().utc(2021, 6, 17), get_timescale().utc(2021, 6, 18), 0)
  283. []
  284. """
  285. events = []
  286. event_time, event_id = almanac.find_discrete(
  287. start_time, end_time, almanac.seasons(get_skf_objects())
  288. )
  289. if len(event_time) == 0:
  290. return []
  291. events.append(
  292. Event(
  293. EventType.SEASON_CHANGE,
  294. [],
  295. translate_to_timezone(event_time.utc_datetime()[0], timezone),
  296. details={"season": SeasonType(event_id[0])},
  297. )
  298. )
  299. return events
  300. def _search_lunar_eclipse(start_time: Time, end_time: Time, timezone: int) -> [Event]:
  301. """Function to detect lunar eclipses.
  302. **Warning:** this is an internal function, not intended for use by end-developers.
  303. Will return a total lunar eclipse for 2021-05-26:
  304. >>> _search_lunar_eclipse(get_timescale().utc(2021, 5, 26), get_timescale().utc(2021, 5, 27), 0)
  305. [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2021-05-26 08:47:54.795821+00:00 end=2021-05-26 13:49:34.353411+00:00 details={'type': <LunarEclipseType.TOTAL: 2>, 'maximum': datetime.datetime(2021, 5, 26, 11, 18, 42, 328842, tzinfo=datetime.timezone.utc)} />]
  306. >>> _search_lunar_eclipse(get_timescale().utc(2019, 7, 16), get_timescale().utc(2019, 7, 17), 0)
  307. [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2019-07-16 18:39:53.391337+00:00 end=2019-07-17 00:21:51.378940+00:00 details={'type': <LunarEclipseType.PARTIAL: 1>, 'maximum': datetime.datetime(2019, 7, 16, 21, 30, 44, 170096, tzinfo=datetime.timezone.utc)} />]
  308. >>> _search_lunar_eclipse(get_timescale().utc(2017, 2, 11), get_timescale().utc(2017, 2, 12), 0)
  309. [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2017-02-10 22:02:59.016572+00:00 end=2017-02-11 03:25:07.627886+00:00 details={'type': <LunarEclipseType.PENUMBRAL: 0>, 'maximum': datetime.datetime(2017, 2, 11, 0, 43, 51, 793786, tzinfo=datetime.timezone.utc)} />]
  310. """
  311. moon = get_aster(ObjectIdentifier.MOON)
  312. events = []
  313. t, y, details = eclipselib.lunar_eclipses(start_time, end_time, get_skf_objects())
  314. for ti, yi in zip(t, y):
  315. penumbra_radius = Angle(radians=details["penumbra_radius_radians"][0])
  316. _, max_lon, _ = (
  317. EARTH.skyfield_object.at(ti)
  318. .observe(moon.skyfield_object)
  319. .apparent()
  320. .ecliptic_latlon()
  321. )
  322. def is_in_penumbra(time: Time):
  323. _, lon, _ = (
  324. EARTH.skyfield_object.at(time)
  325. .observe(moon.skyfield_object)
  326. .apparent()
  327. .ecliptic_latlon()
  328. )
  329. moon_radius = details["moon_radius_radians"]
  330. return (
  331. abs(max_lon.radians - lon.radians)
  332. < penumbra_radius.radians + moon_radius
  333. )
  334. is_in_penumbra.rough_period = 60.0
  335. search_start_time = get_timescale().from_datetime(
  336. start_time.utc_datetime() - timedelta(days=1)
  337. )
  338. search_end_time = get_timescale().from_datetime(
  339. end_time.utc_datetime() + timedelta(days=1)
  340. )
  341. eclipse_start, _ = find_discrete(search_start_time, ti, is_in_penumbra)
  342. eclipse_end, _ = find_discrete(ti, search_end_time, is_in_penumbra)
  343. events.append(
  344. Event(
  345. EventType.LUNAR_ECLIPSE,
  346. [moon],
  347. start_time=translate_to_timezone(
  348. eclipse_start[0].utc_datetime(), timezone
  349. ),
  350. end_time=translate_to_timezone(eclipse_end[0].utc_datetime(), timezone),
  351. details={
  352. "type": LunarEclipseType(yi),
  353. "maximum": translate_to_timezone(ti.utc_datetime(), timezone),
  354. },
  355. )
  356. )
  357. return events
  358. def get_events(for_date: date = date.today(), timezone: int = 0) -> [Event]:
  359. """Calculate and return a list of events for the given date, adjusted to the given timezone if any.
  360. Find events that happen on April 4th, 2020 (show hours in UTC):
  361. >>> get_events(date(2020, 4, 4))
  362. [<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 />]
  363. Find events that happen on April 4th, 2020 (show timezones in UTC+2):
  364. >>> get_events(date(2020, 4, 4), 2)
  365. [<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 />]
  366. Find events that happen on April 3rd, 2020 (show timezones in UTC-2):
  367. >>> get_events(date(2020, 4, 3), -2)
  368. [<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 />]
  369. If there is no events for the given date, then an empty list is returned:
  370. >>> get_events(date(2021, 4, 20))
  371. []
  372. Note that the events can only be found for a date range.
  373. Asking for the events with an out of range date will result in an exception:
  374. >>> get_events(date(1000, 1, 1))
  375. Traceback (most recent call last):
  376. ...
  377. kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-07-28 and 2053-10-08
  378. :param for_date: the date for which the events must be calculated
  379. :param timezone: the timezone to adapt the results to. If not given, defaults to 0.
  380. :return: a list of events found for the given date.
  381. """
  382. start_time = get_timescale().utc(
  383. for_date.year, for_date.month, for_date.day, -timezone
  384. )
  385. end_time = get_timescale().utc(
  386. for_date.year, for_date.month, for_date.day + 1, -timezone
  387. )
  388. try:
  389. found_events = []
  390. for fun in [
  391. _search_oppositions,
  392. _search_conjunctions_occultations,
  393. _search_maximal_elongations,
  394. _search_apogee(ASTERS[1]),
  395. _search_perigee(ASTERS[1]),
  396. _search_apogee(EARTH, from_aster=ASTERS[0]),
  397. _search_perigee(EARTH, from_aster=ASTERS[0]),
  398. _search_earth_season_change,
  399. _search_lunar_eclipse,
  400. ]:
  401. found_events.append(fun(start_time, end_time, timezone))
  402. return sorted(flatten_list(found_events), key=lambda event: event.start_time)
  403. except EphemerisRangeError as error:
  404. start_date = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  405. end_date = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  406. start_date = date(start_date.year, start_date.month, start_date.day)
  407. end_date = date(end_date.year, end_date.month, end_date.day)
  408. raise OutOfRangeDateError(start_date, end_date) from error