mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-19 14:27:52 -04:00
Code Issues Packages Projects Releases Wiki Activity

Reworking Python body codes as enumerated type Body.

Browse Source 
This will help documentation generator pydown.py organize
body codes and other similar enumerations together as classes.
This commit is contained in:
Don Cross
2019-07-08 15:37:56 -04:00
parent 700d3d7870
commit b3e6f185b6
6 changed files with 275 additions and 385 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
source/python/.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
*.pyc

143
source/python/README.md
View File

@@ -1,143 +0,0 @@
# astronomy
Astronomy Engine by Don Cross
See the GitHub project page for full documentation, examples,
and other information:
https://github.com/cosinekitty/astronomy
## BodyCode
```python
BodyCode(name)
```
Finds the integer body code given the name of a body.
Parameters
----------
name: str
The common English name of a supported celestial body.
Returns
-------
int
If `name` is a valid body name, returns the integer value
of the body code associated with that body.
Otherwise, returns `BODY_INVALID`.
Example
-------
>>> astronomy.BodyCode('Mars')
3
## Time
```python
Time(self, ut)
```
Represents a date and time used for performing astronomy calculations.
All calculations performed by Astronomy Engine are based on
dates and times represented by `Time` objects.
Parameters
----------
ut : float
UT1/UTC number of days since noon on January 1, 2000.
See the `ut` attribute of this class for more details.
Attributes
----------
ut : float
The floating point number of days of Universal Time since noon UTC January 1, 2000.
Astronomy Engine approximates UTC and UT1 as being the same thing, although they are
not exactly equivalent; UTC and UT1 can disagree by up to 0.9 seconds.
This approximation is sufficient for the accuracy requirements of Astronomy Engine.
Universal Time Coordinate (UTC) is the international standard for legal and civil
timekeeping and replaces the older Greenwich Mean Time (GMT) standard.
UTC is kept in sync with unpredictable observed changes in the Earth's rotation
by occasionally adding leap seconds as needed.
UT1 is an idealized time scale based on observed rotation of the Earth, which
gradually slows down in an unpredictable way over time, due to tidal drag by the Moon and Sun,
large scale weather events like hurricanes, and internal seismic and convection effects.
Conceptually, UT1 drifts from atomic time continuously and erratically, whereas UTC
is adjusted by a scheduled whole number of leap seconds as needed.
The value in `ut` is appropriate for any calculation involving the Earth's rotation,
such as calculating rise/set times, culumination, and anything involving apparent
sidereal time.
Before the era of atomic timekeeping, days based on the Earth's rotation
were often known as <i>mean solar days</i>.
tt : float
Terrestrial Time days since noon on January 1, 2000.
Terrestrial Time is an atomic time scale defined as a number of days since noon on January 1, 2000.
In this system, days are not based on Earth rotations, but instead by
the number of elapsed [SI seconds](https://physics.nist.gov/cuu/Units/second.html)
divided by 86400. Unlike `ut`, `tt` increases uniformly without adjustments
for changes in the Earth's rotation.
The value in `tt` is used for calculations of movements not involving the Earth's rotation,
such as the orbits of planets around the Sun, or the Moon around the Earth.
Historically, Terrestrial Time has also been known by the term <i>Ephemeris Time</i> (ET).
### Make
```python
Time.Make(year, month, day, hour, minute, second)
```
Creates a `Time` object from a UTC calendar date and time.
Parameters
----------
year : int
The UTC 4-digit year value, e.g. 2019.
month : int
The UTC month in the range 1..12.
day : int
The UTC day of the month, in the range 1..31.
hour : int
The UTC hour, in the range 0..23.
minute : int
The UTC minute, in the range 0..59.
second : float
The real-valued UTC second, in the range [0, 60).
Returns
-------
Time
## Observer
```python
Observer(self, latitude, longitude, height=0)
```
Represents the geographic location of an observer on the surface of the Earth.
:param latitude: Geographic latitude in degrees north of the equator.
:param longitude: Geographic longitude in degrees east of the prime meridian at Greenwich, England.
:param height: Elevation above sea level in meters.
## GeoMoon
```python
GeoMoon(time)
```
Calculates the geocentric position of the Moon at a given time.
Given a time of observation, calculates the Moon's position as a vector.
The vector gives the location of the Moon's center relative to the Earth's center
with x-, y-, and z-components measured in astronomical units.
This algorithm is based on Nautical Almanac Office's <i>Improved Lunar Ephemeris</i> of 1954,
which in turn derives from E. W. Brown's lunar theories from the early twentieth century.
It is adapted from Turbo Pascal code from the book
[Astronomy on the Personal Computer](https://www.springer.com/us/book/9783540672210)
by Montenbruck and Pfleger.
Parameters
----------
time : Time
The date and time for which to calculate the Moon's position.
Returns
-------
Vector
The Moon's position as a vector in J2000 Cartesian equatorial coordinates.

142
source/python/astronomy.py
View File

@@ -33,6 +33,7 @@ https://github.com/cosinekitty/astronomy
import math
import datetime
from enum import IntEnum, unique
_PI2 = 2.0 * math.pi
_EPOCH = datetime.datetime(2000, 1, 1, 12)
@@ -84,40 +85,43 @@ class Vector:
def Length(self):
return math.sqrt(self.x**2 + self.y**2 + self.z**2)
BODY_INVALID = -1
@unique
class Body(IntEnum):
Invalid = -1
"""A placeholder value for an unknown, undefined, or invalid body."""
BODY_MERCURY = 0
"""The body code for the planet Mercury."""
Mercury = 0
"""The body code for the planet Mercury."""
BODY_VENUS = 1
"""The body code for the planet Venus."""
Venus = 1
"""The body code for the planet Venus."""
BODY_EARTH = 2
"""The body code for the planet Earth."""
Earth = 2
"""The body code for the planet Earth."""
BODY_MARS = 3
"""The body code for the planet Mars."""
Mars = 3
"""The body code for the planet Mars."""
BODY_JUPITER = 4
"""The body code for the planet Jupiter."""
Jupiter = 4
"""The body code for the planet Jupiter."""
BODY_SATURN = 5
"""The body code for the planet Saturn."""
Saturn = 5
"""The body code for the planet Saturn."""
BODY_URANUS = 6
"""The body code for the planet Uranus."""
Uranus = 6
"""The body code for the planet Uranus."""
BODY_NEPTUNE = 7
"""The body code for the planet Neptune."""
Neptune = 7
"""The body code for the planet Neptune."""
BODY_PLUTO = 8
"""The body code for the planet Pluto."""
Pluto = 8
"""The body code for the planet Pluto."""
BODY_SUN = 9
"""The body code for the Sun."""
Sun = 9
"""The body code for the Sun."""
BODY_MOON = 10
"""The body code for the Moon."""
Moon = 10
"""The body code for the Moon."""
BodyName = [
'Mercury',
@@ -134,12 +138,12 @@ BodyName = [
]
"""The English names of the supported celestial bodies.
The list `BodyName` is indexed using one of the `BODY_...` constants.
The list `BodyName` is indexed using one of the `Body....` constants.
Example
-------
>>> astronomy.BodyName[astronomy.BODY_JUPITER]
>>> astronomy.BodyName[astronomy.Body.Jupiter]
'Jupiter'
"""
@@ -157,7 +161,7 @@ def BodyCode(name):
int
If `name` is a valid body name, returns the integer value
of the body code associated with that body.
Otherwise, returns `BODY_INVALID`.
Otherwise, returns `Body.Invalid`.
Example
-------
@@ -167,11 +171,11 @@ def BodyCode(name):
"""
if name not in BodyName:
return BODY_INVALID
return BodyName.index(name)
return Body.Invalid
return Body[name]
def _IsSuperiorPlanet(body):
return body in [BODY_MARS, BODY_JUPITER, BODY_SATURN, BODY_URANUS, BODY_NEPTUNE, BODY_PLUTO]
return body in [Body.Mars, Body.Jupiter, Body.Saturn, Body.Uranus, Body.Neptune, Body.Pluto]
_PlanetOrbitalPeriod = [
87.969,
@@ -210,11 +214,11 @@ class NoConvergeError(Error):
Error.__init__(self, 'Numeric solver did not converge - please report issue at https://github.com/cosinekitty/astronomy/issues')
def _SynodicPeriod(body):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
if body < 0 or body >= len(_PlanetOrbitalPeriod):
raise InvalidBodyError()
if body == BODY_MOON:
if body == Body.Moon:
return _MEAN_SYNODIC_MONTH
return abs(_EARTH_ORBITAL_PERIOD / (_EARTH_ORBITAL_PERIOD/_PlanetOrbitalPeriod[body] - 1.0))
@@ -2730,7 +2734,7 @@ def _CalcVsop(model, time):
return Vector(vx, vy, vz, time)
def _CalcEarth(time):
return _CalcVsop(_vsop[BODY_EARTH], time)
return _CalcVsop(_vsop[Body.Earth], time)
# END VSOP
#----------------------------------------------------------------------------
@@ -3025,16 +3029,16 @@ def Search(func, context, t1, t2, dt_tolerance_seconds):
#----------------------------------------------------------------------------
def HelioVector(body, time):
if body == BODY_PLUTO:
if body == Body.Pluto:
return _CalcChebyshev(_pluto, time)
if 0 <= body <= len(_vsop):
return _CalcVsop(_vsop[body], time)
if body == BODY_SUN:
if body == Body.Sun:
return Vector(0.0, 0.0, 0.0, time)
if body == BODY_MOON:
if body == Body.Moon:
e = _CalcEarth(time)
m = GeoMoon(time)
return Vector(e.x+m.x, e.y+m.y, e.z+m.z, time)
@@ -3043,10 +3047,10 @@ def HelioVector(body, time):
def GeoVector(body, time, aberration):
if body == BODY_MOON:
if body == Body.Moon:
return GeoMoon(time)
if body == BODY_EARTH:
if body == Body.Earth:
return Vector(0.0, 0.0, 0.0, time)
if not aberration:
@@ -3072,7 +3076,7 @@ def GeoVector(body, time, aberration):
earth = _CalcEarth(ltime)
geo = Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time)
if body == BODY_SUN:
if body == Body.Sun:
# The Sun's heliocentric coordinates are always (0,0,0). No need to correct.
return geo
@@ -3242,23 +3246,23 @@ def Ecliptic(equ):
return _RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000)
def EclipticLongitude(body, time):
if body == BODY_SUN:
if body == Body.Sun:
raise InvalidBodyError()
hv = HelioVector(body, time)
eclip = Ecliptic(hv)
return eclip.elon
def AngleFromSun(body, time):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
sv = GeoVector(BODY_SUN, time, True)
sv = GeoVector(Body.Sun, time, True)
bv = GeoVector(body, time, True)
return _AngleBetween(sv, bv)
def LongitudeFromSun(body, time):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
sv = GeoVector(BODY_SUN, time, True)
sv = GeoVector(Body.Sun, time, True)
se = Ecliptic(sv)
bv = GeoVector(body, time, True)
be = Ecliptic(bv)
@@ -3284,14 +3288,14 @@ def Elongation(body, time):
def _rlon_offset(body, time, direction, targetRelLon):
plon = EclipticLongitude(body, time)
elon = EclipticLongitude(BODY_EARTH, time)
elon = EclipticLongitude(Body.Earth, time)
diff = direction * (elon - plon)
return _LongitudeOffset(diff - targetRelLon)
def SearchRelativeLongitude(body, targetRelLon, startTime):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
if body == BODY_MOON or body == BODY_SUN:
if body == Body.Moon or body == Body.Sun:
raise InvalidBodyError()
syn = _SynodicPeriod(body)
direction = +1 if _IsSuperiorPlanet(body) else -1
@@ -3331,10 +3335,10 @@ def _neg_elong_slope(body, time):
return (e1 - e2)/dt
def SearchMaxElongation(body, startTime):
if body == BODY_MERCURY:
if body == Body.Mercury:
s1 = 50.0
s2 = 85.0
elif body == BODY_VENUS:
elif body == Body.Venus:
s1 = 40.0
s2 = 50.0
else:
@@ -3343,7 +3347,7 @@ def SearchMaxElongation(body, startTime):
iter = 1
while iter <= 2:
plon = EclipticLongitude(body, startTime)
elon = EclipticLongitude(BODY_EARTH, startTime)
elon = EclipticLongitude(Body.Earth, startTime)
rlon = _LongitudeOffset(plon - elon) # clamp to (-180, +180]
# The slope function is not well-behaved when rlon is near 0 degrees or 180 degrees
@@ -3421,7 +3425,7 @@ def SearchSunLongitude(targetLon, startTime, limitDays):
return Search(_sun_offset, targetLon, startTime, t2, 1.0)
def MoonPhase(time):
return LongitudeFromSun(BODY_MOON, time)
return LongitudeFromSun(Body.Moon, time)
def _moon_offset(targetLon, time):
angle = MoonPhase(time)
@@ -3524,22 +3528,22 @@ def _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time):
def _VisualMagnitude(body, phase, helio_dist, geo_dist):
# For Mercury and Venus, see: https://iopscience.iop.org/article/10.1086/430212
c0 = c1 = c2 = c3 = 0
if body == BODY_MERCURY:
if body == Body.Mercury:
c0 = -0.60; c1 = +4.98; c2 = -4.88; c3 = +3.02
elif body == BODY_VENUS:
elif body == Body.Venus:
if phase < 163.6:
c0 = -4.47; c1 = +1.03; c2 = +0.57; c3 = +0.13
else:
c0 = +0.98; c1 = -1.02
elif body == BODY_MARS:
elif body == Body.Mars:
c0 = -1.52; c1 = +1.60
elif body == BODY_JUPITER:
elif body == Body.Jupiter:
c0 = -9.40; c1 = +0.50
elif body == BODY_URANUS:
elif body == Body.Uranus:
c0 = -7.19; c1 = +0.25
elif body == BODY_NEPTUNE:
elif body == Body.Neptune:
c0 = -6.87
elif body == BODY_PLUTO:
elif body == Body.Pluto:
c0 = -1.00; c1 = +4.00
else:
raise InvalidBodyError()
@@ -3550,15 +3554,15 @@ def _VisualMagnitude(body, phase, helio_dist, geo_dist):
return mag
def Illumination(body, time):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
earth = _CalcEarth(time)
if body == BODY_SUN:
if body == Body.Sun:
gc = Vector(-earth.x, -earth.y, -earth.z, time)
hc = Vector(0.0, 0.0, 0.0, time)
phase = 0.0 # placeholder value; the Sun does not have a phase angle.
else:
if body == BODY_MOON:
if body == Body.Moon:
# For extra numeric precision, use geocentric moon formula directly.
gc = GeoMoon(time)
hc = Vector(earth.x + gc.x, earth.y + gc.y, earth.z + gc.z, time)
@@ -3571,11 +3575,11 @@ def Illumination(body, time):
geo_dist = gc.Length() # distance from body to center of Earth
helio_dist = hc.Length() # distance from body to center of Sun
ring_tilt = None # only reported for Saturn
if body == BODY_SUN:
if body == Body.Sun:
mag = -0.17 + 5.0*math.log10(geo_dist / _AU_PER_PARSEC)
elif body == BODY_MOON:
elif body == Body.Moon:
mag = _MoonMagnitude(phase, helio_dist, geo_dist)
elif body == BODY_SATURN:
elif body == Body.Saturn:
mag, ring_tilt = _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time)
else:
mag = _VisualMagnitude(body, phase, helio_dist, geo_dist)
@@ -3598,7 +3602,7 @@ def SearchPeakMagnitude(body, startTime):
# s1 and s2 are relative longitudes within which peak magnitude of Venus can occur.
s1 = 10.0
s2 = 30.0
if body != BODY_VENUS:
if body != Body.Venus:
raise InvalidBodyError()
iter = 1
@@ -3606,7 +3610,7 @@ def SearchPeakMagnitude(body, startTime):
# Find current heliocentric relative longitude between the
# inferior planet and the Earth.
plon = EclipticLongitude(body, startTime)
elon = EclipticLongitude(BODY_EARTH, startTime)
elon = EclipticLongitude(Body.Earth, startTime)
rlon = _LongitudeOffset(plon - elon)
# The slope function is not well-behaved when rlon is near 0 degrees or 180 degrees
# because there is a cusp there that causes a discontinuity in the derivative.
@@ -3680,7 +3684,7 @@ class HourAngleEvent:
self.hor = hor
def SearchHourAngle(body, observer, hourAngle, startTime):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
if hourAngle < 0.0 or hourAngle >= 24.0:
@@ -3747,11 +3751,11 @@ def _peak_altitude(context, time):
return context.direction * (alt + _REFRACTION_NEAR_HORIZON)
def SearchRiseSet(body, observer, direction, startTime, limitDays):
if body == BODY_EARTH:
if body == Body.Earth:
raise EarthNotAllowedError()
elif body == BODY_SUN:
elif body == Body.Sun:
body_radius = _SUN_RADIUS_AU
elif body == BODY_MOON:
elif body == Body.Moon:
body_radius = _MOON_RADIUS_AU
else:
body_radius = 0.0