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

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