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.
 
 
 
 

146 regels
5.9 KiB

  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 import almanac
  20. from skyfield.searchlib import find_discrete, find_maxima
  21. from skyfield.timelib import Time
  22. from skyfield.constants import tau
  23. from .data import Object, Position, AsterEphemerides, MoonPhase, ASTERS, skyfield_to_moon_phase
  24. from .core import get_skf_objects, get_timescale, get_iau2000b
  25. RISEN_ANGLE = -0.8333
  26. class EphemeridesComputer:
  27. def __init__(self, position: Union[Position, None]):
  28. if position is not None:
  29. position.observation_planet = get_skf_objects()['earth']
  30. self.position = position
  31. def get_sun(self, start_time, end_time) -> dict:
  32. times, is_risen = find_discrete(start_time,
  33. end_time,
  34. almanac.sunrise_sunset(get_skf_objects(), self.position))
  35. sunrise = times[0] if is_risen[0] else times[1]
  36. sunset = times[1] if not is_risen[1] else times[0]
  37. return {'rise': sunrise, 'set': sunset}
  38. @staticmethod
  39. def get_moon_phase(compute_date: datetime.date) -> MoonPhase:
  40. earth = get_skf_objects()['earth']
  41. moon = get_skf_objects()['moon']
  42. sun = get_skf_objects()['sun']
  43. def moon_phase_at(time: Time):
  44. time._nutation_angles = get_iau2000b(time)
  45. current_earth = earth.at(time)
  46. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon('date')
  47. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon('date')
  48. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  49. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  50. today = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day)
  51. time1 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day - 10)
  52. time2 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day + 10)
  53. times, phase = find_discrete(time1, time2, moon_phase_at)
  54. return skyfield_to_moon_phase(times, phase, today)
  55. @staticmethod
  56. def get_asters_ephemerides_for_aster(aster, date: datetime.date, position: Position) -> Object:
  57. skyfield_aster = get_skf_objects()[aster.skyfield_name]
  58. def get_angle(time: Time) -> float:
  59. return position.get_planet_topos().at(time).observe(skyfield_aster).apparent().altaz()[0].degrees
  60. def is_risen(time: Time) -> bool:
  61. return get_angle(time) > RISEN_ANGLE
  62. get_angle.rough_period = 1.0
  63. is_risen.rough_period = 0.5
  64. start_time = get_timescale().utc(date.year, date.month, date.day)
  65. end_time = get_timescale().utc(date.year, date.month, date.day, 23, 59, 59)
  66. rise_times, arr = find_discrete(start_time, end_time, is_risen)
  67. try:
  68. culmination_time, _ = find_maxima(start_time, end_time, f=get_angle, epsilon=1./3600/24, num=12)
  69. culmination_time = culmination_time[0] if len(culmination_time) > 0 else None
  70. except ValueError:
  71. culmination_time = None
  72. if len(rise_times) == 2:
  73. rise_time = rise_times[0 if arr[0] else 1]
  74. set_time = rise_times[1 if not arr[1] else 0]
  75. else:
  76. rise_time = rise_times[0] if arr[0] else None
  77. set_time = rise_times[0] if not arr[0] else None
  78. # Convert the Time instances to Python datetime objects
  79. if rise_time is not None:
  80. rise_time = rise_time.utc_datetime().replace(microsecond=0)
  81. if culmination_time is not None:
  82. culmination_time = culmination_time.utc_datetime().replace(microsecond=0)
  83. if set_time is not None:
  84. set_time = set_time.utc_datetime().replace(microsecond=0) if set_time is not None else None
  85. aster.ephemerides = AsterEphemerides(rise_time, culmination_time, set_time)
  86. return aster
  87. @staticmethod
  88. def is_leap_year(year: int) -> bool:
  89. return (year % 4 == 0 and year % 100 > 0) or (year % 400 == 0)
  90. def compute_ephemerides(self, compute_date: datetime.date) -> dict:
  91. return {'moon_phase': self.get_moon_phase(compute_date),
  92. 'details': [self.get_asters_ephemerides_for_aster(aster, compute_date, self.position)
  93. for aster in ASTERS] if self.position is not None else []}
  94. @staticmethod
  95. def get_seasons(year: int) -> dict:
  96. start_time = get_timescale().utc(year, 1, 1)
  97. end_time = get_timescale().utc(year, 12, 31)
  98. times, almanac_seasons = find_discrete(start_time, end_time, almanac.seasons(get_skf_objects()))
  99. seasons = {}
  100. for time, almanac_season in zip(times, almanac_seasons):
  101. if almanac_season == 0:
  102. season = 'MARCH'
  103. elif almanac_season == 1:
  104. season = 'JUNE'
  105. elif almanac_season == 2:
  106. season = 'SEPTEMBER'
  107. elif almanac_season == 3:
  108. season = 'DECEMBER'
  109. else:
  110. raise AssertionError
  111. seasons[season] = time.utc_iso()
  112. return seasons