#!/usr/bin/env python3 # Kosmorrolib - The Library To Compute Your Ephemerides # Copyright (C) 2021 Jérôme Deuchnord # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as # published by the Free Software Foundation, either version 3 of the # License, or (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . from datetime import date, datetime, timedelta from typing import Union from skyfield.searchlib import find_discrete, find_maxima from skyfield.timelib import Time from skyfield.constants import tau from skyfield.errors import EphemerisRangeError from .model import Position, AsterEphemerides, MoonPhase, Object, ASTERS from .dateutil import translate_to_timezone, normalize_datetime from .core import get_skf_objects, get_timescale, get_iau2000b from .enum import MoonPhaseType from .exceptions import OutOfRangeDateError RISEN_ANGLE = -0.8333 def _get_skyfield_to_moon_phase( times: [Time], vals: [int], now: Time, timezone: int ) -> Union[MoonPhase, None]: tomorrow = get_timescale().utc( now.utc_datetime().year, now.utc_datetime().month, now.utc_datetime().day + 1 ) next_phase_time = None i = 0 # Find the current moon phase: for i, time in enumerate(times): if now.utc_datetime() <= time.utc_datetime(): if time.utc_datetime() >= tomorrow.utc_datetime(): i -= 1 break current_phase = MoonPhaseType(vals[i]) if current_phase in [ MoonPhaseType.NEW_MOON, MoonPhaseType.FIRST_QUARTER, MoonPhaseType.FULL_MOON, MoonPhaseType.LAST_QUARTER, ]: current_phase_time = translate_to_timezone(times[i].utc_datetime(), timezone) else: current_phase_time = None # Find the next moon phase for j in range(i + 1, len(times)): if vals[j] in [0, 2, 4, 6]: next_phase_time = translate_to_timezone(times[j].utc_datetime(), timezone) break return MoonPhase(current_phase, current_phase_time, next_phase_time) def get_moon_phase(for_date: date = date.today(), timezone: int = 0) -> MoonPhase: """Calculate and return the moon phase for the given date, adjusted to the given timezone if any. Get the moon phase for the 27 March, 2021: >>> get_moon_phase(date(2021, 3, 27)) When the moon phase is a new moon, a first quarter, a full moon or a last quarter, you get the exact time of its happening too: >>> get_moon_phase(datetime(2021, 3, 28)) Get the moon phase for the 27 March, 2021, in the UTC+2 timezone: >>> get_moon_phase(date(2021, 3, 27), timezone=2) Note that the moon phase can only be computed for a date range. Asking for the moon phase with an out of range date will result in an exception: >>> get_moon_phase(date(1000, 1, 1)) Traceback (most recent call last): ... kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-08-09 and 2053-09-26 """ earth = get_skf_objects()["earth"] moon = get_skf_objects()["moon"] sun = get_skf_objects()["sun"] def moon_phase_at(time: Time): time._nutation_angles = get_iau2000b(time) current_earth = earth.at(time) _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon("date") _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon("date") return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int) moon_phase_at.rough_period = 7.0 # one lunar phase per week today = get_timescale().utc(for_date.year, for_date.month, for_date.day) start_time = get_timescale().utc(for_date.year, for_date.month, for_date.day - 10) end_time = get_timescale().utc(for_date.year, for_date.month, for_date.day + 10) try: times, phases = find_discrete(start_time, end_time, moon_phase_at) return _get_skyfield_to_moon_phase(times, phases, today, timezone) except EphemerisRangeError as error: start = translate_to_timezone(error.start_time.utc_datetime(), timezone) end = translate_to_timezone(error.end_time.utc_datetime(), timezone) start = date(start.year, start.month, start.day) + timedelta(days=12) end = date(end.year, end.month, end.day) - timedelta(days=12) raise OutOfRangeDateError(start, end) from error def get_ephemerides( position: Position, for_date: date = date.today(), timezone: int = 0 ) -> [AsterEphemerides]: """Compute and return the ephemerides for the given position and date, adjusted to the given timezone if any. Compute the ephemerides for July 7th, 2022: >>> get_ephemerides(Position(36.6794, 4.8555), date(2022, 7, 7)) [>, >, >, >, >, >, >, >, >, >] Timezone can be optionnaly set to adapt the hours to your location: >>> get_ephemerides(Position(36.6794, 4.8555), date(2022, 7, 7), timezone=2) [>, >, >, >, >, >, >, >, >, >] Objects may not rise or set on the given date (e.g. they rise the previous day or set the next day). In this case, you will get `None` values on the rise or set time. If an objet does not rise nor set due to your latitude, then both rise and set will be `None`: >>> north_pole = Position(70, 20) >>> south_pole = Position(-70, 20) >>> get_ephemerides(north_pole, date(2021, 6, 20)) [>, >, >, >, >, >, >, >, >, >] >>> get_ephemerides(north_pole, date(2021, 12, 21)) [>, >, >, >, >, >, >, >, >, >] >>> get_ephemerides(south_pole, date(2021, 6, 20)) [>, >, >, >, >, >, >, >, >, >] >>> get_ephemerides(south_pole, date(2021, 12, 22)) [>, >, >, >, >, >, >, >, >, >] Please note: - The ephemerides can only be computed for a date range. Asking for the ephemerides with an out of range date will result in an exception: >>> get_ephemerides(Position(50.5824, 3.0624), date(1000, 1, 1)) Traceback (most recent call last): ... kosmorrolib.exceptions.OutOfRangeDateError: The date must be between 1899-07-29 and 2053-10-07 - The date given in parameter is considered as being given from the point of view of the given timezone. Using a timezone that does not correspond to the place's actual one can impact the returned times. """ def get_angle(for_aster: Object): def fun(time: Time) -> float: return ( position.get_planet_topos() .at(time) .observe(for_aster.skyfield_object) .apparent() .altaz()[0] .degrees ) fun.rough_period = 1.0 return fun def is_risen(for_aster: Object): def fun(time: Time) -> bool: return get_angle(for_aster)(time) > RISEN_ANGLE fun.rough_period = 0.5 return fun # The date given in argument is supposed to be given in the given timezone (more natural for a human), # but we need it in UTC. Subtracting the timezone to get it in UTC. start_time = get_timescale().utc( for_date.year, for_date.month, for_date.day, -timezone ) end_time = get_timescale().utc( for_date.year, for_date.month, for_date.day + 1, -timezone ) ephemerides = [] try: for aster in ASTERS: times, risen_info = find_discrete(start_time, end_time, is_risen(aster)) culmination_time, _ = find_maxima(start_time, end_time, get_angle(aster)) rise_time, set_time = None, None culmination_time = ( culmination_time[0] if len(culmination_time) == 1 else None ) for i, time in enumerate(times): time_dt = normalize_datetime( translate_to_timezone(time.utc_datetime(), to_tz=timezone) ) if time_dt is not None and time_dt.day != for_date.day: continue if risen_info[i]: rise_time = time_dt else: set_time = time_dt if culmination_time is not None: culmination_time = normalize_datetime( translate_to_timezone( culmination_time.utc_datetime(), to_tz=timezone, ) ) ephemerides.append( AsterEphemerides(rise_time, culmination_time, set_time, aster=aster) ) except EphemerisRangeError as error: start = translate_to_timezone(error.start_time.utc_datetime(), timezone) end = translate_to_timezone(error.end_time.utc_datetime(), timezone) start = date(start.year, start.month, start.day + 1) end = date(end.year, end.month, end.day - 1) raise OutOfRangeDateError(start, end) from error return ephemerides