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.
 
 

219 lines
10 KiB

  1. #!/usr/bin/env python3
  2. import datetime
  3. from typing import Union
  4. from skyfield.searchlib import find_discrete, find_maxima
  5. from skyfield.timelib import Time
  6. from skyfield.constants import tau
  7. from skyfield.errors import EphemerisRangeError
  8. from .model import Position, AsterEphemerides, MoonPhase, Object, ASTERS
  9. from .dateutil import translate_to_timezone, normalize_datetime
  10. from .core import get_skf_objects, get_timescale, get_iau2000b
  11. from .enum import MoonPhaseType
  12. from .exceptions import OutOfRangeDateError
  13. RISEN_ANGLE = -0.8333
  14. def _get_skyfield_to_moon_phase(
  15. times: [Time], vals: [int], now: Time, timezone: int
  16. ) -> Union[MoonPhase, None]:
  17. tomorrow = get_timescale().utc(
  18. now.utc_datetime().year, now.utc_datetime().month, now.utc_datetime().day + 1
  19. )
  20. phases = list(MoonPhaseType)
  21. current_phase = None
  22. current_phase_time = None
  23. next_phase_time = None
  24. i = 0
  25. if len(times) == 0:
  26. return None
  27. for i, time in enumerate(times):
  28. if now.utc_iso() <= time.utc_iso():
  29. if vals[i] in [0, 2, 4, 6]:
  30. if time.utc_datetime() < tomorrow.utc_datetime():
  31. current_phase_time = time
  32. current_phase = phases[vals[i]]
  33. else:
  34. i -= 1
  35. current_phase_time = None
  36. current_phase = phases[vals[i]]
  37. else:
  38. current_phase = phases[vals[i]]
  39. break
  40. for j in range(i + 1, len(times)):
  41. if vals[j] in [0, 2, 4, 6]:
  42. next_phase_time = times[j]
  43. break
  44. return MoonPhase(
  45. current_phase,
  46. translate_to_timezone(current_phase_time.utc_datetime(), timezone)
  47. if current_phase_time is not None
  48. else None,
  49. translate_to_timezone(next_phase_time.utc_datetime(), timezone)
  50. if next_phase_time is not None
  51. else None,
  52. )
  53. def get_moon_phase(
  54. for_date: datetime.date = datetime.date.today(), timezone: int = 0
  55. ) -> MoonPhase:
  56. """Calculate and return the moon phase for the given date, adjusted to the given timezone if any.
  57. Get the moon phase for the 27 March, 2021:
  58. >>> get_moon_phase(datetime.date.fromisoformat("2021-03-27"))
  59. <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 18:48:10.902298+00:00>
  60. Get the moon phase for the 27 March, 2021, in the UTC+2 timezone:
  61. >>> get_moon_phase(datetime.date.fromisoformat("2021-03-27"), timezone=2)
  62. <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 20:48:10.902298+02:00>
  63. """
  64. earth = get_skf_objects()["earth"]
  65. moon = get_skf_objects()["moon"]
  66. sun = get_skf_objects()["sun"]
  67. def moon_phase_at(time: Time):
  68. time._nutation_angles = get_iau2000b(time)
  69. current_earth = earth.at(time)
  70. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon("date")
  71. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon("date")
  72. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  73. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  74. today = get_timescale().utc(for_date.year, for_date.month, for_date.day)
  75. time1 = get_timescale().utc(for_date.year, for_date.month, for_date.day - 10)
  76. time2 = get_timescale().utc(for_date.year, for_date.month, for_date.day + 10)
  77. try:
  78. times, phase = find_discrete(time1, time2, moon_phase_at)
  79. except EphemerisRangeError as error:
  80. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  81. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  82. start = datetime.date(start.year, start.month, start.day) + datetime.timedelta(
  83. days=12
  84. )
  85. end = datetime.date(end.year, end.month, end.day) - datetime.timedelta(days=12)
  86. raise OutOfRangeDateError(start, end) from error
  87. return _get_skyfield_to_moon_phase(times, phase, today, timezone)
  88. def get_ephemerides(
  89. position: Position, date: datetime.date = datetime.date.today(), timezone: int = 0
  90. ) -> [AsterEphemerides]:
  91. """Compute and return the ephemerides for the given position and date, adjusted to the given timezone if any.
  92. Compute the ephemerides for June 9th, 2021:
  93. >>> pos = Position(50.5824, 3.0624)
  94. >>> get_ephemerides(pos, datetime.date(2021, 6, 9))
  95. [<AsterEphemerides rise_time=2021-06-09 03:36:00 culmination_time=2021-06-09 11:47:00 set_time=2021-06-09 19:58:00 aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-09 02:59:00 culmination_time=2021-06-09 11:02:00 set_time=2021-06-09 19:16:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-09 04:06:00 culmination_time=2021-06-09 11:58:00 set_time=2021-06-09 19:49:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-06-09 04:52:00 culmination_time=2021-06-09 13:13:00 set_time=2021-06-09 21:34:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-09 06:38:00 culmination_time=2021-06-09 14:40:00 set_time=2021-06-09 22:41:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-09 23:43:00 culmination_time=2021-06-09 04:54:00 set_time=2021-06-09 10:01:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-09 23:02:00 culmination_time=2021-06-09 03:41:00 set_time=2021-06-09 08:16:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-09 01:56:00 culmination_time=2021-06-09 09:18:00 set_time=2021-06-09 16:40:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-09 00:27:00 culmination_time=2021-06-09 06:13:00 set_time=2021-06-09 11:59:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=2021-06-09 22:22:00 culmination_time=2021-06-09 02:32:00 set_time=2021-06-09 06:38:00 aster=<Object type=PLANET name=PLUTO />>]
  96. Compute the ephemerides for June 9th, 2021:
  97. >>> get_ephemerides(pos, datetime.date(2021, 6, 9), timezone=2)
  98. [<AsterEphemerides rise_time=2021-06-09 05:36:00 culmination_time=2021-06-09 13:47:00 set_time=2021-06-09 21:58:00 aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-09 04:59:00 culmination_time=2021-06-09 13:02:00 set_time=2021-06-09 21:16:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-09 06:06:00 culmination_time=2021-06-09 13:58:00 set_time=2021-06-09 21:49:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-06-09 06:52:00 culmination_time=2021-06-09 15:13:00 set_time=2021-06-09 23:34:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-09 08:38:00 culmination_time=2021-06-09 16:40:00 set_time=2021-06-09 00:44:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-09 01:47:00 culmination_time=2021-06-09 06:54:00 set_time=2021-06-09 12:01:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-09 01:06:00 culmination_time=2021-06-09 05:41:00 set_time=2021-06-09 10:16:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-09 03:56:00 culmination_time=2021-06-09 11:18:00 set_time=2021-06-09 18:40:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-09 02:27:00 culmination_time=2021-06-09 08:13:00 set_time=2021-06-09 13:59:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=2021-06-09 00:26:00 culmination_time=2021-06-09 04:32:00 set_time=2021-06-09 08:38:00 aster=<Object type=PLANET name=PLUTO />>]
  99. """
  100. ephemerides = []
  101. def get_angle(for_aster: Object):
  102. def fun(time: Time) -> float:
  103. return (
  104. position.get_planet_topos()
  105. .at(time)
  106. .observe(for_aster.get_skyfield_object())
  107. .apparent()
  108. .altaz()[0]
  109. .degrees
  110. )
  111. fun.rough_period = 1.0
  112. return fun
  113. def is_risen(for_aster: Object):
  114. def fun(time: Time) -> bool:
  115. return get_angle(for_aster)(time) > RISEN_ANGLE
  116. fun.rough_period = 0.5
  117. return fun
  118. start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)
  119. end_time = get_timescale().utc(
  120. date.year, date.month, date.day, 23 - timezone, 59, 59
  121. )
  122. try:
  123. for aster in ASTERS:
  124. rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))
  125. try:
  126. culmination_time, _ = find_maxima(
  127. start_time,
  128. end_time,
  129. f=get_angle(aster),
  130. epsilon=1.0 / 3600 / 24,
  131. num=12,
  132. )
  133. culmination_time = (
  134. culmination_time[0] if len(culmination_time) > 0 else None
  135. )
  136. except ValueError:
  137. culmination_time = None
  138. if len(rise_times) == 2:
  139. rise_time = rise_times[0 if arr[0] else 1]
  140. set_time = rise_times[1 if not arr[1] else 0]
  141. else:
  142. rise_time = rise_times[0] if arr[0] else None
  143. set_time = rise_times[0] if not arr[0] else None
  144. # Convert the Time instances to Python datetime objects
  145. if rise_time is not None:
  146. rise_time = normalize_datetime(
  147. translate_to_timezone(
  148. rise_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  149. )
  150. )
  151. if culmination_time is not None:
  152. culmination_time = normalize_datetime(
  153. translate_to_timezone(
  154. culmination_time.utc_datetime().replace(microsecond=0),
  155. to_tz=timezone,
  156. )
  157. )
  158. if set_time is not None:
  159. set_time = normalize_datetime(
  160. translate_to_timezone(
  161. set_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  162. )
  163. )
  164. ephemerides.append(
  165. AsterEphemerides(rise_time, culmination_time, set_time, aster=aster)
  166. )
  167. except EphemerisRangeError as error:
  168. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  169. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  170. start = datetime.date(start.year, start.month, start.day + 1)
  171. end = datetime.date(end.year, end.month, end.day - 1)
  172. raise OutOfRangeDateError(start, end) from error
  173. return ephemerides