Kosmorro
/
lib
mirror of https://github.com/Kosmorro/lib
1
0
Fork 0
Code Issues Releases 16 Activity
A library that computes the ephemerides.
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.
67 Commits
6 Branches
781 KiB
 
 
Tree: 59c990ccf9
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from '59c990ccf9'
${ noResults }
lib/kosmorrolib/ephemerides.py

237 lines
17 KiB
Raw Blame History 

  1. #!/usr/bin/env python3

  2. 

  3. import datetime

  4. from typing import Union

  5. 

  6. from skyfield.searchlib import find_discrete, find_maxima

  7. from skyfield.timelib import Time

  8. from skyfield.constants import tau

  9. from skyfield.errors import EphemerisRangeError

  10. 

  11. from .model import Position, AsterEphemerides, MoonPhase, Object, ASTERS

  12. from .dateutil import translate_to_timezone, normalize_datetime

  13. from .core import get_skf_objects, get_timescale, get_iau2000b

  14. from .enum import MoonPhaseType

  15. from .exceptions import OutOfRangeDateError

  16. 

  17. RISEN_ANGLE = -0.8333

  18. 

  19. 

  20. def _get_skyfield_to_moon_phase(

  21.     times: [Time], vals: [int], now: Time, timezone: int

  22. ) -> Union[MoonPhase, None]:

  23.     tomorrow = get_timescale().utc(

  24.         now.utc_datetime().year, now.utc_datetime().month, now.utc_datetime().day + 1

  25.     )

  26. 

  27.     phases = list(MoonPhaseType)

  28.     current_phase = None

  29.     current_phase_time = None

  30.     next_phase_time = None

  31.     i = 0

  32. 

  33.     if len(times) == 0:

  34.         return None

  35. 

  36.     for i, time in enumerate(times):

  37.         if now.utc_iso() <= time.utc_iso():

  38.             if vals[i] in [0, 2, 4, 6]:

  39.                 if time.utc_datetime() < tomorrow.utc_datetime():

  40.                     current_phase_time = time

  41.                     current_phase = phases[vals[i]]

  42.                 else:

  43.                     i -= 1

  44.                     current_phase_time = None

  45.                     current_phase = phases[vals[i]]

  46.             else:

  47.                 current_phase = phases[vals[i]]

  48. 

  49.             break

  50. 

  51.     for j in range(i + 1, len(times)):

  52.         if vals[j] in [0, 2, 4, 6]:

  53.             next_phase_time = times[j]

  54.             break

  55. 

  56.     return MoonPhase(

  57.         current_phase,

  58.         translate_to_timezone(current_phase_time.utc_datetime(), timezone)

  59.         if current_phase_time is not None

  60.         else None,

  61.         translate_to_timezone(next_phase_time.utc_datetime(), timezone)

  62.         if next_phase_time is not None

  63.         else None,

  64.     )

  65. 

  66. 

  67. def get_moon_phase(

  68.     for_date: datetime.date = datetime.date.today(), timezone: int = 0

  69. ) -> MoonPhase:

  70.     """Calculate and return the moon phase for the given date, adjusted to the given timezone if any.

  71. 

  72.     Get the moon phase for the 27 March, 2021:

  73. 

  74.     >>> get_moon_phase(datetime.date.fromisoformat("2021-03-27"))

  75.     <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 18:48:10.902298+00:00>

  76. 

  77.     Get the moon phase for the 27 March, 2021, in the UTC+2 timezone:

  78. 

  79.     >>> get_moon_phase(datetime.date.fromisoformat("2021-03-27"), timezone=2)

  80.     <MoonPhase phase_type=MoonPhaseType.WAXING_GIBBOUS time=None next_phase_date=2021-03-28 20:48:10.902298+02:00>

  81.     """

  82.     earth = get_skf_objects()["earth"]

  83.     moon = get_skf_objects()["moon"]

  84.     sun = get_skf_objects()["sun"]

  85. 

  86.     def moon_phase_at(time: Time):

  87.         time._nutation_angles = get_iau2000b(time)

  88.         current_earth = earth.at(time)

  89.         _, mlon, _ = current_earth.observe(moon).apparent().ecliptic_latlon("date")

  90.         _, slon, _ = current_earth.observe(sun).apparent().ecliptic_latlon("date")

  91.         return (((mlon.radians - slon.radians) // (tau / 8)) % 8).astype(int)

  92. 

  93.     moon_phase_at.rough_period = 7.0  # one lunar phase per week

  94. 

  95.     today = get_timescale().utc(for_date.year, for_date.month, for_date.day)

  96.     time1 = get_timescale().utc(for_date.year, for_date.month, for_date.day - 10)

  97.     time2 = get_timescale().utc(for_date.year, for_date.month, for_date.day + 10)

  98. 

  99.     try:

  100.         times, phase = find_discrete(time1, time2, moon_phase_at)

  101.     except EphemerisRangeError as error:

  102.         start = translate_to_timezone(error.start_time.utc_datetime(), timezone)

  103.         end = translate_to_timezone(error.end_time.utc_datetime(), timezone)

  104. 

  105.         start = datetime.date(start.year, start.month, start.day) + datetime.timedelta(

  106.             days=12

  107.         )

  108.         end = datetime.date(end.year, end.month, end.day) - datetime.timedelta(days=12)

  109. 

  110.         raise OutOfRangeDateError(start, end) from error

  111. 

  112.     return _get_skyfield_to_moon_phase(times, phase, today, timezone)

  113. 

  114. 

  115. def get_ephemerides(

  116.     position: Position, date: datetime.date = datetime.date.today(), timezone: int = 0

  117. ) -> [AsterEphemerides]:

  118.     """Compute and return the ephemerides for the given position and date, adjusted to the given timezone if any.

  119. 

  120.     Compute the ephemerides for June 9th, 2021:

  121. 

  122.     >>> pos = Position(50.5824, 3.0624)

  123.     >>> get_ephemerides(pos, datetime.date(2021, 6, 9))

  124.     [<AsterEphemerides rise_time=2021-06-09 03:36:00 culmination_time=2021-06-09 11:47:00 set_time=2021-06-09 19:58:00 aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-09 02:59:00 culmination_time=2021-06-09 11:02:00 set_time=2021-06-09 19:16:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-09 04:06:00 culmination_time=2021-06-09 11:58:00 set_time=2021-06-09 19:49:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-06-09 04:52:00 culmination_time=2021-06-09 13:13:00 set_time=2021-06-09 21:34:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-09 06:38:00 culmination_time=2021-06-09 14:40:00 set_time=2021-06-09 22:41:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-09 23:43:00 culmination_time=2021-06-09 04:54:00 set_time=2021-06-09 10:01:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-09 23:02:00 culmination_time=2021-06-09 03:41:00 set_time=2021-06-09 08:16:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-09 01:56:00 culmination_time=2021-06-09 09:18:00 set_time=2021-06-09 16:40:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-09 00:27:00 culmination_time=2021-06-09 06:13:00 set_time=2021-06-09 11:59:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=2021-06-09 22:22:00 culmination_time=2021-06-09 02:32:00 set_time=2021-06-09 06:38:00 aster=<Object type=PLANET name=PLUTO />>]

  125. 

  126.     Compute the ephemerides for June 9th, 2021:

  127. 

  128.     >>> get_ephemerides(pos, datetime.date(2021, 6, 9), timezone=2)

  129.     [<AsterEphemerides rise_time=2021-06-09 05:36:00 culmination_time=2021-06-09 13:47:00 set_time=2021-06-09 21:58:00 aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-09 04:59:00 culmination_time=2021-06-09 13:02:00 set_time=2021-06-09 21:16:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-09 06:06:00 culmination_time=2021-06-09 13:58:00 set_time=2021-06-09 21:49:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-06-09 06:52:00 culmination_time=2021-06-09 15:13:00 set_time=2021-06-09 23:34:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-09 08:38:00 culmination_time=2021-06-09 16:40:00 set_time=2021-06-09 00:44:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-09 01:47:00 culmination_time=2021-06-09 06:54:00 set_time=2021-06-09 12:01:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-09 01:06:00 culmination_time=2021-06-09 05:41:00 set_time=2021-06-09 10:16:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-09 03:56:00 culmination_time=2021-06-09 11:18:00 set_time=2021-06-09 18:40:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-09 02:27:00 culmination_time=2021-06-09 08:13:00 set_time=2021-06-09 13:59:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=2021-06-09 00:26:00 culmination_time=2021-06-09 04:32:00 set_time=2021-06-09 08:38:00 aster=<Object type=PLANET name=PLUTO />>]

  130. 

  131.     If an objet does not rise nor set due to your latitude, then both rise and set will be `None`:

  132. 

  133.     >>> north_pole = Position(70, 20)

  134.     >>> south_pole = Position(-70, 20)

  135.     >>> get_ephemerides(north_pole, datetime.date(2021, 6, 20))

  136.     [<AsterEphemerides rise_time=None culmination_time=2021-06-20 10:42:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-20 14:30:00 culmination_time=2021-06-20 18:44:00 set_time=2021-06-20 22:53:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-20 22:56:00 culmination_time=2021-06-20 09:47:00 set_time=2021-06-20 20:34:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 12:20:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 13:17:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-20 23:06:00 culmination_time=2021-06-20 03:04:00 set_time=2021-06-20 06:58:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-20 23:28:00 culmination_time=2021-06-20 01:48:00 set_time=2021-06-20 04:05:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-20 21:53:00 culmination_time=2021-06-20 07:29:00 set_time=2021-06-20 17:02:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-20 22:51:00 culmination_time=2021-06-20 04:22:00 set_time=2021-06-20 09:50:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 00:40:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]

  137. 

  138.     >>> get_ephemerides(north_pole, datetime.date(2021, 12, 21))

  139.     [<AsterEphemerides rise_time=None culmination_time=2021-12-21 10:38:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 00:04:00 set_time=None aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 11:33:00 set_time=None aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=2021-12-21 11:58:00 culmination_time=2021-12-21 12:33:00 set_time=2021-12-21 13:08:00 aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 08:54:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-12-21 11:07:00 culmination_time=2021-12-21 14:43:00 set_time=2021-12-21 18:19:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-12-21 11:32:00 culmination_time=2021-12-21 13:33:00 set_time=2021-12-21 15:33:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-12-21 09:54:00 culmination_time=2021-12-21 19:13:00 set_time=2021-12-21 04:37:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-12-21 10:49:00 culmination_time=2021-12-21 16:05:00 set_time=2021-12-21 21:21:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-12-21 12:31:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]

  140. 

  141.     >>> get_ephemerides(south_pole, datetime.date(2021, 6, 20))

  142.     [<AsterEphemerides rise_time=None culmination_time=2021-06-20 10:42:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=2021-06-20 11:10:00 culmination_time=2021-06-20 19:06:00 set_time=2021-06-20 01:20:00 aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=2021-06-20 07:47:00 culmination_time=2021-06-20 09:47:00 set_time=2021-06-20 11:48:00 aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 12:20:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=2021-06-20 12:14:00 culmination_time=2021-06-20 13:17:00 set_time=2021-06-20 14:21:00 aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-06-20 18:32:00 culmination_time=2021-06-20 03:04:00 set_time=2021-06-20 11:32:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-06-20 15:20:00 culmination_time=2021-06-20 01:48:00 set_time=2021-06-20 12:12:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-06-20 04:32:00 culmination_time=2021-06-20 07:29:00 set_time=2021-06-20 10:26:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-06-20 21:28:00 culmination_time=2021-06-20 04:22:00 set_time=2021-06-20 11:13:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-06-20 00:40:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]

  143. 

  144.     >>> get_ephemerides(south_pole, datetime.date(2021, 12, 22))

  145.     [<AsterEphemerides rise_time=None culmination_time=2021-12-22 10:39:00 set_time=None aster=<Object type=STAR name=SUN />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 01:01:00 set_time=None aster=<Object type=SATELLITE name=MOON />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 11:35:00 set_time=None aster=<Object type=PLANET name=MERCURY />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 12:27:00 set_time=None aster=<Object type=PLANET name=VENUS />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 08:53:00 set_time=None aster=<Object type=PLANET name=MARS />>, <AsterEphemerides rise_time=2021-12-22 05:52:00 culmination_time=2021-12-22 14:40:00 set_time=2021-12-22 23:26:00 aster=<Object type=PLANET name=JUPITER />>, <AsterEphemerides rise_time=2021-12-22 02:41:00 culmination_time=2021-12-22 13:29:00 set_time=2021-12-22 00:21:00 aster=<Object type=PLANET name=SATURN />>, <AsterEphemerides rise_time=2021-12-22 16:01:00 culmination_time=2021-12-22 19:09:00 set_time=2021-12-22 22:17:00 aster=<Object type=PLANET name=URANUS />>, <AsterEphemerides rise_time=2021-12-22 08:59:00 culmination_time=2021-12-22 16:01:00 set_time=2021-12-22 23:04:00 aster=<Object type=PLANET name=NEPTUNE />>, <AsterEphemerides rise_time=None culmination_time=2021-12-22 12:27:00 set_time=None aster=<Object type=PLANET name=PLUTO />>]

  146.     """

  147.     ephemerides = []

  148. 

  149.     def get_angle(for_aster: Object):

  150.         def fun(time: Time) -> float:

  151.             return (

  152.                 position.get_planet_topos()

  153.                 .at(time)

  154.                 .observe(for_aster.get_skyfield_object())

  155.                 .apparent()

  156.                 .altaz()[0]

  157.                 .degrees

  158.             )

  159. 

  160.         fun.rough_period = 1.0

  161.         return fun

  162. 

  163.     def is_risen(for_aster: Object):

  164.         def fun(time: Time) -> bool:

  165.             return get_angle(for_aster)(time) > RISEN_ANGLE

  166. 

  167.         fun.rough_period = 0.5

  168.         return fun

  169. 

  170.     start_time = get_timescale().utc(date.year, date.month, date.day, -timezone)

  171.     end_time = get_timescale().utc(

  172.         date.year, date.month, date.day, 23 - timezone, 59, 59

  173.     )

  174. 

  175.     try:

  176.         for aster in ASTERS:

  177.             rise_times, arr = find_discrete(start_time, end_time, is_risen(aster))

  178.             try:

  179.                 culmination_time, _ = find_maxima(

  180.                     start_time,

  181.                     end_time,

  182.                     f=get_angle(aster),

  183.                     epsilon=1.0 / 3600 / 24,

  184.                     num=12,

  185.                 )

  186.                 culmination_time = (

  187.                     culmination_time[0] if len(culmination_time) > 0 else None

  188.                 )

  189.             except ValueError:

  190.                 culmination_time = None

  191. 

  192.             rise_time, set_time = None, None

  193. 

  194.             if len(rise_times) == 2:

  195.                 rise_time = rise_times[0 if arr[0] else 1]

  196.                 set_time = rise_times[1 if not arr[1] else 0]

  197.             elif len(rise_times) == 1:

  198.                 rise_time = rise_times[0] if arr[0] else None

  199.                 set_time = rise_times[0] if not arr[0] else None

  200. 

  201.             # Convert the Time instances to Python datetime objects

  202.             if rise_time is not None:

  203.                 rise_time = normalize_datetime(

  204.                     translate_to_timezone(

  205.                         rise_time.utc_datetime().replace(microsecond=0), to_tz=timezone

  206.                     )

  207.                 )

  208. 

  209.             if culmination_time is not None:

  210.                 culmination_time = normalize_datetime(

  211.                     translate_to_timezone(

  212.                         culmination_time.utc_datetime().replace(microsecond=0),

  213.                         to_tz=timezone,

  214.                     )

  215.                 )

  216. 

  217.             if set_time is not None:

  218.                 set_time = normalize_datetime(

  219.                     translate_to_timezone(

  220.                         set_time.utc_datetime().replace(microsecond=0), to_tz=timezone

  221.                     )

  222.                 )

  223. 

  224.             ephemerides.append(

  225.                 AsterEphemerides(rise_time, culmination_time, set_time, aster=aster)

  226.             )

  227.     except EphemerisRangeError as error:

  228.         start = translate_to_timezone(error.start_time.utc_datetime(), timezone)

  229.         end = translate_to_timezone(error.end_time.utc_datetime(), timezone)

  230. 

  231.         start = datetime.date(start.year, start.month, start.day + 1)

  232.         end = datetime.date(end.year, end.month, end.day - 1)

  233. 

  234.         raise OutOfRangeDateError(start, end) from error

  235. 

  236.     return ephemerides