|
- #!/usr/bin/env python3
-
- # Kosmorro - Compute The Next Ephemerides
- # Copyright (C) 2019 Jérôme Deuchnord <jerome@deuchnord.fr>
- #
- # 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 <https://www.gnu.org/licenses/>.
-
- import datetime
- from typing import Union
-
- from skyfield import almanac
- from skyfield.timelib import Time
- from skyfield.constants import tau
-
- from .data import Object, Position, AsterEphemerides, MoonPhase
- from .core import get_skf_objects, get_timescale, get_iau2000b, ASTERS, MONTHS, skyfield_to_moon_phase
-
- RISEN_ANGLE = -0.8333
-
-
- class EphemeridesComputer:
- def __init__(self, position: Union[Position, None]):
- if position is not None:
- position.observation_planet = get_skf_objects()['earth']
- self.position = position
-
- def get_sun(self, start_time, end_time) -> dict:
- times, is_risen = almanac.find_discrete(start_time,
- end_time,
- almanac.sunrise_sunset(get_skf_objects(), self.position))
-
- sunrise = times[0] if is_risen[0] else times[1]
- sunset = times[1] if not is_risen[1] else times[0]
-
- return {'rise': sunrise, 'set': sunset}
-
- @staticmethod
- def get_moon_phase(year, month, day) -> MoonPhase:
- 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(year, month, day)
- time1 = get_timescale().utc(year, month, day - 10)
- time2 = get_timescale().utc(year, month, day + 10)
-
- times, phase = almanac.find_discrete(time1, time2, moon_phase_at)
-
- return skyfield_to_moon_phase(times, phase, today)
-
- @staticmethod
- def get_asters_ephemerides_for_aster(aster, date: datetime.date, position: Position) -> Object:
- skyfield_aster = get_skf_objects()[aster.skyfield_name]
-
- def get_angle(time: Time) -> float:
- return position.get_planet_topos().at(time).observe(skyfield_aster).apparent().altaz()[0].degrees
-
- def is_risen(time: Time) -> bool:
- return get_angle(time) > RISEN_ANGLE
-
- get_angle.rough_period = 1.0
- is_risen.rough_period = 0.5
-
- start_time = get_timescale().utc(date.year, date.month, date.day)
- end_time = get_timescale().utc(date.year, date.month, date.day, 23, 59, 59)
-
- rise_times, arr = almanac.find_discrete(start_time, end_time, is_risen)
- try:
- culmination_time, _ = almanac._find_maxima(start_time, end_time, get_angle, epsilon=1./3600/24)
- except ValueError:
- culmination_time = None
-
- if len(rise_times) == 2:
- rise_time = rise_times[0 if arr[0] else 1]
- set_time = rise_times[1 if not arr[1] else 0]
- else:
- rise_time = rise_times[0] if arr[0] else None
- set_time = rise_times[0] if not arr[0] else None
-
- culmination_time = culmination_time[0] if culmination_time is not None else None
-
- aster.ephemerides = AsterEphemerides(rise_time, culmination_time, set_time)
- return aster
-
- @staticmethod
- def is_leap_year(year: int) -> bool:
- return (year % 4 == 0 and year % 100 > 0) or (year % 400 == 0)
-
- def compute_ephemerides_for_day(self, year: int, month: int, day: int) -> dict:
- return {'moon_phase': self.get_moon_phase(year, month, day),
- 'details': [self.get_asters_ephemerides_for_aster(aster, datetime.date(year, month, day), self.position)
- for aster in ASTERS] if self.position is not None else []}
-
- def compute_ephemerides_for_month(self, year: int, month: int) -> [dict]:
- if month == 2:
- max_day = 29 if self.is_leap_year(year) else 28
- elif month < 8:
- max_day = 30 if month % 2 == 0 else 31
- else:
- max_day = 31 if month % 2 == 0 else 30
-
- ephemerides = []
-
- for day in range(1, max_day + 1):
- ephemerides.append(self.compute_ephemerides_for_day(year, month, day))
-
- return ephemerides
-
- def compute_ephemerides_for_year(self, year: int) -> [dict]:
- ephemerides = {'seasons': self.get_seasons(year)}
-
- for month in range(0, 12):
- ephemerides[MONTHS[month]] = self.compute_ephemerides_for_month(year, month + 1)
-
- return ephemerides
-
- @staticmethod
- def get_seasons(year: int) -> dict:
- start_time = get_timescale().utc(year, 1, 1)
- end_time = get_timescale().utc(year, 12, 31)
- times, almanac_seasons = almanac.find_discrete(start_time, end_time, almanac.seasons(get_skf_objects()))
-
- seasons = {}
- for time, almanac_season in zip(times, almanac_seasons):
- if almanac_season == 0:
- season = 'MARCH'
- elif almanac_season == 1:
- season = 'JUNE'
- elif almanac_season == 2:
- season = 'SEPTEMBER'
- elif almanac_season == 3:
- season = 'DECEMBER'
- else:
- raise AssertionError
-
- seasons[season] = time.utc_iso()
-
- return seasons
-
- def compute_ephemerides(self, year: int, month: int, day: int):
- if day is not None:
- return self.compute_ephemerides_for_day(year, month, day)
-
- if month is not None:
- return self.compute_ephemerides_for_month(year, month)
-
- return self.compute_ephemerides_for_year(year)
|