Estoy buscando una herramienta gratuita que pueda generar gráficos de chuletas de cerdo. Hago esto para comparar la salida de mi herramienta con la de un tercero. Además, si alguien tiene casos de prueba específicos para resolver el problema del valor límite de Lambert, estaría más que ansioso por echar un vistazo.
Esto es puramente para la validación.
La siguiente figura es un PCP de la Tierra a Marte que generé justo antes.
Utilizo un solucionador de problemas de Kepler para probar un solucionador de problemas de Lambert. Por ejemplo: si quiero generar una amplia variedad de casos de prueba donde vario la excentricidad y el semieje mayor, simplemente resuelvo:
para una amplia variedad de combinaciones Habiendo resuelto una multitud de problemas de Kepler (que puedo probar de inmediato, porque la ecuación de Kepler debe cumplirse) resuelvo el problema de Lambert correspondiente y compruebo que obtengo la misma anomalía excéntrica.
El siguiente código explica este enfoque.
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
# radial/tangential)
diff = np.sqrt((v[n1] - vs[0])**2 + (v[n2] - vs[1])**2)
if diff > max_diff:
max_diff = diff
# report the maximum difference
print max_diff
Observe cómo puede calcular el tiempo de vuelo dado un rango de anomalías excéntricas, excentricidad y semieje mayor. En este caso, estoy considerando solo una revolución, pero la generalización es sencilla.
Con la anomalía excéntrica, puede calcular la anomalía real y el tiempo de vuelo (el último a través de la ecuación de Kepler, pero no necesita resolver iterativamente ).
Esa información se utiliza para obtener alcance y velocidad en diferentes puntos de la órbita.
Luego puede usar todas esas combinaciones de rango, velocidad, tiempo de vuelo y ángulo de transferencia (diferencia en la anomalía real porque la órbita es Kepleriana) para probar su solucionador de Lambert.
Realicé una prueba rápida con mi propia implementación del algoritmo de Gooding y obtuve una diferencia máxima de 8.76111591946e-14
, perfectamente razonable para el cálculo de doble precisión.
A modo de comparación, generé un diagrama similar: se compara bien con el suyo. Le recomendaría aumentar los ejes en su tiempo de llegada a aproximadamente 420 días después del 1 de enero de 2016 para asegurarse de que su solucionador de Lambert detecte el nodo de transferencia (también, tenga en cuenta nuestra época de referencia diferente para nuestras fechas de llegada).
UH oh
Escualo
UH oh
Escualo
lamont