Inspecting the CR3BP and BCP#
# our stuff
import pykep as pk
import heyoka as hy
# their stuff
import numpy as np
# plotting
import matplotlib.pyplot as plt
%matplotlib ipympl
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
File ~/local/lib/python3.13/site-packages/matplotlib/backends/registry.py:407, in BackendRegistry.resolve_gui_or_backend(self, gui_or_backend)
406 try:
--> 407 return self.resolve_backend(gui_or_backend)
408 except Exception: # KeyError ?
File ~/local/lib/python3.13/site-packages/matplotlib/backends/registry.py:369, in BackendRegistry.resolve_backend(self, backend)
368 if gui is None:
--> 369 raise RuntimeError(f"'{backend}' is not a recognised backend name")
371 return backend, gui if gui != "headless" else None
RuntimeError: 'ipympl' is not a recognised backend name
During handling of the above exception, another exception occurred:
RuntimeError Traceback (most recent call last)
Cell In[1], line 10
8 # plotting
9 import matplotlib.pyplot as plt
---> 10 get_ipython().run_line_magic('matplotlib', 'ipympl')
File ~/local/lib/python3.13/site-packages/IPython/core/interactiveshell.py:2488, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2486 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2487 with self.builtin_trap:
-> 2488 result = fn(*args, **kwargs)
2490 # The code below prevents the output from being displayed
2491 # when using magics with decorator @output_can_be_silenced
2492 # when the last Python token in the expression is a ';'.
2493 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/local/lib/python3.13/site-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
98 print(
99 "Available matplotlib backends: %s"
100 % _list_matplotlib_backends_and_gui_loops()
101 )
102 else:
--> 103 gui, backend = self.shell.enable_matplotlib(args.gui)
104 self._show_matplotlib_backend(args.gui, backend)
File ~/local/lib/python3.13/site-packages/IPython/core/interactiveshell.py:3760, in InteractiveShell.enable_matplotlib(self, gui)
3757 import matplotlib_inline.backend_inline
3759 from IPython.core import pylabtools as pt
-> 3760 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
3762 if gui != None:
3763 # If we have our first gui selection, store it
3764 if self.pylab_gui_select is None:
File ~/local/lib/python3.13/site-packages/IPython/core/pylabtools.py:349, in find_gui_and_backend(gui, gui_select)
347 else:
348 gui = _convert_gui_to_matplotlib(gui)
--> 349 backend, gui = backend_registry.resolve_gui_or_backend(gui)
351 gui = _convert_gui_from_matplotlib(gui)
352 return gui, backend
File ~/local/lib/python3.13/site-packages/matplotlib/backends/registry.py:409, in BackendRegistry.resolve_gui_or_backend(self, gui_or_backend)
407 return self.resolve_backend(gui_or_backend)
408 except Exception: # KeyError ?
--> 409 raise RuntimeError(
410 f"'{gui_or_backend}' is not a recognised GUI loop or backend name")
RuntimeError: 'ipympl' is not a recognised GUI loop or backend name
We start by importing orbital propagators for the Circular Restricted Three-Body Problem (CR3BP) and the Earth-Moon system. We do the same for the Bicircular Problem.
# We get the Taylor adaptive integrator for the Circular Restricted 3-body Problem
ta_cr3bp = pk.ta.get_cr3bp(tol=1e-16)
# We get the Taylor adaptive integrator for the Bi-Circular Problem
ta_bcp = pk.ta.get_bcp(tol=1e-16)
# We set the parameters for the CR3BP (Earth-Moon)
CR3BP_MU = pk.CR3BP_MU_EARTH_MOON
ta_cr3bp.pars[:] = [CR3BP_MU]
# We set the parameters for the Bi-Circular Problem (Earth-Moon-Sun)
ta_bcp.pars[:] = [CR3BP_MU, pk.BCP_MU_S, pk.BCP_RHO_S, pk.BCP_OMEGA_S]
We also define functions to compute the Jacobi constant and the effective potential in the CR3BP.
state = hy.make_vars("x", "y", "z", "vx", "vy", "vz")
jacobi_C = hy.cfunc([pk.ta.cr3bp_jacobi_C()], vars=state)
effective_potential_U = hy.cfunc([pk.ta.cr3bp_effective_potential_U()], vars=state[:3])
Plotting the zero velocity curves#
xx = np.linspace(-1.9, 1.9, 300)
yy = np.linspace(-1.9, 1.9, 300)
x_grid, y_grid = np.meshgrid(xx, yy)
z_grid = np.zeros(np.shape(x_grid))
x_flat = x_grid.flatten()
y_flat = y_grid.flatten()
z_flat = z_grid.flatten()
U_grid = effective_potential_U(
np.vstack((x_flat, y_flat, z_flat)), pars=[[CR3BP_MU] * len(x_flat)]
).reshape(np.shape(x_grid))
fig, ax = plt.subplots()
zv = ax.contour(
x_grid,
y_grid,
2*U_grid,
levels=[3.08, 3.1, 3.12, 3.14, 3.16],
cmap="viridis",
)
ax.set_aspect("equal")
Which we cal also make interactive to add a modern twist …
from ipywidgets import interact
from ipywidgets import FloatSlider, Layout
fig, ax = plt.subplots()
xx = np.linspace(-1.9, 1.9, 300)
yy = np.linspace(-1.9, 1.9, 300)
x_grid, y_grid = np.meshgrid(xx, yy)
z_grid = np.zeros(np.shape(x_grid))
x_flat = x_grid.flatten()
y_flat = y_grid.flatten()
z_flat = z_grid.flatten()
def update_zero_velocity_plot(jacobi):
ax.clear() # <-- this line clears previous contours
U_grid = effective_potential_U(
np.vstack((x_flat, y_flat, z_flat)),
pars=[[CR3BP_MU] * len(x_flat)],
).reshape(np.shape(x_grid))
ax.contour(x_grid, y_grid, 2*U_grid, levels=[jacobi], cmap="viridis")
ax.set_title(f"Jacobi constant: {jacobi:.3f}")
fig.canvas.draw() # force redraw in some environments
ax.set_aspect("equal")
interact(
update_zero_velocity_plot,
jacobi=FloatSlider(
value=1.57*2,
min=1.49*2,
max=1.64*2,
step=0.002,
readout_format=".3f",
layout=Layout(width="900px"),
),
)
<function __main__.update_zero_velocity_plot(jacobi)>
Generating some orbits#
Let us generate some orbits in the CR3BP starting from a circular LEO orbit. We do this “manually” shooting from the corresponding initial conditions and trying to find a trajectory actually intercepting the Moon.
EARTH_MOON_DISTANCE = 384400000
# Let us set the initial conditions in the Earth-Moon system. An equatorial LEO at an altitude of ~1700 km
r = 8000000 / EARTH_MOON_DISTANCE
# In the reference frame of the CR3BP, the initial coordinate is thus:
x = -CR3BP_MU-r
# The initial velocity is obtained from the circular velocity at the distance r from the center of mass.
vy = -np.sqrt((1-CR3BP_MU)/r)+x
# We prform the numerical integration.
ta_cr3bp.time = 0.0
ta_cr3bp.state[:] = [x, 0, 0, 0.0, vy-2.705, 0.0]
ta_bcp.time = 0.0
ta_bcp.state[:] = ta_cr3bp.state
C = jacobi_C(ta_cr3bp.state, pars=[CR3BP_MU])
t_grid = np.linspace(0, 4*np.pi, 1000)
sol_cr3bp = ta_cr3bp.propagate_grid(t_grid)
sol_bcp = ta_bcp.propagate_grid(t_grid)
print("Value of the Jacobi constant C:", C)
Value of the Jacobi constant C: [2.26785855]
We plot the orbit in the CR3BP frame as well in the (more familiar) derotated frame …
def plot(sol):
fig, ax = plt.subplots()
# In case we also plot the zero velocity curves, for many IC they will be completely opened and not show.
zv = ax.contourf(
x_grid,
y_grid,
2*U_grid,
levels=[0.1,C[0]],
colors='green'
)
ax.set_aspect("equal")
# The x,y in the rotating frame
xr = sol[-1][:, 0]
yr = sol[-1][:, 1]
# The x,y, in the derotated frame
x =xr*np.cos(t_grid) - yr*np.sin(t_grid)
y =xr*np.sin(t_grid) + yr*np.cos(t_grid)
ax.plot(xr, yr, "k-", label="CR3BP trajectory")
ax.plot(x,y, "r-", label="CR3BP trajectory")
ax.plot(-CR3BP_MU,0, "bo", label="Earth")
ax.plot(1-CR3BP_MU,0, "yo", label="Moon")
return fig, ax
First without the Sun perturbation ….
plot(sol_cr3bp);
… then with the Sun perturbation. We can see how a fly-by still happens but the overall geometry is significantly changed.
plot(sol_bcp);