Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_propagating_to_certain_nu_is_correct():
# take an elliptic orbit
a = 1.0 * u.AU
ecc = 1.0 / 3.0 * u.one
_a = 0.0 * u.rad
nu = 10 * u.deg
elliptic = Orbit.from_classical(Sun, a, ecc, _a, _a, _a, nu)
elliptic_at_perihelion = elliptic.propagate_to_anomaly(0.0 * u.rad)
r_per, _ = elliptic_at_perihelion.rv()
elliptic_at_aphelion = elliptic.propagate_to_anomaly(np.pi * u.rad)
r_ap, _ = elliptic_at_aphelion.rv()
assert_quantity_allclose(norm(r_per), a * (1.0 - ecc))
assert_quantity_allclose(norm(r_ap), a * (1.0 + ecc))
# TODO: Test specific values
assert elliptic_at_perihelion.epoch > elliptic.epoch
assert elliptic_at_aphelion.epoch > elliptic.epoch
# test 10 random true anomaly values
# TODO: Rework this test
for nu in np.random.uniform(low=-np.pi, high=np.pi, size=10):
elliptic = elliptic.propagate_to_anomaly(nu * u.rad)
r, _ = elliptic.rv()
assert_quantity_allclose(norm(r), a * (1.0 - ecc ** 2) / (1 + ecc * np.cos(nu)))
def test_propagation_hyperbolic():
# Data from Curtis, example 3.5
r0 = [Earth.R.to(u.km).value + 300, 0, 0] * u.km
v0 = [0, 15, 0] * u.km / u.s
expected_r_norm = 163180 * u.km
expected_v_norm = 10.51 * u.km / u.s
ss0 = Orbit.from_vectors(Earth, r0, v0)
tof = 14941 * u.s
ss1 = ss0.propagate(tof)
r, v = ss1.rv()
assert_quantity_allclose(norm(r), expected_r_norm, rtol=1e-4)
assert_quantity_allclose(norm(v), expected_v_norm, rtol=1e-3)
), "Body 1 is not the main body: it has less mass than the 'secondary' body"
# Define local frame of reference:
# Center: main body, NOT the barycenter
# x-axis: points to the secondary body
ux = r2 - r1
r12 = norm(ux)
ux = ux / r12
# y-axis: contained in the orbital plane, perpendicular to x-axis
def cross(x, y):
return np.cross(x, y) * x.unit * y.unit
uy = cross(n, ux)
uy = uy / norm(uy)
# position in x-axis
x1, x2, x3, x4, x5 = lagrange_points(r12, m1, m2)
# position in y-axis
# L1, L2, L3 are located in the x-axis, so y123 = 0
# L4 and L5 are points in the plane of rotation which form an equilateral
# triangle with the two masses (Battin)
# sqrt(3)/2 = sin(60 deg)
y4 = np.sqrt(3) / 2 * r12
y5 = -y4
# Convert L points coordinates (x,y) to original vectorial base [r1 r2]
L1 = r1 + ux * x1
L2 = r1 + ux * x2
num_contours: int
Number of contour lines for flyby speed
N: int
Number of points for flyby angle
Note
----
The algorithm for generating Tisserand plots is the one depicted in
"Preliminary Trajectory Design of a Mission to Enceladus" by David
Falcato Fialho Palma, section 3.6
"""
# Generate mean orbital elements Earth
body_rv = get_mean_elements(body).to_vectors()
R_body, V_body = norm(body_rv.r), norm(body_rv.v)
# Generate non-dimensional velocity and alpha span
vinf_array = np.linspace(vinf_span[0], vinf_span[-1], num_contours)
alpha_array = np.linspace(alpha_lim[0], alpha_lim[-1], N)
vinf_array /= V_body
# Construct the mesh for any configuration
V_INF, ALPHA = np.meshgrid(vinf_array, alpha_array)
# Solving for non-dimensional a_sc and ecc_sc
A_SC = 1 / np.abs(1 - V_INF ** 2 - 2 * V_INF * np.cos(ALPHA))
ECC_SC = np.sqrt(1 - 1 / A_SC * ((3 - 1 / A_SC - V_INF ** 2) / (2)) ** 2)
# Compute main Tisserand variables
RR_P = A_SC * R_body * (1 - ECC_SC)
RR_A = A_SC * R_body * (1 + ECC_SC)
ss_dpt = Orbit.from_vectors(attractor, rr_dpt_body, vv_dpt_body, epoch=t_launch)
ss_arr = Orbit.from_vectors(attractor, rr_arr_body, vv_arr_body, epoch=t_arrival)
# Define time of flight
tof = ss_arr.epoch - ss_dpt.epoch
if tof <= 0:
return None, None, None, None, None
try:
# Lambert is now a Maneuver object
man_lambert = Maneuver.lambert(ss_dpt, ss_arr)
# Get norm delta velocities
dv_dpt = norm(man_lambert.impulses[0][1])
dv_arr = norm(man_lambert.impulses[1][1])
# Compute all the output variables
c3_launch = dv_dpt ** 2
c3_arrival = dv_arr ** 2
return (
dv_dpt.to(u.km / u.s).value,
dv_arr.to(u.km / u.s).value,
c3_launch.to(u.km ** 2 / u.s ** 2).value,
c3_arrival.to(u.km ** 2 / u.s ** 2).value,
tof.jd,
)
except AssertionError:
return None, None, None, None, None
T_vec = np.cross(S_vec, c_vec) * u.one
T_vec = T_vec / norm(T_vec)
R_vec = np.cross(S_vec, T_vec) * u.one
# This vector defines the B-Plane
B_vec = b * (np.cos(theta) * T_vec + np.sin(theta) * R_vec)
# We have to rotate the inbound hyperbolic excess velocity
# an angle delta (turn angle) around a vector that is orthogonal to
# the B-Plane and trajectory plane
rot_v = np.cross(B_vec / b, S_vec) * u.one
# And now we rotate the outbound hyperbolic excess velocity
# u_vec = v_inf_1 / norm(v_inf) = S_vec
v_vec = np.cross(rot_v, v_inf_1) * u.one
v_vec = v_vec / norm(v_vec)
v_inf_2 = v_inf * (np.cos(delta) * S_vec + np.sin(delta) * v_vec)
v_spacecraft_out = v_inf_2 + v_body
return v_spacecraft_out, delta
v_inf = norm(v_inf_1)
ecc = 1 + r_p * v_inf ** 2 / k # Eccentricity of the entry hyperbola
delta = 2 * np.arcsin(1 / ecc) # Turn angle
b = k / v_inf ** 2 * np.sqrt(ecc ** 2 - 1) # Magnitude of the B vector
# Now we compute the unit vectors in which to return the outbound hyperbolic excess velocity:
# * S goes along the hyperbolic excess velocity and is perpendicular to the B-Plane,
# * T goes along the B-Plane and is parallel to _some_ fundamental plane - in this case, the plane in which
# the velocities are computed
# * R completes the orthonormal set
S_vec = v_inf_1 / v_inf
c_vec = np.array([0, 0, 1]) * u.one
T_vec = np.cross(S_vec, c_vec) * u.one
T_vec = T_vec / norm(T_vec)
R_vec = np.cross(S_vec, T_vec) * u.one
# This vector defines the B-Plane
B_vec = b * (np.cos(theta) * T_vec + np.sin(theta) * R_vec)
# We have to rotate the inbound hyperbolic excess velocity
# an angle delta (turn angle) around a vector that is orthogonal to
# the B-Plane and trajectory plane
rot_v = np.cross(B_vec / b, S_vec) * u.one
# And now we rotate the outbound hyperbolic excess velocity
# u_vec = v_inf_1 / norm(v_inf) = S_vec
v_vec = np.cross(rot_v, v_inf_1) * u.one
v_vec = v_vec / norm(v_vec)
v_inf_2 = v_inf * (np.cos(delta) * S_vec + np.sin(delta) * v_vec)
def _lambert(k, r1, r2, tof, M, numiter, rtol):
# Check preconditions
assert tof > 0
assert k > 0
# Chord
c = r2 - r1
c_norm, r1_norm, r2_norm = norm(c), norm(r1), norm(r2)
# Semiperimeter
s = (r1_norm + r2_norm + c_norm) * .5
# Versors
i_r1, i_r2 = r1 / r1_norm, r2 / r2_norm
i_h = np.cross(i_r1, i_r2)
i_h = i_h / norm(i_h) # Fixed from paper
# Geometry of the problem
ll = np.sqrt(1 - c_norm / s)
if i_h[2] < 0:
ll = -ll
i_h = -i_h
def set_frame(self, p_vec, q_vec, w_vec):
"""Sets perifocal frame.
Raises
------
ValueError
If the vectors are not a set of mutually orthogonal unit vectors.
"""
if not np.allclose([norm(v) for v in (p_vec, q_vec, w_vec)], 1):
raise ValueError("Vectors must be unit.")
elif not np.allclose([p_vec.dot(q_vec), q_vec.dot(w_vec), w_vec.dot(p_vec)], 0):
raise ValueError("Vectors must be mutually orthogonal.")
else:
self._frame = p_vec, q_vec, w_vec
if self._orbits:
self._redraw()