From a979f929bb9dfefd88a28dae18ee157f23ea574c Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Fri, 1 Jan 2021 21:33:12 -0500 Subject: [PATCH 01/20] change output shape of `propagate` to (xyz, n, t) --- skyfield/keplerlib.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 08c9957de..0f32f9766 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -611,7 +611,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) From 6cf0faecc02a15ff0c35388bc89e5abccbd8222f Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Fri, 1 Jan 2021 21:34:42 -0500 Subject: [PATCH 02/20] Make `VectorSum._at` work for 3 dimensional arrays --- skyfield/vectorlib.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index 57ee37a0b..cbd8e7847 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -5,7 +5,7 @@ from .constants import C_AUDAY from .descriptorlib import reify from .errors import DeprecationError -from .functions import length_of +from .functions import length_of, _reconcile from .positionlib import build_position from .timelib import Time @@ -215,8 +215,10 @@ def _at(self, t): p2, v2, another_gcrs_position, message = vf._at(t) if gcrs_position is None: # TODO: so bootleg; rework whole idea gcrs_position = another_gcrs_position - p += p2 - v += v2 + p, p2 = _reconcile(p, p2) + p = p + p2 + v, v2 = _reconcile(v, v2) + v = v + v2 return p, v, gcrs_position, message def _correct_for_light_travel_time(observer, target): From 3a5ba1420f31b6f76ac405f6abb21da2e1087507 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Fri, 1 Jan 2021 22:05:55 -0500 Subject: [PATCH 03/20] vectorize `mpc.comet_orbit` --- skyfield/data/mpc.py | 80 ++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 43 deletions(-) diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index 40d9ba73d..18a766db0 100644 --- a/skyfield/data/mpc.py +++ b/skyfield/data/mpc.py @@ -203,50 +203,44 @@ 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'], + ) comet._rotation = inertial_frames['ECLIPJ2000'].T return comet From 6f15401d1ee413eeefeeaf6e780e5af9b0513d2f Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:31:52 -0500 Subject: [PATCH 04/20] Refactor function `n()` --- skyfield/data/mpc.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index 18a766db0..e35df8a57 100644 --- a/skyfield/data/mpc.py +++ b/skyfield/data/mpc.py @@ -245,11 +245,9 @@ def comet_orbit(row, ts, gm_km3_s2): 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')) From 4e775b63c4041cf3e6dde67a92f6720959cbb0fc Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:32:37 -0500 Subject: [PATCH 05/20] Vectorize `mpc.mpc_orbit` --- skyfield/data/mpc.py | 82 +++++++++++++++++++++++++++++--------------- 1 file changed, 54 insertions(+), 28 deletions(-) diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index e35df8a57..ae21d0775 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 From 152fe9c5ecba16304d8c0d8f6eed2796a1247b2f Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:34:31 -0500 Subject: [PATCH 06/20] vectorize `_KeplerOrbit.from_mean_anomaly` --- skyfield/keplerlib.py | 56 ++++++++++++++++++++++++++++++------------- 1 file changed, 39 insertions(+), 17 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 0f32f9766..affc9ced9 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,31 @@ def _from_mean_anomaly( target : int NAIF ID of the secondary body """ + return_scalar = True if isinstance(eccentricity, (float, float64)) else False + eccentricity = atleast_1d(eccentricity) + mean_anomaly_degrees = atleast_1d(mean_anomaly_degrees) + semilatus_rectum_au = atleast_1d(semilatus_rectum_au) + 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) + + closed = eccentricity < 1.0 + hyperbolic = eccentricity > 1.0 + parabolic = ~closed & ~hyperbolic + + v = zeros_like(eccentricity) + + E = eccentric_anomaly(eccentricity[closed], M[closed]) + v[closed] = true_anomaly_closed(eccentricity[closed], E) + + E = eccentric_anomaly(eccentricity[hyperbolic], M[hyperbolic]) + v[hyperbolic] = true_anomaly_hyperbolic(eccentricity[hyperbolic], E) + + v[parabolic] = true_anomaly_parabolic(semilatus_rectum_au[parabolic], gm_au3_d2, M[parabolic]) + + if return_scalar: + v = v[0] + semilatus_rectum_au = semilatus_rectum_au[0] pos, vel = ele_to_vec( semilatus_rectum_au, @@ -278,15 +293,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 + not_done = abs(dE) > 1e-14 + else: - raise ValueError('Failed to converge') + raise ValueError('failed to converge') + return E def true_anomaly_hyperbolic(e, E): From 960447a64554ae700c4c3a878ec32eb7a4f3001c Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:43:17 -0500 Subject: [PATCH 07/20] Pull all true anomaly calcs out into a function --- skyfield/keplerlib.py | 54 ++++++++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 24 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index affc9ced9..60a5aef68 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -187,31 +187,13 @@ def _from_mean_anomaly( target : int NAIF ID of the secondary body """ - return_scalar = True if isinstance(eccentricity, (float, float64)) else False - eccentricity = atleast_1d(eccentricity) - mean_anomaly_degrees = atleast_1d(mean_anomaly_degrees) - semilatus_rectum_au = atleast_1d(semilatus_rectum_au) - - M = DEG2RAD * mean_anomaly_degrees gm_au3_d2 = gm_km3_s2 * _CONVERT_GM - - closed = eccentricity < 1.0 - hyperbolic = eccentricity > 1.0 - parabolic = ~closed & ~hyperbolic - - v = zeros_like(eccentricity) - - E = eccentric_anomaly(eccentricity[closed], M[closed]) - v[closed] = true_anomaly_closed(eccentricity[closed], E) - - E = eccentric_anomaly(eccentricity[hyperbolic], M[hyperbolic]) - v[hyperbolic] = true_anomaly_hyperbolic(eccentricity[hyperbolic], E) - - v[parabolic] = true_anomaly_parabolic(semilatus_rectum_au[parabolic], gm_au3_d2, M[parabolic]) - - if return_scalar: - v = v[0] - semilatus_rectum_au = semilatus_rectum_au[0] + v = true_anomaly( + eccentricity, + DEG2RAD * mean_anomaly_degrees, + semilatus_rectum_au, + gm_au3_d2, + ) pos, vel = ele_to_vec( semilatus_rectum_au, @@ -283,6 +265,30 @@ def __repr__(self): return '<{0}>'.format(str(self)) +def true_anomaly(e, M, p, gm): + return_scalar = True if isinstance(e, (float, float64)) else False + e, M, p = atleast_1d(e, M, p) + + closed = e < 1.0 + hyperbolic = e > 1.0 + parabolic = ~closed & ~hyperbolic + + v = zeros_like(e) + + E = eccentric_anomaly(e[closed], M[closed]) + v[closed] = true_anomaly_closed(e[closed], E) + + E = eccentric_anomaly(e[hyperbolic], M[hyperbolic]) + v[hyperbolic] = true_anomaly_hyperbolic(e[hyperbolic], E) + + v[parabolic] = true_anomaly_parabolic(p[parabolic], gm, M[parabolic]) + + if return_scalar: + return v[0] + else: + return v + + def eccentric_anomaly(e, M): """ Iterates to solve Kepler's equation to find eccentric anomaly From 8e7fae6a79d6479d43f766536c114be7a6ef84fe Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:46:30 -0500 Subject: [PATCH 08/20] Remove `true_anomaly_hyperbolic` --- skyfield/keplerlib.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 60a5aef68..6b6230a62 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -279,7 +279,7 @@ def true_anomaly(e, M, p, gm): v[closed] = true_anomaly_closed(e[closed], E) E = eccentric_anomaly(e[hyperbolic], M[hyperbolic]) - v[hyperbolic] = true_anomaly_hyperbolic(e[hyperbolic], E) + v[hyperbolic] = 2.0 * arctan(sqrt((e[hyperbolic] + 1.0) / (e[hyperbolic] - 1.0)) * tanh(E/2)) v[parabolic] = true_anomaly_parabolic(p[parabolic], gm, M[parabolic]) @@ -317,15 +317,6 @@ def eccentric_anomaly(e, M): return E -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. From 558d682d624fe3d9c118a728267dad916ebdb948 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:47:46 -0500 Subject: [PATCH 09/20] Remove `true_anomaly_closed` --- skyfield/keplerlib.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 6b6230a62..f656a3b55 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -276,7 +276,7 @@ def true_anomaly(e, M, p, gm): v = zeros_like(e) E = eccentric_anomaly(e[closed], M[closed]) - v[closed] = true_anomaly_closed(e[closed], E) + v[closed] = 2.0 * arctan(sqrt((1.0 + e[closed]) / (1.0 - e[closed])) * tan(E/2)) E = eccentric_anomaly(e[hyperbolic], M[hyperbolic]) v[hyperbolic] = 2.0 * arctan(sqrt((e[hyperbolic] + 1.0) / (e[hyperbolic] - 1.0)) * tanh(E/2)) @@ -317,15 +317,6 @@ def eccentric_anomaly(e, M): return E -def true_anomaly_closed(e, E): - """Calculates true anomaly from eccentricity and eccentric anomaly. - - Valid for closed orbits. Equations from the relevant Wikipedia entries. - - """ - return 2.0 * arctan(sqrt((1.0 + e) / (1.0 - e)) * tan(E/2)) - - def true_anomaly_parabolic(p, gm, M): """Calculates true anomaly from semi-latus rectum, gm, and mean anomaly. From 8343debd58e6b7d3724be41531252df11696a318 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 00:49:57 -0500 Subject: [PATCH 10/20] Change variable names --- skyfield/keplerlib.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index f656a3b55..428591817 100644 --- a/skyfield/keplerlib.py +++ b/skyfield/keplerlib.py @@ -275,11 +275,11 @@ def true_anomaly(e, M, p, gm): v = zeros_like(e) - E = eccentric_anomaly(e[closed], M[closed]) - v[closed] = 2.0 * arctan(sqrt((1.0 + e[closed]) / (1.0 - e[closed])) * tan(E/2)) + 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 = eccentric_anomaly(e[hyperbolic], M[hyperbolic]) - v[hyperbolic] = 2.0 * arctan(sqrt((e[hyperbolic] + 1.0) / (e[hyperbolic] - 1.0)) * tanh(E/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]) From c46a849c7e3ff5836152adbfe1a098470504fe44 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 18:55:17 -0500 Subject: [PATCH 11/20] Add tests for multiple comets and minor planets --- skyfield/tests/test_keplerlib.py | 56 ++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) 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) From 59b4937c7a763ae9ccf8273e12a7db477079b5fb Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sat, 2 Jan 2021 19:18:51 -0500 Subject: [PATCH 12/20] Add use of `_reconcile` --- skyfield/vectorlib.py | 1 + 1 file changed, 1 insertion(+) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index cbd8e7847..fa24049f8 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -241,6 +241,7 @@ def _correct_for_light_travel_time(observer, target): tposition, tvelocity, gcrs_position, message = target._at(t) + tposition, cposition = _reconcile(tposition, cposition) distance = length_of(tposition - cposition) light_time0 = 0.0 for i in range(10): From ce963d2880b8312921dc69530212e85a5e22c140 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sun, 3 Jan 2021 03:47:28 -0500 Subject: [PATCH 13/20] Start to vectorize _correct_for_light_travel_time --- skyfield/vectorlib.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index fa24049f8..70ebfbf80 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -1,7 +1,7 @@ """Vector functions and their composition.""" from jplephem.names import target_names as _jpl_code_name_dict -from numpy import max +from numpy import max, newaxis, expand_dims from .constants import C_AUDAY from .descriptorlib import reify from .errors import DeprecationError @@ -236,12 +236,19 @@ def _correct_for_light_travel_time(observer, target): whole = t.whole tdb_fraction = t.tdb_fraction + whole = whole[:, newaxis, newaxis] + tdb_fraction = tdb_fraction[:, newaxis, newaxis] + cposition = observer.position.au cvelocity = observer.velocity.au_per_d + cposition = expand_dims(cposition, 2) + cvelocity = expand_dims(cvelocity, 2) + tposition, tvelocity, gcrs_position, message = target._at(t) + tposition = expand_dims(tposition, 3) + tvelocity = expand_dims(tvelocity, 3) - tposition, cposition = _reconcile(tposition, cposition) distance = length_of(tposition - cposition) light_time0 = 0.0 for i in range(10): From eaae2131d46e2de36b6f466b9c661c4750771421 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Mon, 4 Jan 2021 21:34:03 -0500 Subject: [PATCH 14/20] Finish vectorizing _correct_for_light_travel_time --- skyfield/vectorlib.py | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index 70ebfbf80..8451ec693 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -1,13 +1,14 @@ """Vector functions and their composition.""" from jplephem.names import target_names as _jpl_code_name_dict -from numpy import max, newaxis, expand_dims +from numpy import max, newaxis, expand_dims, zeros, zeros_like, broadcast_to, reshape from .constants import C_AUDAY from .descriptorlib import reify from .errors import DeprecationError from .functions import length_of, _reconcile from .positionlib import build_position from .timelib import Time +from .units import Distance, Velocity class VectorFunction(object): """Given a time, computes a corresponding position.""" @@ -236,9 +237,6 @@ def _correct_for_light_travel_time(observer, target): whole = t.whole tdb_fraction = t.tdb_fraction - whole = whole[:, newaxis, newaxis] - tdb_fraction = tdb_fraction[:, newaxis, newaxis] - cposition = observer.position.au cvelocity = observer.velocity.au_per_d @@ -251,6 +249,13 @@ def _correct_for_light_travel_time(observer, target): distance = length_of(tposition - cposition) light_time0 = 0.0 + + whole = broadcast_to(whole[:, newaxis, newaxis], (len(t.tt), tposition.shape[2], cposition.shape[3])) + tdb_fraction = tdb_fraction[:, newaxis, newaxis] + + tposition_out = zeros((3, len(t.tt), tposition.shape[2], cposition.shape[3])) + tvelocity_out = zeros_like(tposition_out) + for i in range(10): light_time = distance / C_AUDAY delta = light_time - light_time0 @@ -262,12 +267,22 @@ def _correct_for_light_travel_time(observer, target): # fraction, for adding to the whole and fraction of TDB. t2 = ts.tdb_jd(whole, tdb_fraction - light_time) - tposition, tvelocity, gcrs_position, message = target._at(t2) - distance = length_of(tposition - cposition) + for i, (pos, vel, epoch) in enumerate(zip(target.vector_functions[-1].position_at_epoch.au.T, + target.vector_functions[-1].velocity_at_epoch.au_per_d.T, + target.vector_functions[-1].epoch.tt, + )): + k_orbit = type(target.vector_functions[-1]) #avoid import loop? + target_i = target.vector_functions[0] + k_orbit(Distance(au=pos), Velocity(au_per_d=vel), ts.tt_jd(epoch), target.vector_functions[-1].mu_au3_d2, center=10) + t2_i_flat = ts.tt_jd(t2.tt[:, i, :].flatten()) + tposition_i_flat, tvelocity_i_flat, _, _ = target_i._at(t2_i_flat) + tposition_out[:, :, i, :] = reshape(tposition_i_flat, (3,) + t2.tt[:, i, :].shape) + tvelocity_out[:, :, i, :] = reshape(tvelocity_i_flat, (3,) + t2.tt[:, i, :].shape) + + distance = length_of(tposition_out - cposition) light_time0 = light_time else: raise ValueError('light-travel time failed to converge') - return tposition - cposition, tvelocity - cvelocity, t, light_time + return tposition_out - cposition, tvelocity_out - cvelocity, t, light_time def _jpl_name(target): if not isinstance(target, int): From 740b4e728191bef4b9543ef70504186aa3997831 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Tue, 5 Jan 2021 00:22:31 -0500 Subject: [PATCH 15/20] Allow non `int` centers for `Astrometric` --- skyfield/positionlib.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skyfield/positionlib.py b/skyfield/positionlib.py index 31a7887cb..3615be825 100644 --- a/skyfield/positionlib.py +++ b/skyfield/positionlib.py @@ -101,7 +101,7 @@ def __init__(self, position_au, velocity_au_per_d=None, t=None, self.velocity = Velocity(velocity_au_per_d) self.center = self._default_center if center is None else center self.target = target - if center == 0: + if isinstance(center, int) and center == 0: self.center_barycentric = self @classmethod From fb5e7cdf13f55ce618ab1a67e780d0604d4f8719 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Tue, 5 Jan 2021 01:24:44 -0500 Subject: [PATCH 16/20] Fix _KeplerOrbit rather than change ICRF.__init__ --- skyfield/data/mpc.py | 2 -- skyfield/positionlib.py | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index ae21d0775..410cc0dd8 100644 --- a/skyfield/data/mpc.py +++ b/skyfield/data/mpc.py @@ -126,7 +126,6 @@ def d(s): t_epoch, gm_km3_s2, 10, - row.designation.values, ) minor_planet._rotation = inertial_frames['ECLIPJ2000'].T @@ -265,7 +264,6 @@ def comet_orbit(row, ts, gm_km3_s2): t_perihelion, gm_km3_s2, 10, - row['designation'], ) comet._rotation = inertial_frames['ECLIPJ2000'].T return comet diff --git a/skyfield/positionlib.py b/skyfield/positionlib.py index 3615be825..31a7887cb 100644 --- a/skyfield/positionlib.py +++ b/skyfield/positionlib.py @@ -101,7 +101,7 @@ def __init__(self, position_au, velocity_au_per_d=None, t=None, self.velocity = Velocity(velocity_au_per_d) self.center = self._default_center if center is None else center self.target = target - if isinstance(center, int) and center == 0: + if center == 0: self.center_barycentric = self @classmethod From 4c13530ee2e3c938fd518aa9f3d6017d5574a9f4 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Tue, 5 Jan 2021 01:40:34 -0500 Subject: [PATCH 17/20] Add descriptive variable names --- skyfield/vectorlib.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index 8451ec693..df3f11b9a 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -237,23 +237,30 @@ def _correct_for_light_travel_time(observer, target): whole = t.whole tdb_fraction = t.tdb_fraction + num_times = len(t.tt) + cposition = observer.position.au cvelocity = observer.velocity.au_per_d + num_observers = cposition.shape[2] + cposition = expand_dims(cposition, 2) cvelocity = expand_dims(cvelocity, 2) tposition, tvelocity, gcrs_position, message = target._at(t) + + num_targets = tposition.shape[2] + tposition = expand_dims(tposition, 3) tvelocity = expand_dims(tvelocity, 3) distance = length_of(tposition - cposition) light_time0 = 0.0 - whole = broadcast_to(whole[:, newaxis, newaxis], (len(t.tt), tposition.shape[2], cposition.shape[3])) + whole = broadcast_to(whole[:, newaxis, newaxis], (num_times, num_targets, num_observers)) tdb_fraction = tdb_fraction[:, newaxis, newaxis] - tposition_out = zeros((3, len(t.tt), tposition.shape[2], cposition.shape[3])) + tposition_out = zeros((3, num_times, num_targets, num_observers)) tvelocity_out = zeros_like(tposition_out) for i in range(10): @@ -275,8 +282,8 @@ def _correct_for_light_travel_time(observer, target): target_i = target.vector_functions[0] + k_orbit(Distance(au=pos), Velocity(au_per_d=vel), ts.tt_jd(epoch), target.vector_functions[-1].mu_au3_d2, center=10) t2_i_flat = ts.tt_jd(t2.tt[:, i, :].flatten()) tposition_i_flat, tvelocity_i_flat, _, _ = target_i._at(t2_i_flat) - tposition_out[:, :, i, :] = reshape(tposition_i_flat, (3,) + t2.tt[:, i, :].shape) - tvelocity_out[:, :, i, :] = reshape(tvelocity_i_flat, (3,) + t2.tt[:, i, :].shape) + tposition_out[:, :, i, :] = reshape(tposition_i_flat, (3, num_times, num_observers)) + tvelocity_out[:, :, i, :] = reshape(tvelocity_i_flat, (3, num_times, num_observers)) distance = length_of(tposition_out - cposition) light_time0 = light_time From 8aefce4b5f3e518c4085d233c80e40ff2cd40a85 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Wed, 6 Jan 2021 02:02:03 -0500 Subject: [PATCH 18/20] Revert changes outside vectorlib --- skyfield/data/mpc.py | 164 ++++++++++++++----------------- skyfield/keplerlib.py | 86 +++++++--------- skyfield/tests/test_keplerlib.py | 56 ----------- 3 files changed, 112 insertions(+), 194 deletions(-) diff --git a/skyfield/data/mpc.py b/skyfield/data/mpc.py index 410cc0dd8..40d9ba73d 100644 --- a/skyfield/data/mpc.py +++ b/skyfield/data/mpc.py @@ -5,7 +5,6 @@ import re import pandas as pd -import numpy as np from ..data.spice import inertial_frames from ..keplerlib import _KeplerOrbit @@ -76,58 +75,34 @@ def load_mpcorb_dataframe(fobj): return df def mpcorb_orbit(row, ts, gm_km3_s2): - 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, - ) - + 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, + ) minor_planet._rotation = inertial_frames['ECLIPJ2000'].T return minor_planet @@ -228,50 +203,59 @@ def load_comets_dataframe_slow(fobj): return df def comet_orbit(row, ts, gm_km3_s2): - 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'], - ) + e = row.eccentricity + if e == 1.0: + p = row.perihelion_distance_au * 2.0 else: - 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, - ) + 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'], + ) 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, int(s1, 32), int(s[5:7]), s[7].lstrip('0')) + s, n(s1), int(s[5:7]), s[7].lstrip('0')) diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py index 428591817..08c9957de 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, float64, full_like, log, - ndarray, newaxis, pi, power, repeat, sin, sinh, squeeze, - sqrt, sum, tan, tanh, zeros_like) + clip, copy, copyto, cos, cosh, exp, 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,13 +187,16 @@ 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 - v = true_anomaly( - eccentricity, - DEG2RAD * mean_anomaly_degrees, - semilatus_rectum_au, - gm_au3_d2, - ) + 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) pos, vel = ele_to_vec( semilatus_rectum_au, @@ -265,30 +268,6 @@ def __repr__(self): return '<{0}>'.format(str(self)) -def true_anomaly(e, M, p, gm): - return_scalar = True if isinstance(e, (float, float64)) else False - 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]) - - if return_scalar: - return v[0] - else: - return v - - def eccentric_anomaly(e, M): """ Iterates to solve Kepler's equation to find eccentric anomaly @@ -299,22 +278,33 @@ def eccentric_anomaly(e, M): E = M + e*sin(M) max_iters = 100 - dM = M - (E - e*sin(E)) - dE = dM/(1 - e*cos(E)) - not_done = abs(dE) > 1e-14 - iters = 1 + iters = 0 while iters < max_iters: - 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] + dM = M - (E - e*sin(E)) + dE = dM/(1 - e*cos(E)) + E = E + dE + if abs(dE) < 1e-14: return E iters += 1 - not_done = abs(dE) > 1e-14 - else: - raise ValueError('failed to converge') - return E + 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. + + Valid for closed orbits. Equations from the relevant Wikipedia entries. + + """ + return 2.0 * arctan(sqrt((1.0 + e) / (1.0 - e)) * tan(E/2)) def true_anomaly_parabolic(p, gm, M): @@ -621,7 +611,7 @@ def kepler_1d(x, orb_inds): pcdot = -qovr0 / br * x * c1 vcdot = 1 - bq / br * x * x * c2 - position_prop = pc.T[newaxis, :, :]*position[:, newaxis, :] + vc.T[newaxis, :, :]*velocity[:, newaxis, :] - velocity_prop = pcdot.T[newaxis, :, :]*position[:, newaxis, :] + vcdot.T[newaxis, :, :]*velocity[:, newaxis, :] + position_prop = pc[newaxis, :, :]*position[:, :, newaxis] + vc[newaxis, :, :]*velocity[:, :, newaxis] + velocity_prop = pcdot[newaxis, :, :]*position[:, :, newaxis] + vcdot[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 c3411194b..00d31ceff 100644 --- a/skyfield/tests/test_keplerlib.py +++ b/skyfield/tests/test_keplerlib.py @@ -74,31 +74,6 @@ 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)' @@ -128,37 +103,6 @@ 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) From b5f4fd8ae4462e28d9fbb835c6a9fcdbbef1d730 Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sun, 21 Feb 2021 17:23:24 -0500 Subject: [PATCH 19/20] Remove unneccesary stuff --- skyfield/vectorlib.py | 29 +++++------------------------ 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index df3f11b9a..c3ba2832b 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -1,7 +1,7 @@ """Vector functions and their composition.""" from jplephem.names import target_names as _jpl_code_name_dict -from numpy import max, newaxis, expand_dims, zeros, zeros_like, broadcast_to, reshape +from numpy import max, newaxis, expand_dims, zeros, zeros_like, broadcast_to, diagonal from .constants import C_AUDAY from .descriptorlib import reify from .errors import DeprecationError @@ -237,32 +237,23 @@ def _correct_for_light_travel_time(observer, target): whole = t.whole tdb_fraction = t.tdb_fraction - num_times = len(t.tt) - cposition = observer.position.au cvelocity = observer.velocity.au_per_d - num_observers = cposition.shape[2] - cposition = expand_dims(cposition, 2) cvelocity = expand_dims(cvelocity, 2) tposition, tvelocity, gcrs_position, message = target._at(t) - num_targets = tposition.shape[2] - tposition = expand_dims(tposition, 3) tvelocity = expand_dims(tvelocity, 3) distance = length_of(tposition - cposition) light_time0 = 0.0 - whole = broadcast_to(whole[:, newaxis, newaxis], (num_times, num_targets, num_observers)) + whole = broadcast_to(whole[:, newaxis, newaxis], (len(t.tt), tposition.shape[2], cposition.shape[2])) tdb_fraction = tdb_fraction[:, newaxis, newaxis] - tposition_out = zeros((3, num_times, num_targets, num_observers)) - tvelocity_out = zeros_like(tposition_out) - for i in range(10): light_time = distance / C_AUDAY delta = light_time - light_time0 @@ -274,22 +265,12 @@ def _correct_for_light_travel_time(observer, target): # fraction, for adding to the whole and fraction of TDB. t2 = ts.tdb_jd(whole, tdb_fraction - light_time) - for i, (pos, vel, epoch) in enumerate(zip(target.vector_functions[-1].position_at_epoch.au.T, - target.vector_functions[-1].velocity_at_epoch.au_per_d.T, - target.vector_functions[-1].epoch.tt, - )): - k_orbit = type(target.vector_functions[-1]) #avoid import loop? - target_i = target.vector_functions[0] + k_orbit(Distance(au=pos), Velocity(au_per_d=vel), ts.tt_jd(epoch), target.vector_functions[-1].mu_au3_d2, center=10) - t2_i_flat = ts.tt_jd(t2.tt[:, i, :].flatten()) - tposition_i_flat, tvelocity_i_flat, _, _ = target_i._at(t2_i_flat) - tposition_out[:, :, i, :] = reshape(tposition_i_flat, (3, num_times, num_observers)) - tvelocity_out[:, :, i, :] = reshape(tvelocity_i_flat, (3, num_times, num_observers)) - - distance = length_of(tposition_out - cposition) + tposition, tvelocity, gcrs_position, message = diagonal(target._at(t2), axis1=1, axis2=2) + distance = length_of(tposition - cposition) light_time0 = light_time else: raise ValueError('light-travel time failed to converge') - return tposition_out - cposition, tvelocity_out - cvelocity, t, light_time + return tposition - cposition, tvelocity - cvelocity, t, light_time def _jpl_name(target): if not isinstance(target, int): From 98c50b0e46a51ae8e516659679122edde6e974fc Mon Sep 17 00:00:00 2001 From: Josh Paterson Date: Sun, 21 Feb 2021 17:25:29 -0500 Subject: [PATCH 20/20] Squeeze return values --- skyfield/vectorlib.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skyfield/vectorlib.py b/skyfield/vectorlib.py index c3ba2832b..3941a388a 100644 --- a/skyfield/vectorlib.py +++ b/skyfield/vectorlib.py @@ -1,7 +1,7 @@ """Vector functions and their composition.""" from jplephem.names import target_names as _jpl_code_name_dict -from numpy import max, newaxis, expand_dims, zeros, zeros_like, broadcast_to, diagonal +from numpy import max, newaxis, expand_dims, broadcast_to, diagonal, squeeze from .constants import C_AUDAY from .descriptorlib import reify from .errors import DeprecationError @@ -270,7 +270,7 @@ def _correct_for_light_travel_time(observer, target): light_time0 = light_time else: raise ValueError('light-travel time failed to converge') - return tposition - cposition, tvelocity - cvelocity, t, light_time + return squeeze(tposition - cposition), squeeze(tvelocity - cvelocity), t, light_time def _jpl_name(target): if not isinstance(target, int):