1#!/usr/bin/env python3 2 3# Kosmorrolib - The Library To Compute Your Ephemerides 4# Copyright (C) 2021 Jérôme Deuchnord <jerome@deuchnord.fr> 5# 6# This program is free software: you can redistribute it and/or modify 7# it under the terms of the GNU Affero General Public License as 8# published by the Free Software Foundation, either version 3 of the 9# License, or (at your option) any later version. 10# 11# This program is distributed in the hope that it will be useful, 12# but WITHOUT ANY WARRANTY; without even the implied warranty of 13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14# GNU Affero General Public License for more details. 15# 16# You should have received a copy of the GNU Affero General Public License 17# along with this program. If not, see <https://www.gnu.org/licenses/>. 18 19from datetime import date, timedelta 20 21from skyfield.errors import EphemerisRangeError 22from skyfield.timelib import Time 23from skyfield.searchlib import find_discrete, find_maxima, find_minima 24from skyfield.units import Angle 25from skyfield import almanac, eclipselib 26from numpy import pi 27 28from kosmorrolib.model import ( 29 Object, 30 Event, 31 Object, 32 Star, 33 Planet, 34 get_aster, 35 ASTERS, 36 EARTH, 37) 38from kosmorrolib.dateutil import translate_to_timezone 39from kosmorrolib.enum import EventType, ObjectIdentifier, SeasonType, LunarEclipseType 40from kosmorrolib.exceptions import OutOfRangeDateError 41from kosmorrolib.core import get_timescale, get_skf_objects, flatten_list 42 43 44def _search_conjunctions_occultations( 45 start_time: Time, end_time: Time, timezone: int 46) -> [Event]: 47 """Function to search conjunction. 48 49 **Warning:** this is an internal function, not intended for use by end-developers. 50 51 Will return MOON and VENUS opposition on 2021-06-12: 52 53 >>> conjunction = _search_conjunctions_occultations(get_timescale().utc(2021, 6, 12), get_timescale().utc(2021, 6, 13), 0) 54 >>> len(conjunction) 55 1 56 >>> conjunction[0].objects 57 [<Object type=SATELLITE name=MOON />, <Object type=PLANET name=VENUS />] 58 59 Will return nothing if no conjunction happens: 60 61 >>> _search_conjunctions_occultations(get_timescale().utc(2021, 6, 17),get_timescale().utc(2021, 6, 18), 0) 62 [] 63 64 This function detects occultations too: 65 66 >>> _search_conjunctions_occultations(get_timescale().utc(2021, 4, 17),get_timescale().utc(2021, 4, 18), 0) 67 [<Event type=OCCULTATION objects=[<Object type=SATELLITE name=MOON />, <Object type=PLANET name=MARS />] start=2021-04-17 12:08:16.115650+00:00 end=None details=None />] 68 """ 69 earth = get_skf_objects()["earth"] 70 aster1 = None 71 aster2 = None 72 73 def is_in_conjunction(time: Time): 74 earth_pos = earth.at(time) 75 _, aster1_lon, _ = ( 76 earth_pos.observe(aster1.skyfield_object).apparent().ecliptic_latlon() 77 ) 78 _, aster2_lon, _ = ( 79 earth_pos.observe(aster2.skyfield_object).apparent().ecliptic_latlon() 80 ) 81 82 return ((aster1_lon.radians - aster2_lon.radians) / pi % 2.0).astype( 83 "int8" 84 ) == 0 85 86 is_in_conjunction.rough_period = 60.0 87 88 computed = [] 89 events = [] 90 91 for aster1 in ASTERS: 92 # Ignore the Sun 93 if isinstance(aster1, Star): 94 continue 95 96 for aster2 in ASTERS: 97 if isinstance(aster2, Star) or aster2 == aster1 or aster2 in computed: 98 continue 99 100 times, is_conjs = find_discrete(start_time, end_time, is_in_conjunction) 101 102 for i, time in enumerate(times): 103 if is_conjs[i]: 104 aster1_pos = (aster1.skyfield_object - earth).at(time) 105 aster2_pos = (aster2.skyfield_object - earth).at(time) 106 distance = aster1_pos.separation_from(aster2_pos).degrees 107 108 if ( 109 distance - aster2.get_apparent_radius(time).degrees 110 < aster1.get_apparent_radius(time).degrees 111 ): 112 occulting_aster = ( 113 [aster1, aster2] 114 if aster1_pos.distance().km < aster2_pos.distance().km 115 else [aster2, aster1] 116 ) 117 118 events.append( 119 Event( 120 EventType.OCCULTATION, 121 occulting_aster, 122 translate_to_timezone(time.utc_datetime(), timezone), 123 ) 124 ) 125 else: 126 events.append( 127 Event( 128 EventType.CONJUNCTION, 129 [aster1, aster2], 130 translate_to_timezone(time.utc_datetime(), timezone), 131 ) 132 ) 133 134 computed.append(aster1) 135 136 return events 137 138 139def _search_oppositions(start_time: Time, end_time: Time, timezone: int) -> [Event]: 140 """Function to search oppositions. 141 142 **Warning:** this is an internal function, not intended for use by end-developers. 143 144 Will return Mars opposition on 2020-10-13: 145 146 >>> oppositions = _search_oppositions(get_timescale().utc(2020, 10, 13), get_timescale().utc(2020, 10, 14), 0) 147 >>> len(oppositions) 148 1 149 >>> oppositions[0].objects[0] 150 <Object type=PLANET name=MARS /> 151 152 Will return nothing if no opposition happens: 153 154 >>> _search_oppositions(get_timescale().utc(2021, 3, 20), get_timescale().utc(2021, 3, 21), 0) 155 [] 156 """ 157 earth = get_skf_objects()["earth"] 158 sun = get_skf_objects()["sun"] 159 aster = None 160 161 def is_oppositing(time: Time) -> [bool]: 162 diff = get_angle(time) 163 return diff > 180 164 165 def get_angle(time: Time): 166 earth_pos = earth.at(time) 167 168 sun_pos = earth_pos.observe( 169 sun 170 ).apparent() # Never do this without eyes protection! 171 aster_pos = earth_pos.observe(aster.skyfield_object).apparent() 172 173 _, lon1, _ = sun_pos.ecliptic_latlon() 174 _, lon2, _ = aster_pos.ecliptic_latlon() 175 176 return lon1.degrees - lon2.degrees 177 178 is_oppositing.rough_period = 1.0 179 events = [] 180 181 for aster in ASTERS: 182 if not isinstance(aster, Planet) or aster.identifier in [ 183 ObjectIdentifier.MERCURY, 184 ObjectIdentifier.VENUS, 185 ]: 186 continue 187 188 times, _ = find_discrete(start_time, end_time, is_oppositing) 189 for time in times: 190 if get_angle(time) < 0: 191 # If the angle is negative, then it is actually a false positive. 192 # Just ignoring it. 193 continue 194 195 events.append( 196 Event( 197 EventType.OPPOSITION, 198 [aster], 199 translate_to_timezone(time.utc_datetime(), timezone), 200 ) 201 ) 202 203 return events 204 205 206def _search_maximal_elongations( 207 start_time: Time, end_time: Time, timezone: int 208) -> [Event]: 209 """Function to search oppositions. 210 211 **Warning:** this is an internal function, not intended for use by end-developers. 212 213 Will return Mercury maimum elogation for September 14, 2021: 214 215 >>> get_events(date(2021, 9, 14)) 216 [<Event type=MAXIMAL_ELONGATION objects=[<Object type=PLANET name=MERCURY />] start=2021-09-14 04:13:46.664879+00:00 end=None details={'deg': 26.8} />] 217 """ 218 earth = get_skf_objects()["earth"] 219 sun = get_skf_objects()["sun"] 220 221 def get_elongation(planet: Object): 222 def f(time: Time): 223 sun_pos = (sun - earth).at(time) 224 aster_pos = (planet.skyfield_object - earth).at(time) 225 separation = sun_pos.separation_from(aster_pos) 226 return separation.degrees 227 228 f.rough_period = 1.0 229 230 return f 231 232 events = [] 233 234 for identifier in [ObjectIdentifier.MERCURY, ObjectIdentifier.VENUS]: 235 planet = get_aster(identifier) 236 times, elongations = find_maxima( 237 start_time, 238 end_time, 239 f=get_elongation(planet), 240 epsilon=1.0 / 24 / 3600, 241 num=12, 242 ) 243 244 for i, time in enumerate(times): 245 elongation = round(elongations[i], 1) 246 events.append( 247 Event( 248 EventType.MAXIMAL_ELONGATION, 249 [planet], 250 translate_to_timezone(time.utc_datetime(), timezone), 251 details={"deg": elongation}, 252 ) 253 ) 254 255 return events 256 257 258def _get_distance(to_aster: Object, from_aster: Object): 259 def get_distance(time: Time): 260 from_pos = from_aster.skyfield_object.at(time) 261 to_pos = from_pos.observe(to_aster.skyfield_object) 262 263 return to_pos.distance().km 264 265 get_distance.rough_period = 1.0 266 267 return get_distance 268 269 270def _search_apogee(to_aster: Object, from_aster: Object = EARTH) -> callable: 271 """Search for moon apogee 272 273 **Warning:** this is an internal function, not intended for use by end-developers. 274 275 Get the moon apogee: 276 >>> _search_apogee(ASTERS[1])(get_timescale().utc(2021, 6, 8), get_timescale().utc(2021, 6, 9), 0) 277 [<Event type=APOGEE objects=[<Object type=SATELLITE name=MOON />] start=2021-06-08 02:39:40.165271+00:00 end=None details={'distance_km': 406211.04850197025} />] 278 279 Get the Earth's apogee: 280 >>> _search_apogee(EARTH, from_aster=ASTERS[0])(get_timescale().utc(2021, 7, 5), get_timescale().utc(2021, 7, 6), 0) 281 [<Event type=APOGEE objects=[<Object type=PLANET name=EARTH />] start=2021-07-05 22:35:42.148792+00:00 end=None details={'distance_km': 152100521.91712126} />] 282 """ 283 284 def f(start_time: Time, end_time: Time, timezone: int) -> [Event]: 285 events = [] 286 287 times, distances = find_maxima( 288 start_time, 289 end_time, 290 f=_get_distance(to_aster, from_aster), 291 epsilon=1.0 / 24 / 60, 292 ) 293 294 for i, time in enumerate(times): 295 events.append( 296 Event( 297 EventType.APOGEE, 298 [to_aster], 299 translate_to_timezone(time.utc_datetime(), timezone), 300 details={"distance_km": distances[i]}, 301 ) 302 ) 303 304 return events 305 306 return f 307 308 309def _search_perigee(aster: Object, from_aster: Object = EARTH) -> callable: 310 """Search for moon perigee 311 312 **Warning:** this is an internal function, not intended for use by end-developers. 313 314 Get the moon perigee: 315 >>> _search_perigee(ASTERS[1])(get_timescale().utc(2021, 5, 26), get_timescale().utc(2021, 5, 27), 0) 316 [<Event type=PERIGEE objects=[<Object type=SATELLITE name=MOON />] start=2021-05-26 01:56:01.983455+00:00 end=None details={'distance_km': 357313.9680798693} />] 317 318 Get the Earth's perigee: 319 >>> _search_perigee(EARTH, from_aster=ASTERS[0])(get_timescale().utc(2021, 1, 2), get_timescale().utc(2021, 1, 3), 0) 320 [<Event type=PERIGEE objects=[<Object type=PLANET name=EARTH />] start=2021-01-02 13:59:00.495905+00:00 end=None details={'distance_km': 147093166.1686309} />] 321 """ 322 323 def f(start_time: Time, end_time: Time, timezone: int) -> [Event]: 324 events = [] 325 326 times, distances = find_minima( 327 start_time, 328 end_time, 329 f=_get_distance(aster, from_aster), 330 epsilon=1.0 / 24 / 60, 331 ) 332 333 for i, time in enumerate(times): 334 events.append( 335 Event( 336 EventType.PERIGEE, 337 [aster], 338 translate_to_timezone(time.utc_datetime(), timezone), 339 details={"distance_km": distances[i]}, 340 ) 341 ) 342 343 return events 344 345 return f 346 347 348def _search_earth_season_change( 349 start_time: Time, end_time: Time, timezone: int 350) -> [Event]: 351 """Function to find earth season change event. 352 353 **Warning:** this is an internal function, not intended for use by end-developers. 354 355 Will return JUNE SOLSTICE on 2020/06/20: 356 357 >>> season_change = _search_earth_season_change(get_timescale().utc(2020, 6, 20), get_timescale().utc(2020, 6, 21), 0) 358 >>> len(season_change) 359 1 360 >>> season_change[0].event_type 361 <EventType.SEASON_CHANGE: 7> 362 >>> season_change[0].details 363 {'season': <SeasonType.JUNE_SOLSTICE: 1>} 364 365 Will return nothing if there is no season change event in the period of time being calculated: 366 367 >>> _search_earth_season_change(get_timescale().utc(2021, 6, 17), get_timescale().utc(2021, 6, 18), 0) 368 [] 369 """ 370 events = [] 371 event_time, event_id = almanac.find_discrete( 372 start_time, end_time, almanac.seasons(get_skf_objects()) 373 ) 374 if len(event_time) == 0: 375 return [] 376 events.append( 377 Event( 378 EventType.SEASON_CHANGE, 379 [], 380 translate_to_timezone(event_time.utc_datetime()[0], timezone), 381 details={"season": SeasonType(event_id[0])}, 382 ) 383 ) 384 return events 385 386 387def _search_lunar_eclipse(start_time: Time, end_time: Time, timezone: int) -> [Event]: 388 """Function to detect lunar eclipses. 389 390 **Warning:** this is an internal function, not intended for use by end-developers. 391 392 Will return a total lunar eclipse for 2021-05-26: 393 394 >>> _search_lunar_eclipse(get_timescale().utc(2021, 5, 26), get_timescale().utc(2021, 5, 27), 0) 395 [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2021-05-26 08:47:54.795821+00:00 end=2021-05-26 13:49:34.353411+00:00 details={'type': <LunarEclipseType.TOTAL: 2>, 'maximum': datetime.datetime(2021, 5, 26, 11, 18, 42, 328842, tzinfo=datetime.timezone.utc)} />] 396 397 >>> _search_lunar_eclipse(get_timescale().utc(2019, 7, 16), get_timescale().utc(2019, 7, 17), 0) 398 [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2019-07-16 18:39:53.391337+00:00 end=2019-07-17 00:21:51.378940+00:00 details={'type': <LunarEclipseType.PARTIAL: 1>, 'maximum': datetime.datetime(2019, 7, 16, 21, 30, 44, 170096, tzinfo=datetime.timezone.utc)} />] 399 400 >>> _search_lunar_eclipse(get_timescale().utc(2017, 2, 11), get_timescale().utc(2017, 2, 12), 0) 401 [<Event type=LUNAR_ECLIPSE objects=[<Object type=SATELLITE name=MOON />] start=2017-02-10 22:02:59.016572+00:00 end=2017-02-11 03:25:07.627886+00:00 details={'type': <LunarEclipseType.PENUMBRAL: 0>, 'maximum': datetime.datetime(2017, 2, 11, 0, 43, 51, 793786, tzinfo=datetime.timezone.utc)} />] 402 """ 403 moon = get_aster(ObjectIdentifier.MOON) 404 events = [] 405 t, y, details = eclipselib.lunar_eclipses(start_time, end_time, get_skf_objects()) 406 407 for ti, yi in zip(t, y): 408 penumbra_radius = Angle(radians=details["penumbra_radius_radians"][0]) 409 _, max_lon, _ = ( 410 EARTH.skyfield_object.at(ti) 411 .observe(moon.skyfield_object) 412 .apparent() 413 .ecliptic_latlon() 414 ) 415 416 def is_in_penumbra(time: Time): 417 _, lon, _ = ( 418 EARTH.skyfield_object.at(time) 419 .observe(moon.skyfield_object) 420 .apparent() 421 .ecliptic_latlon() 422 ) 423 424 moon_radius = details["moon_radius_radians"] 425 426 return ( 427 abs(max_lon.radians - lon.radians) 428 < penumbra_radius.radians + moon_radius 429 ) 430 431 is_in_penumbra.rough_period = 60.0 432 433 search_start_time = get_timescale().from_datetime( 434 start_time.utc_datetime() - timedelta(days=1) 435 ) 436 search_end_time = get_timescale().from_datetime( 437 end_time.utc_datetime() + timedelta(days=1) 438 ) 439 440 eclipse_start, _ = find_discrete(search_start_time, ti, is_in_penumbra) 441 eclipse_end, _ = find_discrete(ti, search_end_time, is_in_penumbra) 442 443 events.append( 444 Event( 445 EventType.LUNAR_ECLIPSE, 446 [moon], 447 start_time=translate_to_timezone( 448 eclipse_start[0].utc_datetime(), timezone 449 ), 450 end_time=translate_to_timezone(eclipse_end[0].utc_datetime(), timezone), 451 details={ 452 "type": LunarEclipseType(yi), 453 "maximum": translate_to_timezone(ti.utc_datetime(), timezone), 454 }, 455 ) 456 ) 457 458 return events 459 460 461def get_events(for_date: date = date.today(), timezone: int = 0) -> [Event]: 462 """Calculate and return a list of events for the given date, adjusted to the given timezone if any. 463 464 Find events that happen on April 4th, 2020 (show hours in UTC): 465 466 >>> get_events(date(2020, 4, 4)) 467 [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-04 01:14:39.063308+00:00 end=None details=None />] 468 469 Find events that happen on April 4th, 2020 (show timezones in UTC+2): 470 471 >>> get_events(date(2020, 4, 4), 2) 472 [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-04 03:14:39.063267+02:00 end=None details=None />] 473 474 Find events that happen on April 3rd, 2020 (show timezones in UTC-2): 475 476 >>> get_events(date(2020, 4, 3), -2) 477 [<Event type=CONJUNCTION objects=[<Object type=PLANET name=MERCURY />, <Object type=PLANET name=NEPTUNE />] start=2020-04-03 23:14:39.063388-02:00 end=None details=None />] 478 479 If there is no events for the given date, then an empty list is returned: 480 481 >>> get_events(date(2021, 4, 20)) 482 [] 483 484 Note that the events can only be found for a date range. 485 Asking for the events with an out of range date will result in an exception: 486 487 >>> get_events(date(1000, 1, 1)) 488 Traceback (most recent call last): 489 ... 490 kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-07-28 and 2053-10-08 491 492 :param for_date: the date for which the events must be calculated 493 :param timezone: the timezone to adapt the results to. If not given, defaults to 0. 494 :return: a list of events found for the given date. 495 """ 496 497 start_time = get_timescale().utc( 498 for_date.year, for_date.month, for_date.day, -timezone 499 ) 500 end_time = get_timescale().utc( 501 for_date.year, for_date.month, for_date.day + 1, -timezone 502 ) 503 504 try: 505 found_events = [] 506 507 for fun in [ 508 _search_oppositions, 509 _search_conjunctions_occultations, 510 _search_maximal_elongations, 511 _search_apogee(ASTERS[1]), 512 _search_perigee(ASTERS[1]), 513 _search_apogee(EARTH, from_aster=ASTERS[0]), 514 _search_perigee(EARTH, from_aster=ASTERS[0]), 515 _search_earth_season_change, 516 _search_lunar_eclipse, 517 ]: 518 found_events.append(fun(start_time, end_time, timezone)) 519 520 return sorted(flatten_list(found_events), key=lambda event: event.start_time) 521 except EphemerisRangeError as error: 522 start_date = translate_to_timezone(error.start_time.utc_datetime(), timezone) 523 end_date = translate_to_timezone(error.end_time.utc_datetime(), timezone) 524 525 start_date = date(start_date.year, start_date.month, start_date.day) 526 end_date = date(end_date.year, end_date.month, end_date.day) 527 528 raise OutOfRangeDateError(start_date, end_date) from error 529