您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符

ephemerides.py 4.0 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  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 skyfield.searchlib import find_discrete, find_maxima
  19. from skyfield.timelib import Time
  20. from skyfield.constants import tau
  21. from .data import Position, AsterEphemerides, MoonPhase, Object, ASTERS, skyfield_to_moon_phase
  22. from .core import get_skf_objects, get_timescale, get_iau2000b
  23. RISEN_ANGLE = -0.8333
  24. def get_moon_phase(compute_date: datetime.date) -> MoonPhase:
  25. earth = get_skf_objects()['earth']
  26. moon = get_skf_objects()['moon']
  27. sun = get_skf_objects()['sun']
  28. def moon_phase_at(time: Time):
  29. time._nutation_angles = get_iau2000b(time)
  30. current_earth = earth.at(time)
  31. _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon('date')
  32. _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon('date')
  33. return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)
  34. moon_phase_at.rough_period = 7.0 # one lunar phase per week
  35. today = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day)
  36. time1 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day - 10)
  37. time2 = get_timescale().utc(compute_date.year, compute_date.month, compute_date.day + 10)
  38. times, phase = find_discrete(time1, time2, moon_phase_at)
  39. return skyfield_to_moon_phase(times, phase, today)
  40. def get_ephemerides(date: datetime.date, position: Position) -> [AsterEphemerides]:
  41. ephemerides = []
  42. def get_angle(for_aster: Object):
  43. def fun(time: Time) -> float:
  44. return position.get_planet_topos().at(time).observe(for_aster.get_skyfield_object()).apparent().altaz()[0]\
  45. .degrees
  46. fun.rough_period = 1.0
  47. return fun
  48. def is_risen(for_aster: Object):
  49. def fun(time: Time) -> bool:
  50. return get_angle(for_aster)(time) > RISEN_ANGLE
  51. fun.rough_period = 0.5
  52. return fun
  53. start_time = get_timescale().utc(date.year, date.month, date.day)
  54. end_time = get_timescale().utc(date.year, date.month, date.day, 23, 59, 59)
  55. for aster in ASTERS:
  56. rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))
  57. try:
  58. culmination_time, _ = find_maxima(start_time, end_time, f=get_angle(aster), epsilon=1./3600/24, num=12)
  59. culmination_time = culmination_time[0] if len(culmination_time) > 0 else None
  60. except ValueError:
  61. culmination_time = None
  62. if len(rise_times) == 2:
  63. rise_time = rise_times[0 if arr[0] else 1]
  64. set_time = rise_times[1 if not arr[1] else 0]
  65. else:
  66. rise_time = rise_times[0] if arr[0] else None
  67. set_time = rise_times[0] if not arr[0] else None
  68. # Convert the Time instances to Python datetime objects
  69. if rise_time is not None:
  70. rise_time = rise_time.utc_datetime().replace(microsecond=0)
  71. if culmination_time is not None:
  72. culmination_time = culmination_time.utc_datetime().replace(microsecond=0)
  73. if set_time is not None:
  74. set_time = set_time.utc_datetime().replace(microsecond=0) if set_time is not None else None
  75. ephemerides.append(AsterEphemerides(rise_time, culmination_time, set_time, aster=aster))
  76. return ephemerides