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.
 
 

188 lines
6.2 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 .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
  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. current_phase_time.utc_datetime() if current_phase_time is not None else None,
  47. next_phase_time.utc_datetime() if next_phase_time is not None else None,
  48. )
  49. def get_moon_phase(
  50. compute_date: datetime.date = datetime.date.today(), timezone: int = 0
  51. ) -> MoonPhase:
  52. earth = get_skf_objects()["earth"]
  53. moon = get_skf_objects()["moon"]
  54. sun = get_skf_objects()["sun"]
  55. def moon_phase_at(time: Time):
  56. time._nutation_angles = get_iau2000b(time)
  57. current_earth = earth.at(time)
  58. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon("date")
  59. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon("date")
  60. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  61. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  62. today = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day)
  63. time1 = get_timescale().utc(
  64. compute_date.year, compute_date.month, compute_date.day - 10
  65. )
  66. time2 = get_timescale().utc(
  67. compute_date.year, compute_date.month, compute_date.day + 10
  68. )
  69. try:
  70. times, phase = find_discrete(time1, time2, moon_phase_at)
  71. except EphemerisRangeError as error:
  72. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  73. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  74. start = datetime.date(start.year, start.month, start.day) + datetime.timedelta(
  75. days=12
  76. )
  77. end = datetime.date(end.year, end.month, end.day) - datetime.timedelta(days=12)
  78. raise OutOfRangeDateError(start, end) from error
  79. return _get_skyfield_to_moon_phase(times, phase, today)
  80. def get_ephemerides(
  81. position: Position, date: datetime.date = datetime.date.today(), timezone: int = 0
  82. ) -> [AsterEphemerides]:
  83. ephemerides = []
  84. def get_angle(for_aster: Object):
  85. def fun(time: Time) -> float:
  86. return (
  87. position.get_planet_topos()
  88. .at(time)
  89. .observe(for_aster.get_skyfield_object())
  90. .apparent()
  91. .altaz()[0]
  92. .degrees
  93. )
  94. fun.rough_period = 1.0
  95. return fun
  96. def is_risen(for_aster: Object):
  97. def fun(time: Time) -> bool:
  98. return get_angle(for_aster)(time) > RISEN_ANGLE
  99. fun.rough_period = 0.5
  100. return fun
  101. start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)
  102. end_time = get_timescale().utc(
  103. date.year, date.month, date.day, 23 - timezone, 59, 59
  104. )
  105. try:
  106. for aster in ASTERS:
  107. rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))
  108. try:
  109. culmination_time, _ = find_maxima(
  110. start_time,
  111. end_time,
  112. f=get_angle(aster),
  113. epsilon=1.0 / 3600 / 24,
  114. num=12,
  115. )
  116. culmination_time = (
  117. culmination_time[0] if len(culmination_time) > 0 else None
  118. )
  119. except ValueError:
  120. culmination_time = None
  121. if len(rise_times) == 2:
  122. rise_time = rise_times[0 if arr[0] else 1]
  123. set_time = rise_times[1 if not arr[1] else 0]
  124. else:
  125. rise_time = rise_times[0] if arr[0] else None
  126. set_time = rise_times[0] if not arr[0] else None
  127. # Convert the Time instances to Python datetime objects
  128. if rise_time is not None:
  129. rise_time = translate_to_timezone(
  130. rise_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  131. )
  132. if culmination_time is not None:
  133. culmination_time = translate_to_timezone(
  134. culmination_time.utc_datetime().replace(microsecond=0),
  135. to_tz=timezone,
  136. )
  137. if set_time is not None:
  138. set_time = translate_to_timezone(
  139. set_time.utc_datetime().replace(microsecond=0), to_tz=timezone
  140. )
  141. ephemerides.append(
  142. AsterEphemerides(rise_time, culmination_time, set_time, aster=aster)
  143. )
  144. except EphemerisRangeError as error:
  145. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  146. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  147. start = datetime.date(start.year, start.month, start.day + 1)
  148. end = datetime.date(end.year, end.month, end.day - 1)
  149. raise OutOfRangeDateError(start, end) from error
  150. return ephemerides