Generador de gráficos de chuleta de cerdo y casos de prueba del solucionador de Lambert

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.PCP Tierra a Marte

Respuestas (1)

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:

METRO 0 + m a 3 Δ t = mi mi pecado ( mi )

para una amplia variedad de ( mi , a , Δ t ) 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).

Transferencia Tierra-Marte Porkchop

¿Puede proporcionar un enlace a una "herramienta gratuita"? ¿Qué ha pedido el OP en la pregunta? ¿O unas pocas líneas de código, un guión, un esquema de un algoritmo, o al menos un enlace a un algoritmo? No todo el mundo tiene a mano un solucionador de ecuaciones trascendentales. En este momento, esto en realidad no ayuda al OP, solo dice lo que haría usted mismo si estuviera en una situación similar. Intenta ponerte en la posición del OP.
No necesita un "solucionador de ecuaciones trascendentales" a mano. Puede hacer el cálculo explícitamente porque la anomalía excéntrica es la variable independiente. En cualquier caso, adjunto un código.
Fantástico - ¡Nunca supe esto! si tu escoges a y mi y establecer m = 1 y paso mi de 0 a π , puedes obtener r y θ directamente, sin necesidad de un solucionador elegante o incluso el método de Newton . Me he caído de la silla. Si lo desea, también podría escribir una buena explicación para esto como una respuesta complementaria aquí . ¡Gracias por aclarármelo!
Cosa segura. Creo que la ética en torno a la ecuación de Kepler nos ha hecho creer a todos que siempre debemos resolverla para E . Me enamoro de eso muy a menudo.
Me siento terrible por señalar esto, pero si desea un punto por grado al incluir los puntos finales, necesita 361 puntos, el código existente avanza 1,00278 grados. (Estoy tratando de validar estrictamente mi propio código de validación contra su código de validación antes de intentar depurar mi propio Gooding Solver con errores...)