# Analyzing the Parker Solar Probe flybys¶

## 1. Modulus of the exit velocity, some features of Orbit #2¶

First, using the data available in the reports, we try to compute some of the properties of orbit #2. This is not enough to completely define the trajectory, but will give us information later on in the process.

[1]:

from astropy import units as u

[2]:

T_ref = 150 * u.day
T_ref

[2]:

$150 \; \mathrm{d}$
[3]:

from poliastro.bodies import Earth, Sun, Venus

[4]:

k = Sun.k
k

[4]:

$1.3271244 \times 10^{20} \; \mathrm{\frac{m^{3}}{s^{2}}}$
[5]:

import numpy as np

$T = 2 \pi \sqrt{\frac{a^3}{\mu}} \Rightarrow a = \sqrt[3]{\frac{\mu T^2}{4 \pi^2}}$
[6]:

a_ref = np.cbrt(k * T_ref**2 / (4 * np.pi**2)).to(u.km)
a_ref.to(u.au)

[6]:

$0.55249526 \; \mathrm{AU}$
$\varepsilon = -\frac{\mu}{r} + \frac{v^2}{2} = -\frac{\mu}{2a} \Rightarrow v = +\sqrt{\frac{2\mu}{r} - \frac{\mu}{a}}$
[7]:

energy_ref = (-k / (2 * a_ref)).to(u.J / u.kg)
energy_ref

[7]:

$-8.0283755 \times 10^{8} \; \mathrm{\frac{J}{kg}}$
[8]:

from poliastro.twobody import Orbit
from poliastro.util import norm

from astropy.time import Time

[9]:

flyby_1_time = Time("2018-09-28", scale="tdb")
flyby_1_time

[9]:

<Time object: scale='tdb' format='iso' value=2018-09-28 00:00:00.000>

[10]:

r_mag_ref = norm(Orbit.from_body_ephem(Venus, epoch=flyby_1_time).r)
r_mag_ref.to(u.au)

[10]:

$0.72573132 \; \mathrm{AU}$
[11]:

v_mag_ref = np.sqrt(2 * k / r_mag_ref - k / a_ref)
v_mag_ref.to(u.km / u.s)

[11]:

$28.967364 \; \mathrm{\frac{km}{s}}$

## 2. Lambert arc between #0 and #1¶

To compute the arrival velocity to Venus at flyby #1, we have the necessary data to solve the boundary value problem.

[12]:

d_launch = Time("2018-08-11", scale="tdb")
d_launch

[12]:

<Time object: scale='tdb' format='iso' value=2018-08-11 00:00:00.000>

[13]:

ss0 = Orbit.from_body_ephem(Earth, d_launch)
ss1 = Orbit.from_body_ephem(Venus, epoch=flyby_1_time)

[14]:

tof = flyby_1_time - d_launch

[15]:

from poliastro import iod

[16]:

(v0, v1_pre), = iod.lambert(Sun.k, ss0.r, ss1.r, tof.to(u.s))

[17]:

v0

[17]:

$[9.5993373,~11.298552,~2.9244933] \; \mathrm{\frac{km}{s}}$
[18]:

v1_pre

[18]:

$[-16.980821,~23.307528,~9.1312908] \; \mathrm{\frac{km}{s}}$
[19]:

norm(v1_pre)

[19]:

$30.248465 \; \mathrm{\frac{km}{s}}$

## 3. Flyby #1 around Venus¶

We compute a flyby using poliastro with the default value of the entry angle, just to discover that the results do not match what we expected.

[20]:

from poliastro.threebody.flybys import compute_flyby

[21]:

V = Orbit.from_body_ephem(Venus, epoch=flyby_1_time).v
V

[21]:

$[648499.74,~2695078.4,~1171563.7] \; \mathrm{\frac{km}{d}}$
[22]:

h = 2548 * u.km

[23]:

d_flyby_1 = Venus.R + h
d_flyby_1.to(u.km)

[23]:

$8599.8 \; \mathrm{km}$
[24]:

V_2_v_, delta_ = compute_flyby(v1_pre, V, Venus.k, d_flyby_1)

[25]:

norm(V_2_v_)

[25]:

$27.755339 \; \mathrm{\frac{km}{s}}$

## 4. Optimization¶

Now we will try to find the value of $$\theta$$ that satisfies our requirements.

[26]:

def func(theta):
V_2_v, _ = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta * u.rad)
ss_1 = Orbit.from_vectors(Sun, ss1.r, V_2_v, epoch=flyby_1_time)
return (ss_1.period - T_ref).to(u.day).value


There are two solutions:

[27]:

import matplotlib.pyplot as plt

[28]:

theta_range = np.linspace(0, 2 * np.pi)
plt.plot(theta_range, [func(theta) for theta in theta_range])
plt.axhline(0, color='k', linestyle="dashed")

[28]:

<matplotlib.lines.Line2D at 0x7f5f3bddd2b0>

[29]:

func(0)

[29]:

-9.142672330001195

[30]:

func(1)

[30]:

7.09811543934556

[31]:

from scipy.optimize import brentq

[32]:

theta_opt_a = brentq(func, 0, 1) * u.rad
theta_opt_a.to(u.deg)

[32]:

$38.598709 \; \mathrm{{}^{\circ}}$
[33]:

theta_opt_b = brentq(func, 4, 5) * u.rad
theta_opt_b.to(u.deg)

[33]:

$279.3477 \; \mathrm{{}^{\circ}}$
[34]:

V_2_v_a, delta_a = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta_opt_a)
V_2_v_b, delta_b = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta_opt_b)

[35]:

norm(V_2_v_a)

[35]:

$28.967364 \; \mathrm{\frac{km}{s}}$
[36]:

norm(V_2_v_b)

[36]:

$28.967364 \; \mathrm{\frac{km}{s}}$

## 5. Exit orbit¶

And finally, we compute orbit #2 and check that the period is the expected one.

[37]:

ss01 = Orbit.from_vectors(Sun, ss1.r, v1_pre, epoch=flyby_1_time)
ss01

[37]:

0 x 1 AU x 18.8 deg (HCRS) orbit around Sun (☉)


The two solutions have different inclinations, so we still have to find out which is the good one. We can do this by computing the inclination over the ecliptic - however, as the original data was in the International Celestial Reference Frame (ICRF), whose fundamental plane is parallel to the Earth equator of a reference epoch, we have change the plane to the Earth ecliptic, which is what the original reports use.

[38]:

ss_1_a = Orbit.from_vectors(Sun, ss1.r, V_2_v_a, epoch=flyby_1_time)
ss_1_a

[38]:

0 x 1 AU x 25.0 deg (HCRS) orbit around Sun (☉)

[39]:

ss_1_b = Orbit.from_vectors(Sun, ss1.r, V_2_v_b, epoch=flyby_1_time)
ss_1_b

[39]:

0 x 1 AU x 13.1 deg (HCRS) orbit around Sun (☉)


Let’s define a function to do that quickly for us, using the get_frame <https://docs.poliastro.space/en/latest/safe.html#poliastro.frames.get_frame>__ function from poliastro.frames:

[40]:

from astropy.coordinates import CartesianRepresentation
from poliastro.frames import Planes, get_frame

def change_plane(ss_orig, plane):
"""Changes the plane of the Orbit.

"""
ss_orig_rv = ss_orig.frame.realize_frame(
ss_orig.represent_as(CartesianRepresentation)
)

dest_frame = get_frame(ss_orig.attractor, plane, obstime=ss_orig.epoch)

ss_dest_rv = ss_orig_rv.transform_to(dest_frame)
ss_dest_rv.representation_type = CartesianRepresentation

ss_dest = Orbit.from_vectors(
ss_orig.attractor,
r=ss_dest_rv.data.xyz,
v=ss_dest_rv.data.differentials['s'].d_xyz,
epoch=ss_orig.epoch,
plane=plane,
)
return ss_dest

[41]:

change_plane(ss_1_a, Planes.EARTH_ECLIPTIC)

[41]:

0 x 1 AU x 3.5 deg (HeliocentricEclipticJ2000) orbit around Sun (☉)

[42]:

change_plane(ss_1_b, Planes.EARTH_ECLIPTIC)

[42]:

0 x 1 AU x 13.1 deg (HeliocentricEclipticJ2000) orbit around Sun (☉)


Therefore, the correct option is the first one.

[43]:

ss_1_a.period.to(u.day)

[43]:

$150 \; \mathrm{d}$
[44]:

ss_1_a.a

[44]:

$82652115 \; \mathrm{km}$

And, finally, we plot the solution:

[45]:

from poliastro.plotting import OrbitPlotter

[46]:

frame = OrbitPlotter()

frame.plot(ss0, label=Earth)
frame.plot(ss1, label=Venus)
frame.plot(ss01, label="#0 to #1")
frame.plot(ss_1_a, label="#1 to #2");

/home/juanlu/Personal/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:434: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned