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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141
  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. from datetime import date as date_type
  18. from skyfield.timelib import Time
  19. from skyfield.searchlib import find_discrete, find_maxima
  20. from numpy import pi
  21. from .data import Event, Star, Planet, ASTERS
  22. from .core import get_timescale, get_skf_objects, flatten_list
  23. def _search_conjunction(start_time: Time, end_time: Time) -> [Event]:
  24. earth = get_skf_objects()['earth']
  25. aster1 = None
  26. aster2 = None
  27. def is_in_conjunction(time: Time):
  28. earth_pos = earth.at(time)
  29. _, aster1_lon, _ = earth_pos.observe(aster1.get_skyfield_object()).apparent().ecliptic_latlon()
  30. _, aster2_lon, _ = earth_pos.observe(aster2.get_skyfield_object()).apparent().ecliptic_latlon()
  31. return ((aster1_lon.radians - aster2_lon.radians) / pi % 2.0).astype('int8') == 0
  32. is_in_conjunction.rough_period = 60.0
  33. computed = []
  34. conjunctions = []
  35. for aster1 in ASTERS:
  36. # Ignore the Sun
  37. if isinstance(aster1, Star):
  38. continue
  39. for aster2 in ASTERS:
  40. if isinstance(aster2, Star) or aster2 == aster1 or aster2 in computed:
  41. continue
  42. times, is_conjs = find_discrete(start_time, end_time, is_in_conjunction)
  43. for i, time in enumerate(times):
  44. if is_conjs[i]:
  45. aster1_pos = (aster1.get_skyfield_object() - earth).at(time)
  46. aster2_pos = (aster2.get_skyfield_object() - earth).at(time)
  47. distance = aster1_pos.separation_from(aster2_pos).degrees
  48. if distance - aster2.get_apparent_radius(time, earth) < aster1.get_apparent_radius(time, earth):
  49. occulting_aster = [aster1,
  50. aster2] if aster1_pos.distance().km < aster2_pos.distance().km else [aster2,
  51. aster1]
  52. conjunctions.append(Event('OCCULTATION', occulting_aster, time.utc_datetime()))
  53. else:
  54. conjunctions.append(Event('CONJUNCTION', [aster1, aster2], time.utc_datetime()))
  55. computed.append(aster1)
  56. return conjunctions
  57. def _search_oppositions(start_time: Time, end_time: Time) -> [Event]:
  58. earth = get_skf_objects()['earth']
  59. sun = get_skf_objects()['sun']
  60. aster = None
  61. def is_oppositing(time: Time) -> [bool]:
  62. earth_pos = earth.at(time)
  63. sun_pos = earth_pos.observe(sun).apparent() # Never do this without eyes protection!
  64. aster_pos = earth_pos.observe(get_skf_objects()[aster.skyfield_name]).apparent()
  65. _, lon1, _ = sun_pos.ecliptic_latlon()
  66. _, lon2, _ = aster_pos.ecliptic_latlon()
  67. return (lon1.degrees - lon2.degrees) > 180
  68. is_oppositing.rough_period = 1.0
  69. events = []
  70. for aster in ASTERS:
  71. if not isinstance(aster, Planet) or aster.skyfield_name in ['MERCURY', 'VENUS']:
  72. continue
  73. times, _ = find_discrete(start_time, end_time, is_oppositing)
  74. for time in times:
  75. events.append(Event('OPPOSITION', [aster], time.utc_datetime()))
  76. return events
  77. def _search_maximal_elongations(start_time: Time, end_time: Time) -> [Event]:
  78. earth = get_skf_objects()['earth']
  79. sun = get_skf_objects()['sun']
  80. aster = None
  81. def get_elongation(time: Time):
  82. sun_pos = (sun - earth).at(time)
  83. aster_pos = (aster.get_skyfield_object() - earth).at(time)
  84. separation = sun_pos.separation_from(aster_pos)
  85. return separation.degrees
  86. get_elongation.rough_period = 1.0
  87. events = []
  88. for aster in ASTERS:
  89. if aster.skyfield_name not in ['MERCURY', 'VENUS']:
  90. continue
  91. times, elongations = find_maxima(start_time, end_time, f=get_elongation, epsilon=1./24/3600, num=12)
  92. for i, time in enumerate(times):
  93. elongation = elongations[i]
  94. events.append(Event('MAXIMAL_ELONGATION', [aster], time.utc_datetime(),
  95. details='{:.3n}°'.format(elongation)))
  96. return events
  97. def search_events(date: date_type) -> [Event]:
  98. start_time = get_timescale().utc(date.year, date.month, date.day)
  99. end_time = get_timescale().utc(date.year, date.month, date.day + 1)
  100. return sorted(flatten_list([
  101. _search_oppositions(start_time, end_time),
  102. _search_conjunction(start_time, end_time),
  103. _search_maximal_elongations(start_time, end_time)
  104. ]), key=lambda event: event.start_time)