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.

ephemerides.py 18 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276
  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, datetime, timedelta
  18. from typing import Union
  19. from skyfield.searchlib import find_discrete, find_maxima
  20. from skyfield.timelib import Time
  21. from skyfield.constants import tau
  22. from skyfield.errors import EphemerisRangeError
  23. from .model import Position, AsterEphemerides, MoonPhase, Object, ASTERS
  24. from .dateutil import translate_to_timezone, normalize_datetime
  25. from .core import get_skf_objects, get_timescale, get_iau2000b
  26. from .enum import MoonPhaseType
  27. from .exceptions import OutOfRangeDateError
  28. RISEN_ANGLE = -0.8333
  29. def _get_skyfield_to_moon_phase(
  30. times: [Time], vals: [int], now: Time, timezone: int
  31. ) -> Union[MoonPhase, None]:
  32. tomorrow = get_timescale().utc(
  33. now.utc_datetime().year, now.utc_datetime().month, now.utc_datetime().day + 1
  34. )
  35. next_phase_time = None
  36. i = 0
  37. # Find the current moon phase:
  38. for i, time in enumerate(times):
  39. if now.utc_datetime() <= time.utc_datetime():
  40. if time.utc_datetime() >= tomorrow.utc_datetime():
  41. i -= 1
  42. break
  43. current_phase = MoonPhaseType(vals[i])
  44. if current_phase in [
  45. MoonPhaseType.NEW_MOON,
  46. MoonPhaseType.FIRST_QUARTER,
  47. MoonPhaseType.FULL_MOON,
  48. MoonPhaseType.LAST_QUARTER,
  49. ]:
  50. current_phase_time = translate_to_timezone(times[i].utc_datetime(), timezone)
  51. else:
  52. current_phase_time = None
  53. # Find the next moon phase
  54. for j in range(i + 1, len(times)):
  55. if vals[j] in [0, 2, 4, 6]:
  56. next_phase_time = translate_to_timezone(times[j].utc_datetime(), timezone)
  57. break
  58. return MoonPhase(current_phase, current_phase_time, next_phase_time)
  59. def get_moon_phase(for_date: date = date.today(), timezone: int = 0) -> MoonPhase:
  60. """Calculate and return the moon phase for the given date, adjusted to the given timezone if any.
  61. Get the moon phase for the 27 March, 2021:
  62. >>> get_moon_phase(date(2021, 3, 27))
  63. <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 18:48:10.902298+00:00>
  64. When the moon phase is a new moon, a first quarter, a full moon or a last quarter, you get the exact time
  65. of its happening too:
  66. >>> get_moon_phase(datetime(2021, 3, 28))
  67. <MoonPhase phase_type=MoonPhaseType.FULL_MOON time=2021-03-28 18:48:10.902298+00:00 next_phase_date=2021-04-04 10:02:27.393689+00:00>
  68. Get the moon phase for the 27 March, 2021, in the UTC+2 timezone:
  69. >>> get_moon_phase(date(2021, 3, 27), timezone=2)
  70. <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 20:48:10.902298+02:00>
  71. Note that the moon phase can only be computed for a date range.
  72. Asking for the moon phase with an out of range date will result in an exception:
  73. >>> get_moon_phase(date(1000, 1, 1))
  74. Traceback (most recent call last):
  75. ...
  76. kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-08-09 and 2053-09-26
  77. """
  78. earth = get_skf_objects()["earth"]
  79. moon = get_skf_objects()["moon"]
  80. sun = get_skf_objects()["sun"]
  81. def moon_phase_at(time: Time):
  82. time._nutation_angles = get_iau2000b(time)
  83. current_earth = earth.at(time)
  84. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon("date")
  85. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon("date")
  86. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  87. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  88. today = get_timescale().utc(for_date.year, for_date.month, for_date.day)
  89. start_time = get_timescale().utc(for_date.year, for_date.month, for_date.day - 10)
  90. end_time = get_timescale().utc(for_date.year, for_date.month, for_date.day + 10)
  91. try:
  92. times, phases = find_discrete(start_time, end_time, moon_phase_at)
  93. return _get_skyfield_to_moon_phase(times, phases, today, timezone)
  94. except EphemerisRangeError as error:
  95. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  96. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  97. start = date(start.year, start.month, start.day) + timedelta(days=12)
  98. end = date(end.year, end.month, end.day) - timedelta(days=12)
  99. raise OutOfRangeDateError(start, end) from error
  100. def get_ephemerides(
  101. position: Position, for_date: date = date.today(), timezone: int = 0
  102. ) -> [AsterEphemerides]:
  103. """Compute and return the ephemerides for the given position and date, adjusted to the given timezone if any.
  104. Compute the ephemerides for July 7th, 2022:
  105. >>> get_ephemerides(Position(36.6794, 4.8555), date(2022, 7, 7))
  106. [<AsterEphemerides rise_time=2022-07-07 04:29:00 culmination_time=2022-07-07 11:46:00 set_time=2022-07-07 19:02:00 aster=<Object type=STAR name=SUN />>,
  107. <AsterEphemerides rise_time=2022-07-07 12:16:00 culmination_time=2022-07-07 18:06:00 set_time=2022-07-07 23:54:00 aster=<Object type=SATELLITE name=MOON />>,
  108. <AsterEphemerides rise_time=2022-07-07 03:36:00 culmination_time=2022-07-07 10:58:00 set_time=2022-07-07 18:21:00 aster=<Object type=PLANET name=MERCURY />>,
  109. <AsterEphemerides rise_time=2022-07-07 02:30:00 culmination_time=2022-07-07 09:44:00 set_time=2022-07-07 16:58:00 aster=<Object type=PLANET name=VENUS />>,
  110. <AsterEphemerides rise_time=2022-07-07 00:05:00 culmination_time=2022-07-07 06:39:00 set_time=2022-07-07 13:14:00 aster=<Object type=PLANET name=MARS />>,
  111. <AsterEphemerides rise_time=2022-07-07 22:59:00 culmination_time=2022-07-07 05:11:00 set_time=2022-07-07 11:20:00 aster=<Object type=PLANET name=JUPITER />>,
  112. <AsterEphemerides rise_time=2022-07-07 21:06:00 culmination_time=2022-07-07 02:29:00 set_time=2022-07-07 07:48:00 aster=<Object type=PLANET name=SATURN />>,
  113. <AsterEphemerides rise_time=2022-07-07 00:47:00 culmination_time=2022-07-07 07:42:00 set_time=2022-07-07 14:38:00 aster=<Object type=PLANET name=URANUS />>,
  114. <AsterEphemerides rise_time=2022-07-07 22:27:00 culmination_time=2022-07-07 04:25:00 set_time=2022-07-07 10:20:00 aster=<Object type=PLANET name=NEPTUNE />>,
  115. <AsterEphemerides rise_time=2022-07-07 19:46:00 culmination_time=2022-07-07 00:41:00 set_time=2022-07-07 05:33:00 aster=<Object type=PLANET name=PLUTO />>]
  116. Timezone can be optionnaly set to adapt the hours to your location:
  117. >>> get_ephemerides(Position(36.6794, 4.8555), date(2022, 7, 7), timezone=2)
  118. [<AsterEphemerides rise_time=2022-07-07 06:29:00 culmination_time=2022-07-07 13:46:00 set_time=2022-07-07 21:02:00 aster=<Object type=STAR name=SUN />>,
  119. <AsterEphemerides rise_time=2022-07-07 14:16:00 culmination_time=2022-07-07 20:06:00 set_time=2022-07-08 01:54:00 aster=<Object type=SATELLITE name=MOON />>,
  120. <AsterEphemerides rise_time=2022-07-07 05:36:00 culmination_time=2022-07-07 12:58:00 set_time=2022-07-07 20:21:00 aster=<Object type=PLANET name=MERCURY />>,
  121. <AsterEphemerides rise_time=2022-07-07 04:30:00 culmination_time=2022-07-07 11:44:00 set_time=2022-07-07 18:58:00 aster=<Object type=PLANET name=VENUS />>,
  122. <AsterEphemerides rise_time=2022-07-08 02:03:00 culmination_time=2022-07-07 08:39:00 set_time=2022-07-07 15:14:00 aster=<Object type=PLANET name=MARS />>,
  123. <AsterEphemerides rise_time=2022-07-08 00:59:00 culmination_time=2022-07-07 07:11:00 set_time=2022-07-07 13:20:00 aster=<Object type=PLANET name=JUPITER />>,
  124. <AsterEphemerides rise_time=2022-07-07 23:06:00 culmination_time=2022-07-07 04:29:00 set_time=2022-07-07 09:48:00 aster=<Object type=PLANET name=SATURN />>,
  125. <AsterEphemerides rise_time=2022-07-08 02:43:00 culmination_time=2022-07-07 09:42:00 set_time=2022-07-07 16:38:00 aster=<Object type=PLANET name=URANUS />>,
  126. <AsterEphemerides rise_time=2022-07-08 00:27:00 culmination_time=2022-07-07 06:25:00 set_time=2022-07-07 12:20:00 aster=<Object type=PLANET name=NEPTUNE />>,
  127. <AsterEphemerides rise_time=2022-07-07 21:46:00 culmination_time=2022-07-08 02:37:00 set_time=2022-07-07 07:33:00 aster=<Object type=PLANET name=PLUTO />>]
  128. Objects may not rise or set on the given date (e.g. they rise the previous day or set the next day).
  129. In this case, you will get `None` values on the rise or set time.
  130. If an objet does not rise nor set due to your latitude, then both rise and set will be `None`:
  131. >>> north_pole = Position(70, 20)
  132. >>> south_pole = Position(-70, 20)
  133. >>> get_ephemerides(north_pole, date(2021, 6, 20))
  134. [<AsterEphemerides rise_time=None culmination_time=2021-06-20 10:42:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-20 14:30:00 culmination_time=2021-06-20 18:44:00 set_time=2021-06-20 22:53:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-20 22:56:00 culmination_time=2021-06-20 09:47:00 set_time=2021-06-20 20:34:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 12:20:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 13:17:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-20 23:06:00 culmination_time=2021-06-20 03:04:00 set_time=2021-06-20 06:58:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-20 23:28:00 culmination_time=2021-06-20 01:48:00 set_time=2021-06-20 04:05:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-20 21:53:00 culmination_time=2021-06-20 07:29:00 set_time=2021-06-20 17:02:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-20 22:51:00 culmination_time=2021-06-20 04:22:00 set_time=2021-06-20 09:50:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 00:40:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]
  135. >>> get_ephemerides(north_pole, date(2021, 12, 21))
  136. [<AsterEphemerides rise_time=None culmination_time=2021-12-21 10:38:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 00:04:00 set_time=None aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 11:33:00 set_time=None aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-12-21 11:58:00 culmination_time=2021-12-21 12:33:00 set_time=2021-12-21 13:08:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 08:54:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-12-21 11:07:00 culmination_time=2021-12-21 14:43:00 set_time=2021-12-21 18:19:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-12-21 11:32:00 culmination_time=2021-12-21 13:33:00 set_time=2021-12-21 15:33:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-12-21 09:54:00 culmination_time=2021-12-21 19:13:00 set_time=2021-12-21 04:37:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-12-21 10:49:00 culmination_time=2021-12-21 16:05:00 set_time=2021-12-21 21:21:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 12:31:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]
  137. >>> get_ephemerides(south_pole, date(2021, 6, 20))
  138. [<AsterEphemerides rise_time=None culmination_time=2021-06-20 10:42:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-20 11:10:00 culmination_time=2021-06-20 19:06:00 set_time=2021-06-20 01:20:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-20 07:47:00 culmination_time=2021-06-20 09:47:00 set_time=2021-06-20 11:48:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 12:20:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-20 12:14:00 culmination_time=2021-06-20 13:17:00 set_time=2021-06-20 14:21:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-20 18:32:00 culmination_time=2021-06-20 03:04:00 set_time=2021-06-20 11:32:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-20 15:20:00 culmination_time=2021-06-20 01:48:00 set_time=2021-06-20 12:12:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-20 04:32:00 culmination_time=2021-06-20 07:29:00 set_time=2021-06-20 10:26:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-20 21:28:00 culmination_time=2021-06-20 04:22:00 set_time=2021-06-20 11:13:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 00:40:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]
  139. >>> get_ephemerides(south_pole, date(2021, 12, 22))
  140. [<AsterEphemerides rise_time=None culmination_time=2021-12-22 10:39:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 01:01:00 set_time=None aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 11:35:00 set_time=None aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 12:27:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 08:53:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-12-22 05:52:00 culmination_time=2021-12-22 14:40:00 set_time=2021-12-22 23:26:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-12-22 02:41:00 culmination_time=2021-12-22 13:29:00 set_time=2021-12-22 00:21:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-12-22 16:01:00 culmination_time=2021-12-22 19:09:00 set_time=2021-12-22 22:17:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-12-22 08:59:00 culmination_time=2021-12-22 16:01:00 set_time=2021-12-22 23:04:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 12:27:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]
  141. Note that the ephemerides can only be computed for a date range.
  142. Asking for the ephemerides with an out of range date will result in an exception:
  143. >>> get_ephemerides(Position(50.5824, 3.0624), date(1000, 1, 1))
  144. Traceback (most recent call last):
  145. ...
  146. kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-07-29 and 2053-10-07
  147. """
  148. ephemerides = []
  149. def get_angle(for_aster: Object):
  150. def fun(time: Time) -> float:
  151. return (
  152. position.get_planet_topos()
  153. .at(time)
  154. .observe(for_aster.skyfield_object)
  155. .apparent()
  156. .altaz()[0]
  157. .degrees
  158. )
  159. fun.rough_period = 1.0
  160. return fun
  161. def is_risen(for_aster: Object):
  162. def fun(time: Time) -> bool:
  163. return get_angle(for_aster)(time) > RISEN_ANGLE
  164. fun.rough_period = 0.5
  165. return fun
  166. start_time = get_timescale().utc(
  167. for_date.year, for_date.month, for_date.day, timezone
  168. )
  169. end_time = get_timescale().utc(
  170. for_date.year, for_date.month, for_date.day, 23 + timezone, 59, 59
  171. )
  172. try:
  173. for aster in ASTERS:
  174. times, risen_info = find_discrete(start_time, end_time, is_risen(aster))
  175. culmination_time, _ = find_maxima(start_time, end_time, get_angle(aster))
  176. rise_time, set_time = None, None
  177. culmination_time = (
  178. culmination_time[0] if len(culmination_time) == 1 else None
  179. )
  180. if len(times) > 0:
  181. rise_time = times[0] if risen_info[0] else None
  182. set_time = times[0] if not risen_info[0] else None
  183. if len(times) == 2:
  184. rise_time = times[0 if risen_info[0] else 1]
  185. set_time = times[1 if not risen_info[1] else 0]
  186. # Convert the Time instances to Python datetime objects
  187. if rise_time is not None:
  188. rise_time = normalize_datetime(
  189. translate_to_timezone(
  190. rise_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  191. )
  192. )
  193. if culmination_time is not None:
  194. culmination_time = normalize_datetime(
  195. translate_to_timezone(
  196. culmination_time.utc_datetime().replace(microsecond=0),
  197. to_tz=timezone,
  198. )
  199. )
  200. if set_time is not None:
  201. set_time = normalize_datetime(
  202. translate_to_timezone(
  203. set_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  204. )
  205. )
  206. ephemerides.append(
  207. AsterEphemerides(rise_time, culmination_time, set_time, aster=aster)
  208. )
  209. except EphemerisRangeError as error:
  210. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  211. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  212. start = date(start.year, start.month, start.day + 1)
  213. end = date(end.year, end.month, end.day - 1)
  214. raise OutOfRangeDateError(start, end) from error
  215. return ephemerides