Lambert Solver¶
Given an initial and target position vectors, a desired time-of-flight, as well as the gravitational parameter of the host planet, godot's lambert solver calculates--by default--the shortest elliptic arc that solves the equations of the two-body problem between the initial and target states.
Optionally, godot users can force the Lambert solver to calculate the clockwise or counter-clockwise arc or to account for multiple revolutions.
This tutorial illustrates some of the key capabilities of godot's lambert solver
We start by importing linear algebra and plotting libraries, along with the godot astro module
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import godot.core.astro as astro
Next, we define the initial and target orbits. As an example, we consider a bi-impulsive transfer from a nearly circular Sun-synchronous orbit (a = 7000, e = 0.01, I = 97.5 deg, RAAN = 15 deg, argP = 0 deg) to a piece of debris orbiting at higher altitudes (a = 8000, e = 0.1, I = 99.23 deg, RAAN = 18 deg, argP = 23 deg). The spacecraft fires its thrusters at the periapsis of its initial orbit.
# Earth gravitational parameter
GM = 398600.4418
# conversion factor from degrees 2 radians
d2r = np.pi/180
# Initial Orbit at periapsis
coe1 = np.array([7000.0, 0.01, 97.5*d2r, 15*d2r, 0.0, 0.0])
xyz1 = astro.cartFromKep(coe1, GM)
r1 = xyz1[0:3]
v1 = xyz1[3:6]
# Target Orbit without true anomaly
coe2 = np.array([7123.0, 0.1, 99.23*d2r, 18*d2r, 23*d2r])
So far, we have not specified the location within the target orbit for which the impulsive transfer needs to be computed, as well as the total time-of-flight of the transfer.
To this end, we create a two-dimensional grid of equally spaced points in the domain ta2: [0, 2*pi] x tof: [600, 6600] and search for near-optimal bi-impulsive burns between the initial and target orbit. This is where godot's lambert solver is first invoked. Both clockwise and counter-clockwise are calculated here, by simply changing the corresponding enum.
# Create a 2D grid of equally spaced points in the true anomaly x time-of-flight domains
nta, ntof = (500, 500)
dta2 = 2*np.pi*np.linspace(0, 1, nta)
dtof = 600 + 6000*np.linspace(0, 1, ntof)
ta2, tof = np.meshgrid(dta2, dtof)
# Initialize Arrays
long_arc = True
dv1_la = np.zeros(ta2.shape)
dv2_la = np.zeros(ta2.shape)
dv1_sa = np.zeros(ta2.shape)
dv2_sa = np.zeros(ta2.shape)
for i in range(nta):
for j in range(ntof):
# append arrival true anomaly to coe2 and convert to cartesian coordinates (pos & vel)
coe2t = np.append(coe2, ta2[i,j])
xyz2 = astro.cartFromKep(coe2t, GM)
r2 = xyz2[0:3]
v2 = xyz2[3:6]
branch = astro.LambertBranch.Left
# Lambert solver - long arc
sol_la = astro.lambert(r1, r2, tof[i,j], GM, 0, branch, astro.LambertDirection.Clockwise)
# calculate initial and final maneuvre costs for the short arc transfer
dv1_la[i, j] = np.linalg.norm(sol_la.v1 - v1)
dv2_la[i, j] = np.linalg.norm(v2 - sol_la.v2)
# Lambert solver - short arc
sol_sa = astro.lambert(r1, r2, tof[i,j], GM, 0, branch, astro.LambertDirection.CounterClockwise)
# calculate initial and final maneuvre costs for the short arc transfer
dv1_sa[i, j] = np.linalg.norm(sol_sa.v1 - v1)
dv2_sa[i, j] = np.linalg.norm(v2 - sol_sa.v2)
Having calculated the impulsive burns, we can analyse the data and find minimum fuel clockwise and counter-clockwise solutions
dvt_la = dv1_la + dv2_la
idx_la = np.where(dvt_la == np.min(dvt_la))
dvt_la_min = dvt_la[idx_la]
dvt_sa = dv1_sa + dv2_sa
idx_sa = np.where(dvt_sa == np.min(dvt_sa))
dvt_sa_min = dvt_sa[idx_sa]
Next, we can compare the results and determine the most fuel-efficient transfer
if dvt_la_min < dvt_sa_min:
print('Long arc solution is more fuel efficient')
print('delta-v tot: ', dvt_la_min, ' (km/s)')
long_arc = True
ta2_min = ta2[idx_la]
tof_min = tof[idx_la]
else:
print('Short arc solution is more fuel efficient')
print('delta-v tot: ', dvt_sa_min, ' (km/s)')
long_arc = False
ta2_min = ta2[idx_sa]
tof_min = tof[idx_sa]
print(ta2_min, tof_min, long_arc)
At last, we can make some plots to display the total delta-v contour lines as well as the location of the minimum-fuel clockwise and counter-clockwise solutions. In this case, since the orbits are retrograde (inc > 90 deg), the optimal solution will always be the clockwise.
fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
lapc = ax1.contour(ta2, tof, dv1_la + dv2_la, np.linspace(0, 8, 21))
ax1.plot(ta2[idx_la], tof[idx_la], '*r')
ax1.set_xlabel('Arrival True Anomaly (rad)')
ax1.set_ylabel('Time of flight (s)')
ax1.title.set_text('Long Arcs')
sapc = ax2.contour(ta2, tof, dv1_sa + dv2_sa, np.linspace(0, 8, 21))
ax2.plot(ta2[idx_sa], tof[idx_sa], '*r')
ax2.set_xlabel('Arrival True Anomaly (rad)')
ax2.set_ylabel('Time of flight (s)')
ax2.title.set_text('Short Arcs')
cbar = plt.colorbar(sapc, ax=ax2)
cbar.set_label('Total delta-v (km/s)')
Recycling the results of the comparative analysis, we can assess the performance and accuracy of godot's Lambert Solver by re-calculating the transfer for the identified true anomaly at arrival and minimum-delta v time-of-flight:
# Set the true anomaly at arrival
coe2f = np.append(coe2, ta2_min)
xyz2 = astro.cartFromKep(coe2f, GM)
r2 = xyz2[0:3]
v2 = xyz2[3:6]
# Solve the Lambert Problem from r0 to r1 in time TOF=tof_min
sol = astro.lambert(r1, r2, tof_min, GM, 0, branch, astro.LambertDirection.Clockwise)
print(sol.a)
# print(sol.p)
print(sol.v1)
print(sol.v2)
print('dv1: ', sol.v1 - v1)
print('dv2: ', v2 - sol.v2)
print('Total delta-v: ', np.linalg.norm(sol.v1 - v1) + np.linalg.norm(v2 - sol.v2), ' (km/s)')
We can also make same plots to double-check the appropriatedeness of the solution.
First, we need to propagate the spacecraft states over time. We use godot's astro.keplerPropagate to propagate the initial, transfer, and final orbits of the spacecraft using Kepler's equation. The method astro.cartFromKep is used in parallel to convert from keplerian orbit elements to cartesian coordinates and store the spacecraft's xyz coordinates for plotting purposes
N = 1000
dt = np.linspace(0, 1, N)
# calculate initial orbit's period
n1 = np.sqrt(GM/coe1[0]**3)
P1 = 2*np.pi/n1
# calculate target orbit's period
n2 = np.sqrt(GM/coe2[0]**3)
P2 = 2*np.pi/n2
# convert transfer orbit from Cart to Kep
coet = astro.kepFromCart(np.concatenate((r1, sol.v1)), GM)
# propagate over time to generate trajectories
X1t = np.zeros((N, 6))
Xtt = np.zeros((N, 6))
X2t = np.zeros((N, 6))
for idx, t in enumerate(dt):
# propagate orbits using Kepler's equation
coe1t = astro.keplerPropagate(coe1, GM, P1*t)
coett = astro.keplerPropagate(coet, GM, tof_min*t)
coe2t = astro.keplerPropagate(coe2f, GM, P2*t)
# convert and store cartesian coordinates
X1t[idx] = astro.cartFromKep(coe1t, GM)
Xtt[idx] = astro.cartFromKep(coett, GM)
X2t[idx] = astro.cartFromKep(coe2t, GM)
3D Plots of the spacecraft trajectories can be now added
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot3D(X1t[:, 0], X1t[:, 1], X1t[:, 2], '--k')
ax.plot3D(X1t[0, 0], X1t[0, 1], X1t[0, 2], 'ok')
ax.plot3D(X2t[:, 0], X2t[:, 1], X2t[:, 2], '--r')
ax.plot3D(X2t[0, 0], X2t[0, 1], X2t[0, 2], 'or')
ax.plot3D(Xtt[:, 0], Xtt[:, 1], Xtt[:, 2], 'b')
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Z (km)')
ax.set_xlim(-1e4, 1e4)
ax.set_ylim(-1e4, 1e4)
ax.set_zlim(-1e4, 1e4)
ax.view_init(30, 30)
At last, we can assess the accuracy of the method by comparing the final state of the transfer trajectory with the target state:
# Check position and velocity errors
err_pos = Xtt[-1, 0:3] - r2
print('Pos. Error: ', np.linalg.norm(err_pos))