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 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. >>> _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:49:13.314888+00:00 end=2021-05-26 13:48:15.757096+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:41:24.400419+00:00 end=2019-07-17 00:20:20.079536+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:04:24.192412+00:00 end=2017-02-11 03:23:42.046415+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