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.

ephemerides.py 6.9 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  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 .data import Position, AsterEphemerides, MoonPhase, Object, ASTERS
  9. from .dateutil import translate_to_timezone
  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. ephemerides = []
  92. def get_angle(for_aster: Object):
  93. def fun(time: Time) -> float:
  94. return (
  95. position.get_planet_topos()
  96. .at(time)
  97. .observe(for_aster.get_skyfield_object())
  98. .apparent()
  99. .altaz()[0]
  100. .degrees
  101. )
  102. fun.rough_period = 1.0
  103. return fun
  104. def is_risen(for_aster: Object):
  105. def fun(time: Time) -> bool:
  106. return get_angle(for_aster)(time) > RISEN_ANGLE
  107. fun.rough_period = 0.5
  108. return fun
  109. start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)
  110. end_time = get_timescale().utc(
  111. date.year, date.month, date.day, 23 - timezone, 59, 59
  112. )
  113. try:
  114. for aster in ASTERS:
  115. rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))
  116. try:
  117. culmination_time, _ = find_maxima(
  118. start_time,
  119. end_time,
  120. f=get_angle(aster),
  121. epsilon=1.0 / 3600 / 24,
  122. num=12,
  123. )
  124. culmination_time = (
  125. culmination_time[0] if len(culmination_time) > 0 else None
  126. )
  127. except ValueError:
  128. culmination_time = None
  129. if len(rise_times) == 2:
  130. rise_time = rise_times[0 if arr[0] else 1]
  131. set_time = rise_times[1 if not arr[1] else 0]
  132. else:
  133. rise_time = rise_times[0] if arr[0] else None
  134. set_time = rise_times[0] if not arr[0] else None
  135. # Convert the Time instances to Python datetime objects
  136. if rise_time is not None:
  137. rise_time = translate_to_timezone(
  138. rise_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  139. )
  140. if culmination_time is not None:
  141. culmination_time = translate_to_timezone(
  142. culmination_time.utc_datetime().replace(microsecond=0),
  143. to_tz=timezone,
  144. )
  145. if set_time is not None:
  146. set_time = translate_to_timezone(
  147. set_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  148. )
  149. ephemerides.append(
  150. AsterEphemerides(rise_time, culmination_time, set_time, aster=aster)
  151. )
  152. except EphemerisRangeError as error:
  153. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  154. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  155. start = datetime.date(start.year, start.month, start.day + 1)
  156. end = datetime.date(end.year, end.month, end.day - 1)
  157. raise OutOfRangeDateError(start, end) from error
  158. return ephemerides