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 6.6 KiB

3 years ago
3 years ago
3 years ago
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  1. #!/usr/bin/env python3
  2. # Kosmorro - Compute The Next Ephemerides
  3. # Copyright (C) 2019 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. import datetime
  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 .data import Position, AsterEphemerides, MoonPhase, Object, ASTERS
  24. from .dateutil import translate_to_timezone
  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(times: [Time], vals: [int], now: Time) -> Union[MoonPhase, None]:
  30. tomorrow = get_timescale().utc(now.utc_datetime().year, now.utc_datetime().month, now.utc_datetime().day + 1)
  31. phases = list(MoonPhaseType)
  32. current_phase = None
  33. current_phase_time = None
  34. next_phase_time = None
  35. i = 0
  36. if len(times) == 0:
  37. return None
  38. for i, time in enumerate(times):
  39. if now.utc_iso() <= time.utc_iso():
  40. if vals[i] in [0, 2, 4, 6]:
  41. if time.utc_datetime() < tomorrow.utc_datetime():
  42. current_phase_time = time
  43. current_phase = phases[vals[i]]
  44. else:
  45. i -= 1
  46. current_phase_time = None
  47. current_phase = phases[vals[i]]
  48. else:
  49. current_phase = phases[vals[i]]
  50. break
  51. for j in range(i + 1, len(times)):
  52. if vals[j] in [0, 2, 4, 6]:
  53. next_phase_time = times[j]
  54. break
  55. return MoonPhase(current_phase,
  56. current_phase_time.utc_datetime() if current_phase_time is not None else None,
  57. next_phase_time.utc_datetime() if next_phase_time is not None else None)
  58. def get_moon_phase(compute_date: datetime.date, timezone: int = 0) -> MoonPhase:
  59. earth = get_skf_objects()['earth']
  60. moon = get_skf_objects()['moon']
  61. sun = get_skf_objects()['sun']
  62. def moon_phase_at(time: Time):
  63. time._nutation_angles = get_iau2000b(time)
  64. current_earth = earth.at(time)
  65. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon('date')
  66. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon('date')
  67. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  68. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  69. today = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day)
  70. time1 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day - 10)
  71. time2 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day + 10)
  72. try:
  73. times, phase = find_discrete(time1, time2, moon_phase_at)
  74. except EphemerisRangeError as error:
  75. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  76. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  77. start = datetime.date(start.year, start.month, start.day) + datetime.timedelta(days=12)
  78. end = datetime.date(end.year, end.month, end.day) - datetime.timedelta(days=12)
  79. raise OutOfRangeDateError(start, end) from error
  80. return _get_skyfield_to_moon_phase(times, phase, today)
  81. def get_ephemerides(date: datetime.date, position: Position, timezone: int = 0) -> [AsterEphemerides]:
  82. ephemerides = []
  83. def get_angle(for_aster: Object):
  84. def fun(time: Time) -> float:
  85. return position.get_planet_topos().at(time).observe(for_aster.get_skyfield_object()).apparent().altaz()[0]\
  86. .degrees
  87. fun.rough_period = 1.0
  88. return fun
  89. def is_risen(for_aster: Object):
  90. def fun(time: Time) -> bool:
  91. return get_angle(for_aster)(time) > RISEN_ANGLE
  92. fun.rough_period = 0.5
  93. return fun
  94. start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)
  95. end_time = get_timescale().utc(date.year, date.month, date.day, 23 - timezone, 59, 59)
  96. try:
  97. for aster in ASTERS:
  98. rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))
  99. try:
  100. culmination_time, _ = find_maxima(start_time, end_time, f=get_angle(aster), epsilon=1./3600/24, num=12)
  101. culmination_time = culmination_time[0] if len(culmination_time) > 0 else None
  102. except ValueError:
  103. culmination_time = None
  104. if len(rise_times) == 2:
  105. rise_time = rise_times[0 if arr[0] else 1]
  106. set_time = rise_times[1 if not arr[1] else 0]
  107. else:
  108. rise_time = rise_times[0] if arr[0] else None
  109. set_time = rise_times[0] if not arr[0] else None
  110. # Convert the Time instances to Python datetime objects
  111. if rise_time is not None:
  112. rise_time = translate_to_timezone(rise_time.utc_datetime().replace(microsecond=0),
  113. to_tz=timezone)
  114. if culmination_time is not None:
  115. culmination_time = translate_to_timezone(culmination_time.utc_datetime().replace(microsecond=0),
  116. to_tz=timezone)
  117. if set_time is not None:
  118. set_time = translate_to_timezone(set_time.utc_datetime().replace(microsecond=0),
  119. to_tz=timezone)
  120. ephemerides.append(AsterEphemerides(rise_time, culmination_time, set_time, aster=aster))
  121. except EphemerisRangeError as error:
  122. start = translate_to_timezone(error.start_time.utc_datetime(), timezone)
  123. end = translate_to_timezone(error.end_time.utc_datetime(), timezone)
  124. start = datetime.date(start.year, start.month, start.day + 1)
  125. end = datetime.date(end.year, end.month, end.day - 1)
  126. raise OutOfRangeDateError(start, end) from error
  127. return ephemerides