diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index 40d9ba73d..68b49164e 100644 --- a/skyfield/data/mpc.py +++ b/skyfield/data/mpc.py @@ -5,6 +5,7 @@ import re import pandas as pd +import numpy as np from ..data.spice import inertial_frames from ..keplerlib import _KeplerOrbit @@ -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 @@ -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')) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 08c9957de..9ec83b68c 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -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 @@ -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, @@ -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 @@ -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): @@ -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) diff --git a/skyfield/tests/test_keplerlib.py b/skyfield/tests/test_keplerlib.py index 00d31ceff..c3411194b 100644 --- a/skyfield/tests/test_keplerlib.py +++ b/skyfield/tests/test_keplerlib.py @@ -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)' @@ -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)