Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Vectorize mpc.mpc_orbit, mpc.comet_orbit #532

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Vectorize mpc.mpc_orbit, mpc.comet_orbit #532

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1c10961
Vectorize `mpc.mpc_orbit`, `mpc.comet_orbit`
JoshPaterson Jan 6, 2021
e717fc5
Remove unnecessary if-else
JoshPaterson Feb 24, 2021
2b8dbfa
Condense simple if-else statement to one line
JoshPaterson Feb 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
166 changes: 92 additions & 74 deletions skyfield/data/mpc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import re

import pandas as pd
import numpy as np

from ..data.spice import inertial_frames
from ..keplerlib import _KeplerOrbit
Expand Down Expand Up @@ -75,34 +76,59 @@ def load_mpcorb_dataframe(fobj):
return df

def mpcorb_orbit(row, ts, gm_km3_s2):
a = row.semimajor_axis_au
e = row.eccentricity
p = a * (1.0 - e*e)

def n(c):
return ord(c) - (48 if c.isdigit() else 55)

def d(s):
year = 100 * n(s[0]) + int(s[1:3])
month = n(s[3])
day = n(s[4])
return julian_day(year, month, day) - 0.5

epoch_jd = d(row.epoch_packed)
t_epoch = ts.tt_jd(epoch_jd)

minor_planet = _KeplerOrbit._from_mean_anomaly(
p,
e,
row.inclination_degrees,
row.longitude_of_ascending_node_degrees,
row.argument_of_perihelion_degrees,
row.mean_anomaly_degrees,
t_epoch,
gm_km3_s2,
10,
row.designation,
)
if isinstance(row, pd.Series):
a = row.semimajor_axis_au
e = row.eccentricity
p = a * (1.0 - e*e)

def d(s):
year = 100 * int(s[0], 32) + int(s[1:3])
month = int(s[3], 32)
day = int(s[4], 32)
return julian_day(year, month, day) - 0.5

epoch_jd = d(row.epoch_packed)
t_epoch = ts.tt_jd(epoch_jd)

minor_planet = _KeplerOrbit._from_mean_anomaly(
p,
e,
row.inclination_degrees,
row.longitude_of_ascending_node_degrees,
row.argument_of_perihelion_degrees,
row.mean_anomaly_degrees,
t_epoch,
gm_km3_s2,
10,
row.designation,
)
else:
a = row.semimajor_axis_au.values
e = row.eccentricity.values
p = a * (1.0 - e*e)

def d(s):
year = 100 * int(s[0], 32) + int(s[1:3])
month = int(s[3], 32)
day = int(s[4], 32)
return julian_day(year, month, day) - 0.5

epoch_jd = np.array([d(epoch) for epoch in row.epoch_packed.values])
t_epoch = ts.tt_jd(epoch_jd)

minor_planet = _KeplerOrbit._from_mean_anomaly(
p,
e,
row.inclination_degrees.values,
row.longitude_of_ascending_node_degrees.values,
row.argument_of_perihelion_degrees.values,
row.mean_anomaly_degrees.values,
t_epoch,
gm_km3_s2,
10,
row.designation.values,
)

minor_planet._rotation = inertial_frames['ECLIPJ2000'].T
return minor_planet

Expand Down Expand Up @@ -203,59 +229,51 @@ def load_comets_dataframe_slow(fobj):
return df

def comet_orbit(row, ts, gm_km3_s2):
e = row.eccentricity
if e == 1.0:
p = row.perihelion_distance_au * 2.0
if isinstance(row, pd.Series):
e = row.eccentricity
if e == 1:
p = row.perihelion_distance_au * 2.0
else:
p = row.perihelion_distance_au / (1.0 - e) * (1.0 - e*e)
t_perihelion = ts.tt(row.perihelion_year, row.perihelion_month,
row.perihelion_day)
comet = _KeplerOrbit._from_periapsis(
p,
e,
row.inclination_degrees,
row.longitude_of_ascending_node_degrees,
row.argument_of_perihelion_degrees,
t_perihelion,
gm_km3_s2,
10,
row['designation'],
)
else:
a = row.perihelion_distance_au / (1.0 - e)
p = a * (1.0 - e*e)
t_perihelion = ts.tt(row.perihelion_year, row.perihelion_month,
row.perihelion_day)

comet = _KeplerOrbit._from_periapsis(
p,
e,
row.inclination_degrees,
row.longitude_of_ascending_node_degrees,
row.argument_of_perihelion_degrees,
t_perihelion,
gm_km3_s2,
10,
row['designation'],
)
comet._rotation = inertial_frames['ECLIPJ2000'].T
return comet

def _comet_orbits(rows, ts, gm_km3_s2):
e = rows.eccentricity.values
parabolic = (e == 1.0)
p = (1 - e*e) / (1.0 - e + parabolic)
p[parabolic] += 2.0
p *= rows.perihelion_distance_au.values

t_perihelion = ts.tt(rows.perihelion_year.values, rows.perihelion_month.values,
rows.perihelion_day.values)

comet = _KeplerOrbit._from_periapsis(
p,
e,
rows.inclination_degrees.values,
rows.longitude_of_ascending_node_degrees.values,
rows.argument_of_perihelion_degrees.values,
t_perihelion,
gm_km3_s2,
10,
rows['designation'],
)
e = row.eccentricity.values
parabolic = (e == 1.0)
p = (1 - e*e) / (1.0 - e + parabolic)
p[parabolic] += 2.0
p *= row.perihelion_distance_au.values
t_perihelion = ts.tt(row.perihelion_year.values, row.perihelion_month.values,
row.perihelion_day.values)
comet = _KeplerOrbit._from_periapsis(
p,
e,
row.inclination_degrees.values,
row.longitude_of_ascending_node_degrees.values,
row.argument_of_perihelion_degrees.values,
t_perihelion,
gm_km3_s2,
10,
row['designation'].values,
)
comet._rotation = inertial_frames['ECLIPJ2000'].T
return comet

def unpack(designation_packed):
def n(c):
return ord(c) - (48 if c.isdigit() else 55)
s = designation_packed
s1 = s[1]
if s1 == '/':
return s
return '{0[0]}/{1}{0[2]}{0[3]} {0[4]}{2}{3}'.format(
s, n(s1), int(s[5:7]), s[7].lstrip('0'))
s, int(s1, 32), int(s[5:7]), s[7].lstrip('0'))
83 changes: 45 additions & 38 deletions skyfield/keplerlib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
import sys
import math
from numpy import(abs, amax, amin, arange, arccos, arctan, array, atleast_1d,
clip, copy, copyto, cos, cosh, exp, full_like, log, ndarray,
newaxis, pi, power, repeat, sin, sinh, squeeze, sqrt, sum,
tan, tanh, zeros_like)
clip, copy, copyto, cos, cosh, exp, float64, full_like, log,
ndarray, newaxis, pi, power, repeat, sin, sinh, squeeze,
sqrt, sum, tan, tanh, zeros_like)

from skyfield.constants import AU_KM, DAY_S, DEG2RAD
from skyfield.functions import dots, length_of, mxv
Expand Down Expand Up @@ -187,16 +187,13 @@ def _from_mean_anomaly(
target : int
NAIF ID of the secondary body
"""
M = DEG2RAD * mean_anomaly_degrees
gm_au3_d2 = gm_km3_s2 * _CONVERT_GM
if eccentricity < 1.0:
E = eccentric_anomaly(eccentricity, M)
v = true_anomaly_closed(eccentricity, E)
elif eccentricity > 1.0:
E = eccentric_anomaly(eccentricity, M)
v = true_anomaly_hyperbolic(eccentricity, E)
else:
v = true_anomaly_parabolic(semilatus_rectum_au, gm_au3_d2, M)
v = true_anomaly(
eccentricity,
DEG2RAD * mean_anomaly_degrees,
semilatus_rectum_au,
gm_au3_d2,
)

pos, vel = ele_to_vec(
semilatus_rectum_au,
Expand Down Expand Up @@ -268,6 +265,27 @@ def __repr__(self):
return '<{0}>'.format(str(self))


def true_anomaly(e, M, p, gm):
return_scalar = isinstance(e, (float, float64))
e, M, p = atleast_1d(e, M, p)

closed = e < 1.0
hyperbolic = e > 1.0
parabolic = ~closed & ~hyperbolic

v = zeros_like(e)

E_closed = eccentric_anomaly(e[closed], M[closed])
v[closed] = 2.0 * arctan(sqrt((1.0 + e[closed]) / (1.0 - e[closed])) * tan(E_closed/2))

E_hyperbolic = eccentric_anomaly(e[hyperbolic], M[hyperbolic])
v[hyperbolic] = 2.0 * arctan(sqrt((e[hyperbolic] + 1.0) / (e[hyperbolic] - 1.0)) * tanh(E_hyperbolic/2))

v[parabolic] = true_anomaly_parabolic(p[parabolic], gm, M[parabolic])

return v[0] if return_scalar else v


def eccentric_anomaly(e, M):
""" Iterates to solve Kepler's equation to find eccentric anomaly

Expand All @@ -278,33 +296,22 @@ def eccentric_anomaly(e, M):
E = M + e*sin(M)

max_iters = 100
iters = 0
dM = M - (E - e*sin(E))
dE = dM/(1 - e*cos(E))
not_done = abs(dE) > 1e-14
iters = 1
while iters < max_iters:
dM = M - (E - e*sin(E))
dE = dM/(1 - e*cos(E))
E = E + dE
if abs(dE) < 1e-14: return E
if not not_done.any():
break
dM[not_done] = M[not_done] - (E[not_done] - e[not_done]*sin(E[not_done]))
dE[not_done] = dM[not_done]/(1 - e[not_done]*cos(E[not_done]))
E[not_done] += dE[not_done]
iters += 1
else:
raise ValueError('Failed to converge')


def true_anomaly_hyperbolic(e, E):
"""Calculates true anomaly from eccentricity and eccentric anomaly.

Valid for hyperbolic orbits. Equations from the relevant Wikipedia entries.

"""
return 2.0 * arctan(sqrt((e + 1.0) / (e - 1.0)) * tanh(E/2))


def true_anomaly_closed(e, E):
"""Calculates true anomaly from eccentricity and eccentric anomaly.
not_done = abs(dE) > 1e-14

Valid for closed orbits. Equations from the relevant Wikipedia entries.

"""
return 2.0 * arctan(sqrt((1.0 + e) / (1.0 - e)) * tan(E/2))
else:
raise ValueError('failed to converge')
return E


def true_anomaly_parabolic(p, gm, M):
Expand Down Expand Up @@ -611,7 +618,7 @@ def kepler_1d(x, orb_inds):
pcdot = -qovr0 / br * x * c1
vcdot = 1 - bq / br * x * x * c2

position_prop = pc[newaxis, :, :]*position[:, :, newaxis] + vc[newaxis, :, :]*velocity[:, :, newaxis]
velocity_prop = pcdot[newaxis, :, :]*position[:, :, newaxis] + vcdot[newaxis, :, :]*velocity[:, :, newaxis]
position_prop = pc.T[newaxis, :, :]*position[:, newaxis, :] + vc.T[newaxis, :, :]*velocity[:, newaxis, :]
velocity_prop = pcdot.T[newaxis, :, :]*position[:, newaxis, :] + vcdot.T[newaxis, :, :]*velocity[:, newaxis, :]

return squeeze(position_prop), squeeze(velocity_prop)
56 changes: 56 additions & 0 deletions skyfield/tests/test_keplerlib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,31 @@ def test_minor_planet():
assert abs(ra.hours - 23.1437) < 0.00005
assert abs(dec.degrees - -17.323) < 0.0005

def test_minor_planets():
text = (b'00001 3.4 0.15 K205V 162.68631 73.73161 80.28698'
b' 10.58862 0.0775571 0.21406009 2.7676569 0 MPO492748'
b' 6751 115 1801-2019 0.60 M-v 30h Williams 0000 '
b'(1) Ceres 20190915\n'
b'00001 3.4 0.15 K205V 162.68631 73.73161 80.28698'
b' 10.58862 0.0775571 0.21406009 2.7676569 0 MPO492748'
b' 6751 115 1801-2019 0.60 M-v 30h Williams 0000 '
b'(1) Ceres 20190915\n')

ts = load.timescale()
t = ts.utc(2020, 6, 17)
eph = load('de421.bsp')
df = mpc.load_mpcorb_dataframe(BytesIO(text))

assert (df.designation_packed.values == '00001').all()
assert (df.designation.values == '(1) Ceres').all()

ceres = mpc.mpcorb_orbit(df, ts, GM_SUN)
ra, dec, distance = eph['earth'].at(t).observe(eph['sun'] + ceres).radec()

assert (ceres.target == '(1) Ceres').all()
assert (abs(ra.hours - 23.1437) < 0.00005).all()
assert (abs(dec.degrees - -17.323) < 0.0005).all()

def test_comet():
text = (b' CJ95O010 1997 03 29.6333 0.916241 0.994928 130.6448'
b' 283.3593 88.9908 20200224 -2.0 4.0 C/1995 O1 (Hale-Bopp)'
Expand Down Expand Up @@ -103,6 +128,37 @@ def test_comet():

assert k.target == 'C/1995 O1 (Hale-Bopp)'

def test_comets():
text = (b' CJ95O010 1997 03 29.6333 0.916241 0.994928 130.6448'
b' 283.3593 88.9908 20200224 -2.0 4.0 C/1995 O1 (Hale-Bopp)'
b' MPC106342\n'
b' CJ95O010 1997 03 29.6333 0.916241 0.994928 130.6448'
b' 283.3593 88.9908 20200224 -2.0 4.0 C/1995 O1 (Hale-Bopp)'
b' MPC106342\n')

ts = load.timescale()
t = ts.utc(2020, 5, 31)
eph = load('de421.bsp')
e = eph['earth'].at(t)

for loader in mpc.load_comets_dataframe, mpc.load_comets_dataframe_slow:
df = loader(BytesIO(text))
k = mpc.comet_orbit(df, ts, GM_SUN)
p = e.observe(eph['sun'] + k)
ra, dec, distance = p.radec()

# The file authorities/mpc-hale-bopp in the repository is the
# source of these angles. TODO: can we tighten this bound and
# drive it to fractions of an arcsecond?

ra_want = Angle(hours=(23, 59, 16.6))
dec_want = Angle(degrees=(-84, 46, 58))
assert (abs(ra_want.arcseconds() - ra.arcseconds()) < 2.0).all()
assert (abs(dec_want.arcseconds() - dec.arcseconds()) < 0.2).all()
assert (abs(distance.au - 43.266) < 0.0005).all()

assert (k.target == 'C/1995 O1 (Hale-Bopp)').all()

def test_comet_with_eccentricity_of_exactly_one():
ts = load.timescale()
t = ts.utc(2020, 8, 13)
Expand Down