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.

model.py 11 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353
  1. #!/usr/bin/env python3
  2. # Kosmorrolib - The Library To Compute Your Ephemerides
  3. # Copyright (C) 2021 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 abc import ABC, abstractmethod
  18. from typing import Union
  19. from datetime import datetime, timezone
  20. from math import asin
  21. from skyfield.api import Topos, Time, Angle
  22. from skyfield.vectorlib import VectorSum as SkfPlanet
  23. from .core import get_skf_objects, get_timescale
  24. from .enum import MoonPhaseType, EventType, ObjectIdentifier, ObjectType
  25. class Serializable(ABC):
  26. @abstractmethod
  27. def serialize(self) -> dict:
  28. pass
  29. class MoonPhase(Serializable):
  30. def __init__(
  31. self,
  32. phase_type: MoonPhaseType,
  33. time: datetime = None,
  34. next_phase_date: datetime = None,
  35. ):
  36. self.phase_type = phase_type
  37. self.time = time
  38. self.next_phase_date = next_phase_date
  39. def __repr__(self):
  40. return "<MoonPhase phase_type=%s time=%s next_phase_date=%s>" % (
  41. self.phase_type,
  42. self.time,
  43. self.next_phase_date,
  44. )
  45. def get_next_phase(self):
  46. """Helper to get the Moon phase that follows the one described by the object.
  47. If the current Moon phase is New Moon or Waxing crescent, the next one will be First Quarter:
  48. >>> moon_phase = MoonPhase(MoonPhaseType.NEW_MOON)
  49. >>> moon_phase.get_next_phase()
  50. <MoonPhaseType.FIRST_QUARTER: 2>
  51. >>> moon_phase = MoonPhase(MoonPhaseType.NEW_MOON)
  52. >>> moon_phase.get_next_phase()
  53. <MoonPhaseType.FIRST_QUARTER: 2>
  54. If the current Moon phase is First Quarter or Waxing gibbous, the next one will be Full Moon:
  55. >>> moon_phase = MoonPhase(MoonPhaseType.FIRST_QUARTER)
  56. >>> moon_phase.get_next_phase()
  57. <MoonPhaseType.FULL_MOON: 4>
  58. >>> moon_phase = MoonPhase(MoonPhaseType.WAXING_GIBBOUS)
  59. >>> moon_phase.get_next_phase()
  60. <MoonPhaseType.FULL_MOON: 4>
  61. If the current Moon phase is Full Moon or Waning gibbous, the next one will be Last Quarter:
  62. >>> moon_phase = MoonPhase(MoonPhaseType.FULL_MOON)
  63. >>> moon_phase.get_next_phase()
  64. <MoonPhaseType.LAST_QUARTER: 6>
  65. >>> moon_phase = MoonPhase(MoonPhaseType.WANING_GIBBOUS)
  66. >>> moon_phase.get_next_phase()
  67. <MoonPhaseType.LAST_QUARTER: 6>
  68. If the current Moon phase is Last Quarter Moon or Waning crescent, the next one will be New Moon:
  69. >>> moon_phase = MoonPhase(MoonPhaseType.LAST_QUARTER)
  70. >>> moon_phase.get_next_phase()
  71. <MoonPhaseType.NEW_MOON: 0>
  72. >>> moon_phase = MoonPhase(MoonPhaseType.WANING_CRESCENT)
  73. >>> moon_phase.get_next_phase()
  74. <MoonPhaseType.NEW_MOON: 0>
  75. """
  76. if self.phase_type in [MoonPhaseType.NEW_MOON, MoonPhaseType.WAXING_CRESCENT]:
  77. return MoonPhaseType.FIRST_QUARTER
  78. if self.phase_type in [
  79. MoonPhaseType.FIRST_QUARTER,
  80. MoonPhaseType.WAXING_GIBBOUS,
  81. ]:
  82. return MoonPhaseType.FULL_MOON
  83. if self.phase_type in [MoonPhaseType.FULL_MOON, MoonPhaseType.WANING_GIBBOUS]:
  84. return MoonPhaseType.LAST_QUARTER
  85. return MoonPhaseType.NEW_MOON
  86. def serialize(self) -> dict:
  87. return {
  88. "phase": self.phase_type.name,
  89. "time": self.time.isoformat() if self.time is not None else None,
  90. "next": {
  91. "phase": self.get_next_phase().name,
  92. "time": self.next_phase_date.isoformat(),
  93. },
  94. }
  95. class Object(Serializable):
  96. """
  97. An astronomical object.
  98. """
  99. def __init__(
  100. self,
  101. identifier: ObjectIdentifier,
  102. skyfield_object: SkfPlanet,
  103. radius: float = None,
  104. ):
  105. """
  106. Initialize an astronomical object
  107. :param ObjectIdentifier identifier: the official name of the object (may be internationalized)
  108. :param str skyfield_object: the object from Skyfield library
  109. :param float radius: the radius (in km) of the object
  110. """
  111. self.identifier = identifier
  112. self.skyfield_object = skyfield_object
  113. self.radius = radius
  114. def __repr__(self):
  115. return "<Object type=%s name=%s />" % (
  116. self.get_type().name,
  117. self.identifier.name,
  118. )
  119. @abstractmethod
  120. def get_type(self) -> ObjectType:
  121. pass
  122. def get_apparent_radius(self, for_date: Union[Time, datetime]) -> Angle:
  123. """Calculate the apparent radius, in degrees, of the object from the given place at a given time.
  124. **Warning:** this is an internal function, not intended for use by end-developers.
  125. For an easier usage, this method accepts datetime and Skyfield's Time objects:
  126. >>> sun = ASTERS[0]
  127. >>> sun.get_apparent_radius(datetime(2021, 6, 9, tzinfo=timezone.utc))
  128. <Angle 00deg 31' 31.6">
  129. >>> sun.get_apparent_radius(get_timescale().utc(2021, 6, 9))
  130. <Angle 00deg 31' 31.6">
  131. Source of the algorithm: https://rhodesmill.org/skyfield/examples.html#what-is-the-angular-diameter-of-a-planet-given-its-radius
  132. :param for_date: the date for which the apparent radius has to be returned
  133. :return: an object representing a Skyfield angle
  134. """
  135. if isinstance(for_date, datetime):
  136. for_date = get_timescale().from_datetime(for_date)
  137. ra, dec, distance = (
  138. EARTH.skyfield_object.at(for_date)
  139. .observe(self.skyfield_object)
  140. .apparent()
  141. .radec()
  142. )
  143. return Angle(radians=asin(self.radius / distance.km) * 2.0)
  144. def serialize(self) -> dict:
  145. """Serialize the given object
  146. >>> planet = Planet(ObjectIdentifier.MARS, "MARS")
  147. >>> planet.serialize()
  148. {'identifier': 'MARS', 'type': 'PLANET', 'radius': None}
  149. """
  150. return {
  151. "identifier": self.identifier.name,
  152. "type": self.get_type().name,
  153. "radius": self.radius,
  154. }
  155. class Star(Object):
  156. def get_type(self) -> ObjectType:
  157. return ObjectType.STAR
  158. class Planet(Object):
  159. def get_type(self) -> ObjectType:
  160. return ObjectType.PLANET
  161. class DwarfPlanet(Planet):
  162. def get_type(self) -> ObjectType:
  163. return ObjectType.DWARF_PLANET
  164. class Satellite(Object):
  165. def get_type(self) -> ObjectType:
  166. return ObjectType.SATELLITE
  167. class Event(Serializable):
  168. def __init__(
  169. self,
  170. event_type: EventType,
  171. objects: [Object],
  172. start_time: datetime,
  173. end_time: Union[datetime, None] = None,
  174. details: {str: any} = None,
  175. ):
  176. self.event_type = event_type
  177. self.objects = objects
  178. self.start_time = start_time
  179. self.end_time = end_time
  180. self.details = details
  181. def __repr__(self):
  182. return "<Event type=%s objects=%s start=%s end=%s details=%s />" % (
  183. self.event_type.name,
  184. self.objects,
  185. self.start_time,
  186. self.end_time,
  187. self.details,
  188. )
  189. def get_description(self, show_details: bool = True) -> str:
  190. description = self.event_type.value % self._get_objects_name()
  191. if show_details and self.details is not None:
  192. description += " ({:s})".format(self.details)
  193. return description
  194. def _get_objects_name(self):
  195. if len(self.objects) == 1:
  196. return self.objects[0].name
  197. return tuple(object.name for object in self.objects)
  198. def serialize(self) -> dict:
  199. return {
  200. "objects": [object.serialize() for object in self.objects],
  201. "EventType": self.event_type.name,
  202. "starts_at": self.start_time.isoformat(),
  203. "ends_at": self.end_time.isoformat() if self.end_time is not None else None,
  204. "details": self.details,
  205. }
  206. class AsterEphemerides(Serializable):
  207. def __init__(
  208. self,
  209. rise_time: Union[datetime, None],
  210. culmination_time: Union[datetime, None],
  211. set_time: Union[datetime, None],
  212. aster: Object,
  213. ):
  214. self.rise_time = rise_time
  215. self.culmination_time = culmination_time
  216. self.set_time = set_time
  217. self.object = aster
  218. def __repr__(self):
  219. return (
  220. "<AsterEphemerides rise_time=%s culmination_time=%s set_time=%s aster=%s>"
  221. % (self.rise_time, self.culmination_time, self.set_time, self.object)
  222. )
  223. def serialize(self) -> dict:
  224. return {
  225. "object": self.object.serialize(),
  226. "rise_time": (
  227. self.rise_time.isoformat() if self.rise_time is not None else None
  228. ),
  229. "culmination_time": (
  230. self.culmination_time.isoformat()
  231. if self.culmination_time is not None
  232. else None
  233. ),
  234. "set_time": (
  235. self.set_time.isoformat() if self.set_time is not None else None
  236. ),
  237. }
  238. EARTH = Planet(ObjectIdentifier.EARTH, get_skf_objects()["EARTH"])
  239. ASTERS = [
  240. Star(ObjectIdentifier.SUN, get_skf_objects()["SUN"], radius=696342),
  241. Satellite(ObjectIdentifier.MOON, get_skf_objects()["MOON"], radius=1737.4),
  242. Planet(ObjectIdentifier.MERCURY, get_skf_objects()["MERCURY"], radius=2439.7),
  243. Planet(ObjectIdentifier.VENUS, get_skf_objects()["VENUS"], radius=6051.8),
  244. Planet(ObjectIdentifier.MARS, get_skf_objects()["MARS"], radius=3396.2),
  245. Planet(
  246. ObjectIdentifier.JUPITER, get_skf_objects()["JUPITER BARYCENTER"], radius=71492
  247. ),
  248. Planet(
  249. ObjectIdentifier.SATURN, get_skf_objects()["SATURN BARYCENTER"], radius=60268
  250. ),
  251. Planet(
  252. ObjectIdentifier.URANUS, get_skf_objects()["URANUS BARYCENTER"], radius=25559
  253. ),
  254. Planet(
  255. ObjectIdentifier.NEPTUNE, get_skf_objects()["NEPTUNE BARYCENTER"], radius=24764
  256. ),
  257. Planet(ObjectIdentifier.PLUTO, get_skf_objects()["PLUTO BARYCENTER"], radius=1185),
  258. ]
  259. def get_aster(identifier: ObjectIdentifier) -> Object:
  260. """Return the aster with the given identifier
  261. >>> get_aster(ObjectIdentifier.SATURN)
  262. <Object type=PLANET name=SATURN />
  263. You can also use it to get the `EARTH` object, even though it has its own constant:
  264. <Object type=PLANET name=EARTH />
  265. """
  266. if identifier == ObjectIdentifier.EARTH:
  267. return EARTH
  268. for aster in ASTERS:
  269. if aster.identifier == identifier:
  270. return aster
  271. class Position:
  272. def __init__(self, latitude: float, longitude: float):
  273. self.latitude = latitude
  274. self.longitude = longitude
  275. self._topos = None
  276. def get_planet_topos(self) -> Topos:
  277. if self._topos is None:
  278. self._topos = EARTH.skyfield_object + Topos(
  279. latitude_degrees=self.latitude, longitude_degrees=self.longitude
  280. )
  281. return self._topos