Coordinates examples¶
Let’s define a small helper function:
def print_me(msg, val):
print("{}: {}".format(msg, val))
It follows a series of important parameters related to the angle between Earth’s rotation axis and the ecliptic.
The mean angle between Earth rotation axis and ecliptic axis is the mean obliquity:
e0 = Earth.mean_obliquity(1987, 4, 10) print_me("Mean obliquity for 1987/4/10", e0.dms_str()) # Mean obliquity for 1987/4/10: 23d 26' 27.4066466278''
If we take into account the nutation effect on the obliquity, we get the true obliquity:
epsilon = Earth.true_obliquity(1987, 4, 10) print_me("True obliquity for 1987/4/10", epsilon.dms_str()) # True obliquity for 1987/4/10: 23d 26' 36.8491882378'' epsilon = Earth.true_obliquity(2018, 7, 29) print_me("True obliquity for 2018/7/29", epsilon.dms_str()) # True obliquity for 2018/7/29: 23d 26' 7.21570241139''
The nutation effect is separated in two components: One parallel to the ecliptic (nutation in longitude) and other perpendicular to the ecliptic (nutation in obliquity):
dpsi = Earth.nutation_longitude(1987, 4, 10) print_me("Nutation in longitude for 1987/4/10", dpsi.dms_str(n_dec=3)) # Nutation in longitude for 1987/4/10: -3.788'' depsilon = Earth.nutation_obliquity(1987, 4, 10) print_me("Nutation in obliquity for 1987/4/10", depsilon.dms_str(n_dec=3)) # Nutation in obliquity for 1987/4/10: 9.443''
We can compute the effects of precession on the equatorial coordinates of a given star, taking also into account its proper motion:
start_epoch = JDE2000 final_epoch = Epoch(2028, 11, 13.19) alpha0 = Angle(2, 44, 11.986, ra=True) delta0 = Angle(49, 13, 42.48) print_me("Initial right ascension", alpha0.ra_str(n_dec=3)) # Initial right ascension: 2h 44' 11.986'' print_me("Initial declination", delta0.dms_str(n_dec=2)) # Initial declination: 49d 13' 42.48'' pm_ra = Angle(0, 0, 0.03425, ra=True) pm_dec = Angle(0, 0, -0.0895) alpha, delta = Earth.precession_equatorial(start_epoch, final_epoch, alpha0, delta0, pm_ra, pm_dec) print_me("Final right ascension", alpha.ra_str(n_dec=3)) # Final right ascension: 2h 46' 11.331'' print_me("Final declination", delta.dms_str(n_dec=2)) # Final declination: 49d 20' 54.54''
Something similar can also be done with the ecliptical coordinates:
start_epoch = JDE2000
final_epoch = Epoch(-214, 6, 30.0)
lon0 = Angle(149.48194)
lat0 = Angle(1.76549)
print_me("Initial ecliptical longitude", round(lon0(), 5))
# Initial ecliptical longitude: 149.48194
print_me("Initial ecliptical latitude", round(lat0(), 5))
# Initial ecliptical latitude: 1.76549
lon, lat = Earth.precession_ecliptical(start_epoch, final_epoch, lon0, lat0)
print_me("Final ecliptical longitude", round(lon(), 3))
# Final ecliptical longitude: 118.704
print_me("Final ecliptical latitude", round(lat(), 3))
# Final ecliptical latitude: 1.615
Additionally, module Coordinates
provides a function to compute the true movement of a star through the sky relative to the Sun:
ra = Angle(6, 45, 8.871, ra=True)
dec = Angle(-16.716108)
pm_ra = Angle(0, 0, -0.03847, ra=True)
pm_dec = Angle(0, 0, -1.2053)
dist = 2.64
vel = -7.6
alpha, delta = motion_in_space(ra, dec, dist, vel, pm_ra, pm_dec, -4000.0)
print(alpha.ra_str(False, 2))
# 6:47:39.91
print(delta.dms_str(False, 1))
# -15:23:30.6
This module Coordinates
also provides a series of functions to convert between equatorial, ecliptical, horizontal and galactic coordinates.
Equatorial to ecliptical coordinates:
ra = Angle(7, 45, 18.946, ra=True) dec = Angle(28, 1, 34.26) epsilon = Angle(23.4392911) lon, lat = equatorial2ecliptical(ra, dec, epsilon) print_me("Equatorial to ecliptical. Longitude", round(lon(), 5)) # Equatorial to ecliptical. Longitude: 113.21563 print_me("Equatorial to ecliptical. Latitude", round(lat(), 5)) # Equatorial to ecliptical. Latitude: 6.68417
Ecliptical to equatorial coordinates:
lon = Angle(113.21563) lat = Angle(6.68417) epsilon = Angle(23.4392911) ra, dec = ecliptical2equatorial(lon, lat, epsilon) print_me("Ecliptical to equatorial. Right ascension", ra.ra_str(n_dec=3)) # Ecliptical to equatorial. Right ascension: 7h 45' 18.946'' print_me("Ecliptical to equatorial. Declination", dec.dms_str(n_dec=2)) # Ecliptical to equatorial. Declination: 28d 1' 34.26''
Equatorial to horizontal coordinates:
lon = Angle(77, 3, 56) lat = Angle(38, 55, 17) ra = Angle(23, 9, 16.641, ra=True) dec = Angle(-6, 43, 11.61) theta0 = Angle(8, 34, 57.0896, ra=True) eps = Angle(23, 26, 36.87) # Compute correction to convert from mean to apparent sidereal time delta = Angle(0, 0, ((-3.868*cos(eps.rad()))/15.0), ra=True) theta0 += delta h = theta0 - lon - ra azi, ele = equatorial2horizontal(h, dec, lat) print_me("Equatorial to horizontal: Azimuth", round(azi, 3)) # Equatorial to horizontal: Azimuth: 68.034 print_me("Equatorial to horizontal: Elevation", round(ele, 3)) # Equatorial to horizontal: Elevation: 15.125
Horizontal to equatorial coordinates:
azi = Angle(68.0337) ele = Angle(15.1249) lat = Angle(38, 55, 17) h, dec = horizontal2equatorial(azi, ele, lat) print_me("Horizontal to equatorial. Hour angle", round(h, 4)) # Horizontal to equatorial. Hour angle: 64.3521 print_me("Horizontal to equatorial. Declination", dec.dms_str(n_dec=0)) # Horizontal to equatorial. Declination: -6d 43' 12.0''
Equatorial to galactic coordinates:
ra = Angle(17, 48, 59.74, ra=True) dec = Angle(-14, 43, 8.2) lon, lat = equatorial2galactic(ra, dec) print_me("Equatorial to galactic. Longitude", round(lon, 4)) # Equatorial to galactic. Longitude: 12.9593 print_me("Equatorial to galactic. Latitude", round(lat, 4)) # Equatorial to galactic. Latitude: 6.0463
Galactic to equatorial coordinates:
lon = Angle(12.9593) lat = Angle(6.0463) ra, dec = galactic2equatorial(lon, lat) print_me("Galactic to equatorial. Right ascension", ra.ra_str(n_dec=1)) # Galactic to equatorial. Right ascension: 17h 48' 59.7'' print_me("Galactic to equatorial. Declination", dec.dms_str(n_dec=0)) # Galactic to equatorial. Declination: -14d 43' 8.0''
In addition, there is a function to compute the ecliptic longitudes of the two points of the ecliptic which are on the horizon, as well as the angle between the ecliptic and the horizon:
sidereal_time = Angle(5.0, ra=True)
lat = Angle(51.0)
epsilon = Angle(23.44)
lon1, lon2, i = ecliptic_horizon(sidereal_time, lat, epsilon)
print_me("Longitude of ecliptic point #1 on the horizon", lon1.dms_str(n_dec=1))
# Longitude of ecliptic point #1 on the horizon: 169d 21' 29.9''
print_me("Longitude of ecliptic point #2 on the horizon", lon2.dms_str(n_dec=1))
# Longitude of ecliptic point #2 on the horizon: 349d 21' 29.9''
print_me("Angle between the ecliptic and the horizon", round(i, 0))
# Angle between the ecliptic and the horizon: 62.0
Also, it is possible to compute the angle of the diurnal path of a celestial body relative to the horizon at the time of rising and setting:
dec = Angle(23.44)
lat = Angle(40.0)
j = diurnal_path_horizon(dec, lat)
print_me("Diurnal path vs. horizon angle at time of rising and setting", j.dms_str(n_dec=1))
# Diurnal path vs. horizon angle at time of rising and setting: 45d 31' 28.4''
The times (in hours of the day) of rising, transit and setting of a given celestial body can be computed with the appropriate function:
longitude = Angle(71, 5, 0.0)
latitude = Angle(42, 20, 0.0)
alpha1 = Angle(2, 42, 43.25, ra=True)
delta1 = Angle(18, 2, 51.4)
alpha2 = Angle(2, 46, 55.51, ra=True)
delta2 = Angle(18, 26, 27.3)
alpha3 = Angle(2, 51, 7.69, ra=True)
delta3 = Angle(18, 49, 38.7)
h0 = Angle(-0.5667)
delta_t = 56.0
theta0 = Angle(11, 50, 58.1, ra=True)
rising, transit, setting = times_rise_transit_set(longitude, latitude,alpha1, delta1, \
alpha2, delta2, alpha3, delta3, h0, \
delta_t, theta0)
print_me("Time of rising (hours of day)", round(rising, 4))
# Time of rising (hours of day): 12.4238
print_me("Time of transit (hours of day)", round(transit, 3))
# Time of transit (hours of day): 19.675
print_me("Time of setting (hours of day, next day)", round(setting, 3))
# Time of setting (hours of day, next day): 2.911
The air in the atmosphere introduces an error in the elevation due to the refraction. We can compute the true (airless) elevation from the apparent elevation, and viceversa.
Apparent elevation to true (airless) elevation:
apparent_elevation = Angle(0, 30, 0.0) true_elevation = refraction_apparent2true(apparent_elevation) print_me("True elevation for an apparent elevation of 30'", true_elevation.dms_str(n_dec=1)) # True elevation for an apparent elevation of 30': 1' 14.7''
True elevation to apparent elevation:
true_elevation = Angle(0, 33, 14.76) apparent_elevation = refraction_true2apparent(true_elevation) print_me("Apparent elevation for a true elevation of 33' 14.76''", apparent_elevation.dms_str(n_dec=2)) # Apparent elevation for a true elevation of 33' 14.76'': 57' 51.96''
This module provides a function to compute the angular separation between two celestial bodies:
alpha1 = Angle(14, 15, 39.7, ra=True)
delta1 = Angle(19, 10, 57.0)
alpha2 = Angle(13, 25, 11.6, ra=True)
delta2 = Angle(-11, 9, 41.0)
sep_ang = angular_separation(alpha1, delta1, alpha2, delta2)
print_me("Angular separation between two given celestial bodies, in degrees",
round(sep_ang, 3))
# Angular separation between two given celestial bodies, in degrees: 32.793
We can compute the minimum angular separation achieved between two celestial objects. For that, we must provide the positions at three equidistant epochs:
# EPOCH: Sep 13th, 1978, 0h TT:
alpha1_1 = Angle(10, 29, 44.27, ra=True)
delta1_1 = Angle(11, 2, 5.9)
alpha2_1 = Angle(10, 33, 29.64, ra=True)
delta2_1 = Angle(10, 40, 13.2)
# EPOCH: Sep 14th, 1978, 0h TT:
alpha1_2 = Angle(10, 36, 19.63, ra=True)
delta1_2 = Angle(10, 29, 51.7)
alpha2_2 = Angle(10, 33, 57.97, ra=True)
delta2_2 = Angle(10, 37, 33.4)
# EPOCH: Sep 15th, 1978, 0h TT:
alpha1_3 = Angle(10, 43, 1.75, ra=True)
delta1_3 = Angle(9, 55, 16.7)
alpha2_3 = Angle(10, 34, 26.22, ra=True)
delta2_3 = Angle(10, 34, 53.9)
a = minimum_angular_separation(alpha1_1, delta1_1, alpha1_2, delta1_2,
alpha1_3, delta1_3, alpha2_1, delta2_1,
alpha2_2, delta2_2, alpha2_3, delta2_3)
print_me("Minimum angular separation, epoch fraction", round(a[0], 6))
# Minimum angular separation, epoch fraction: -0.370726
# NOTE: Given that 'n' is negative, and Sep 14th is the middle epoch (n=0),
# then the minimum angular separation is achieved on Sep 13th, specifically
# at: 1.0 - 0.370726 = 0.629274 => Sep 13.629274 = Sep 13th, 15h 6' 9''
print_me("Minimum angular separation", a[1].dms_str(n_dec=0))
# Minimum angular separation: 3' 44.0''
There is a function to compute the relative position angle P of a body with respect to another body. In this example, given that the two bodies have the same right ascension, then the relative position angle between them must be 0 (or 180):
alpha1 = Angle(14, 15, 39.7, ra=True)
delta1 = Angle(19, 10, 57.0)
alpha2 = Angle(14, 15, 39.7, ra=True) # Same as alpha1
delta2 = Angle(-11, 9, 41.0)
pos_ang = relative_position_angle(alpha1, delta1, alpha2, delta2)
print_me("Relative position angle", round(pos_ang, 1))
# Relative position angle: 0.0
Planetary conjunctions may be computed with the appropriate function:
alpha1_1 = Angle(10, 24, 30.125, ra=True)
delta1_1 = Angle(6, 26, 32.05)
alpha1_2 = Angle(10, 25, 0.342, ra=True)
delta1_2 = Angle(6, 10, 57.72)
alpha1_3 = Angle(10, 25, 12.515, ra=True)
delta1_3 = Angle(5, 57, 33.08)
alpha1_4 = Angle(10, 25, 6.235, ra=True)
delta1_4 = Angle(5, 46, 27.07)
alpha1_5 = Angle(10, 24, 41.185, ra=True)
delta1_5 = Angle(5, 37, 48.45)
alpha2_1 = Angle(10, 27, 27.175, ra=True)
delta2_1 = Angle(4, 4, 41.83)
alpha2_2 = Angle(10, 26, 32.410, ra=True)
delta2_2 = Angle(3, 55, 54.66)
alpha2_3 = Angle(10, 25, 29.042, ra=True)
delta2_3 = Angle(3, 48, 3.51)
alpha2_4 = Angle(10, 24, 17.191, ra=True)
delta2_4 = Angle(3, 41, 10.25)
alpha2_5 = Angle(10, 22, 57.024, ra=True)
delta2_5 = Angle(3, 35, 16.61)
alpha1_list = [alpha1_1, alpha1_2, alpha1_3, alpha1_4, alpha1_5]
delta1_list = [delta1_1, delta1_2, delta1_3, delta1_4, delta1_5]
alpha2_list = [alpha2_1, alpha2_2, alpha2_3, alpha2_4, alpha2_5]
delta2_list = [delta2_1, delta2_2, delta2_3, delta2_4, delta2_5]
pc = planetary_conjunction(alpha1_list, delta1_list, alpha2_list, delta2_list)
print_me("Epoch fraction 'n' for planetary conjunction", round(pc[0], 5))
# Epoch fraction 'n' for planetary conjunction: 0.23797
print_me("Difference in declination at conjunction", pc[1].dms_str(n_dec=1))
# Difference in declination at conjunction: 2d 8' 21.8''
If the planetary conjunction is with a star, it is a little bit simpler:
alpha_1 = Angle(15, 3, 51.937, ra=True)
delta_1 = Angle(-8, 57, 34.51)
alpha_2 = Angle(15, 9, 57.327, ra=True)
delta_2 = Angle(-9, 9, 3.88)
alpha_3 = Angle(15, 15, 37.898, ra=True)
delta_3 = Angle(-9, 17, 37.94)
alpha_4 = Angle(15, 20, 50.632, ra=True)
delta_4 = Angle(-9, 23, 16.25)
alpha_5 = Angle(15, 25, 32.695, ra=True)
delta_5 = Angle(-9, 26, 1.01)
alpha_star = Angle(15, 17, 0.446, ra=True)
delta_star = Angle(-9, 22, 58.47)
alpha_list = [alpha_1, alpha_2, alpha_3, alpha_4, alpha_5]
delta_list = [delta_1, delta_2, delta_3, delta_4, delta_5]
pc = planet_star_conjunction(alpha_list, delta_list, alpha_star, delta_star)
print_me("Epoch fraction 'n' for planetary conjunction with star", round(pc[0], 4))
# Epoch fraction 'n' for planetary conjunction with star: 0.2551
print_me("Difference in declination with star at conjunction", pc[1].dms_str(n_dec=0))
# Difference in declination with star at conjunction: 3' 38.0''
It is possible to compute when a planet and two other stars will be in a straight line:
alpha_1 = Angle(7, 55, 55.36, ra=True)
delta_1 = Angle(21, 41, 3.0)
alpha_2 = Angle(7, 58, 22.55, ra=True)
delta_2 = Angle(21, 35, 23.4)
alpha_3 = Angle(8, 0, 48.99, ra=True)
delta_3 = Angle(21, 29, 38.2)
alpha_4 = Angle(8, 3, 14.66, ra=True)
delta_4 = Angle(21, 23, 47.5)
alpha_5 = Angle(8, 5, 39.54, ra=True)
delta_5 = Angle(21, 17, 51.4)
alpha_star1 = Angle(7, 34, 16.40, ra=True)
delta_star1 = Angle(31, 53, 51.2)
alpha_star2 = Angle(7, 45, 0.10, ra=True)
delta_star2 = Angle(28, 2, 12.5)
alpha_list = [alpha_1, alpha_2, alpha_3, alpha_4, alpha_5]
delta_list = [delta_1, delta_2, delta_3, delta_4, delta_5]
n = planet_stars_in_line(alpha_list, delta_list, alpha_star1, delta_star1,
alpha_star2, delta_star2)
print_me("Epoch fraction 'n' when bodies are in a straight line", round(n, 4))
# Epoch fraction 'n' when bodies are in a straight line: 0.2233
The function ‘straight_line()’ computes if three celestial bodies are in line providing the angle with which the bodies differ from a great circle:
alpha1 = Angle(5, 32, 0.40, ra=True)
delta1 = Angle(0, -17, 56.9)
alpha2 = Angle(5, 36, 12.81, ra=True)
delta2 = Angle(-1, 12, 7.0)
alpha3 = Angle(5, 40, 45.52, ra=True)
delta3 = Angle(-1, 56, 33.3)
psi, omega = straight_line(alpha1, delta1, alpha2, delta2, alpha3, delta3)
print_me("Angle deviation from a straight line", psi.dms_str(n_dec=0))
# Angle deviation from a straight line: 7d 31' 1.0''
print_me("Angular distance of central point to the straight line", omega.dms_str(n_dec=0))
# Angular distance of central point to the straight line: -5' 24.0''
Now let’s compute the size of the smallest circle that contains three given celestial bodies:
alpha1 = Angle(12, 41, 8.63, ra=True)
delta1 = Angle(-5, 37, 54.2)
alpha2 = Angle(12, 52, 5.21, ra=True)
delta2 = Angle(-4, 22, 26.2)
alpha3 = Angle(12, 39, 28.11, ra=True)
delta3 = Angle(-1, 50, 3.7)
d = circle_diameter(alpha1, delta1, alpha2, delta2, alpha3, delta3)
print_me("Diameter of smallest circle containing three celestial bodies", d.dms_str(n_dec=0))
# Diameter of smallest circle containing three celestial bodies: 4d 15' 49.0''
Let’s find the apparent position of a star (Theta Persei) for a given epoch:
epoch = Epoch(2028, 11, 13.19)
alpha = Angle(2, 46, 11.331, ra=True)
delta = Angle(49, 20, 54.54)
sun_lon = Angle(231.328)
app_alpha, app_delta = apparent_position(epoch, alpha, delta, sun_lon)
print_me("Apparent right ascension", app_alpha.ra_str(n_dec=2))
# Apparent right ascension: 2h 46' 14.39''
print_me("Apparent declination", app_delta.dms_str(n_dec=2))
# Apparent declination: 49d 21' 7.45''
Convert orbital elements of a celestial object from one equinox to another:
epoch0 = Epoch(2358042.5305)
epoch = Epoch(2433282.4235)
i0 = Angle(47.122)
arg0 = Angle(151.4486)
lon0 = Angle(45.7481)
i1, arg1, lon1 = orbital_equinox2equinox(epoch0, epoch, i0, arg0, lon0)
print_me("New inclination", round(i1(), 3))
# New inclination: 47.138
print_me("New argument of perihelion", round(arg1(), 4))
# New argument of perihelion: 151.4782
print_me("New longitude of ascending node", round(lon1(), 4))
# New longitude of ascending node: 48.6037
Compute the eccentric and true anomalies using Kepler’s equation:
eccentricity = 0.1
mean_anomaly = Angle(5.0)
e, v = kepler_equation(eccentricity, mean_anomaly)
print_me("Eccentric anomaly, Case #1", round(e(), 6))
# Eccentric anomaly, Case #1: 5.554589
print_me("True anomaly, Case #1", round(v(), 6))
# True anomaly, Case #1: 6.139762
e, v = kepler_equation(0.99, Angle(0.2, radians=True))
print_me("Eccentric anomaly, Case #2", round(e(), 8))
# Eccentric anomaly, Case #2: 61.13444578
print_me("True anomaly, Case #2", round(v(), 6))
# True anomaly, Case #2: 166.311977
Compute the velocity of a body in a given point of its (unperturbated elliptic) orbit, in this case the comet Halley in 1986:
r = 1.0
a = 17.9400782
v = velocity(r, a)
print_me("Velocity ar 1 AU", round(v, 2))
# Velocity at 1 AU: 41.53
Compute the velocity at perihelion:
a = 17.9400782
e = 0.96727426
vp = velocity_perihelion(e, a)
print_me("Velocity at perihelion", round(vp, 2))
# Velocity at perihelion: 54.52
And now compute the velocity at aphelion:
a = 17.9400782
e = 0.96727426
va = velocity_aphelion(e, a)
print_me("Velocity at aphelion", round(va, 2))
# Velocity at aphelion: 0.91
Calculate the length of the orbit for the same comet Halley:
a = 17.9400782
e = 0.96727426
length = length_orbit(e, a)
print_me("Length of the orbit (AU)", round(length, 2))
# Length of the orbit (AU): 77.06
Compute passage through the nodes of an elliptic orbit:
omega = Angle(111.84644)
e = 0.96727426
a = 17.9400782
t = Epoch(1986, 2, 9.45891)
time, r = passage_nodes_elliptic(omega, e, a, t)
y, m, d = time.get_date()
d = round(d, 2)
print("Time of passage through ascending node: {}/{}/{}".format(y, m, d))
# Time of passage through ascending node: 1985/11/9.16
print("Radius vector at ascending node: {}".format(round(r, 4)))
# Radius vector at ascending node: 1.8045
Passage through the nodes of a parabolic orbit:
omega = Angle(154.9103)
q = 1.324502
t = Epoch(1989, 8, 20.291)
time, r = passage_nodes_parabolic(omega, q, t, ascending=False)
y, m, d = time.get_date()
d = round(d, 2)
print("Time of passage through descending node: {}/{}/{}".format(y, m, d))
# Time of passage through descending node: 1989/9/17.64
print("Radius vector at descending node: {}".format(round(r, 4)))
# Radius vector at descending node: 1.3901
Compute the phase angle:
sun_dist = 0.724604
earth_dist = 0.910947
sun_earth_dist = 0.983824
angle = phase_angle(sun_dist, earth_dist, sun_earth_dist)
print_me("Phase angle", round(angle, 2))
# Phase angle: 72.96
Now, let’s compute the illuminated fraction of the disk:
k = illuminated_fraction(sun_dist, earth_dist, sun_earth_dist)
print_me("Illuminated fraction of planet disk", round(k, 3))
# Illuminated fraction of planet disk: 0.647