Pork chop plot generator and Lambert solver test cases

I am looking for a free tool which can generates pork chop plots. I'm doing this to compare the output of my tool with one of a third party. Additionally, if anyone has specific test cases on solving Lambert's boundary value problem, I would be more than eager to have a look.

This is purely for validation.

The following figure is an Earth to Mars PCP I generated just earlier. I use a Kepler's Problem solver to test a Lambert's Problem solver. For example: if I want to generate a wide variety of test cases where I vary eccentricity and semimajor axis, I simply solve:

$$M_0 + \sqrt{\frac{\mu}{a^3}}\Delta t = E - e \sin(E)$$

for a wide variety of $(e, a, \Delta t)$ combinations. Having solved a multitude of Kepler's problems (which I can test immediately, because Kepler's equation must hold) I solve the corresponding Lambert's problem and verify that I obtain the same eccentric anomaly.

The following code explains this approach

import numpy as np

def kepler(a, e, mu=1.0):
# EXPLICIT calculation of time of flight given eccentric anomaly
# [in 0, 2 pi]
n = np.sqrt(mu / pow(a, 3))
E = np.linspace(0, 2 * np.pi, 360) # one point per degree
t = (E - e * np.sin(E)) / n
return (t, E)

def cartesian(sma, ecc, ean, mu=1.0):
smp = sma * (1 - ecc**2)
energy = -mu / (2 * sma)
f = 2 * np.arctan(np.sqrt((1 + ecc) / (1 - ecc)) * np.tan(ean / 2))
r = smp / (1 + ecc * np.cos(f))
v = np.sqrt(2 * (energy + mu / r))
if f < 0:
f = 2 * np.pi + f
return np.array([r, v, f])

# Let's test our algorithm for a given combination of sma/ecc. In a
# real test you would add more values, and you would consider multiple
# rev (you will get two solutions for most of the multirev cases)
sma = 1.0
ecc = 0.1

# this computes time of flight and eccentric anomaly NOTICE THIS IS AN
# EXPLICIT CALCULATION
tof, ean = kepler(sma, ecc)

# this transforms the quantities to cartesian vector magnitudes and
# unrestricted true anomaly (to compute transfer angle)
rvf = np.array([cartesian(sma, ecc, E) for E in ean])

r = rvf[:, 0]
v = rvf[:, 1]
f = rvf[:, 2]

# Here I'll just call my own Lambert solver (implemented in C++)
import ctypes
lib = ctypes.CDLL("libsolve.so")
vs = (ctypes.c_double * 2)()
N = len(tof)
max_diff = 0.0

# This ugly loop simply compares every staged-pairwise combination for
# SINGLE REV (you can do multiple rev as well, you just need to
# increase TOF and f to be urestricted). You would attach your own
# Lambert solver here.
for n1 in xrange(N):
for n2 in xrange(n1 + 1, N):
lib.solve(ctypes.c_double(1.0), # gm
ctypes.c_double(r[n1]), # norm of first position
ctypes.c_double(r[n2]), # norm of second position
ctypes.c_double(f[n2] - f[n1]), # angle traveled
ctypes.c_double(tof[n2] - tof[n1]), # time of flight
ctypes.byref(vs)) # output

# compare the norm of velocities (you could decompose in
diff = np.sqrt((v[n1] - vs)**2 + (v[n2] - vs)**2)

if diff > max_diff:
max_diff = diff

# report the maximum difference
print max_diff

Notice how you can calculate the time of flight given a range of eccentric anomalies, eccentricity, and semimajor axis. In this case I'm considering single-rev only, but the generalization is straightforward.

With the eccentric anomaly you can compute true anomaly and time of flight (the later via Kepler's equation, but you do not need to solve iteratively).

That information is used to obtain range and speed at different points in the orbit.

You can then use all those combinations of range, speed, time of flight, and transfer angle (difference in true anomaly because the orbit is Keplerian) to test your Lambert solver.

I ran a quick test with my own implementation of Gooding's algorithm and obtained a maximum difference of 8.76111591946e-14, perfectly reasonable for double-precision calculation.

For comparison, I generated a similar diagram - it compares well with yours. I would recommend you increase the axes on your arrival time to about 420 days past Jan 01, 2016 to ensure that your Lambert solver detects the transfer node (also, please notice our different baseline epoch for our arrival dates). • Can you supply a link to a "free tool" - what the OP has asked for in the question? Or a few lines of code, or script, or an outline of an algorithm, or at least a link to an algorithm? Not everyone has a transcendental equation solver handy. Right now this doesn't actually help the OP - it just says what you would do yourself if you were in a similar situation. Try to put yourself in the OP's position. – uhoh Mar 18 '17 at 4:06
• You don't need a "trascendental equation solver" handy. You can do the calculation explicitly because eccentric anomaly is the independent variable. At any rate, I'm attaching some code. – Escualo Mar 20 '17 at 1:37
• Fantastic - I never knew this! if you choose $a$ and $e$ and set $\mu=1$ and step $E$ from $0$ to $\pi$, you can get $r$ and $\theta$ directly, without ever needing a fancy solver or even Newton's method. I have fallen off my chair. If you like, you could write a nice explanation for this as a supplemental answer here as well. Thank you for clearing that up for me! – uhoh Mar 20 '17 at 3:03
• Sure thing. I think the ethos around Kepler's equation has made us all believe that we always need to solve it for E. I fall for it very often. – Escualo Mar 20 '17 at 3:21
• I feel terrible for pointing this out, but if you want one point per degree when including the endpoints, you need 361 points, the existing code steps by 1.00278 degrees. (I'm trying to strictly validate my own validation code against your validation code before trying to debug my own buggy Gooding Solver...) – lamont Mar 10 at 20:10