En los comentarios, que ya estaban informados de que los datos numéricos se tenía que utilizar, así que voy a explicar a partir de ese punto.
La figura 8 diseño no es porque es una ruta óptima. Esto ocurre debido a la gravedad de la luna y la Tierra. Cuando la nave entra dentro de la esfera de influencia (SOI) de la luna, la nave espacial se tira hacia él. Si la nave espacial se mueve a la velocidad de escape, la luna va a perturbar el vuelo, pero la nave no una mosca. Con la velocidad actual, la gravedad de la luna es suficiente para causar que un orbital volar. Al salir de la luna, del SOL, la nave es de ser tragados por la Tierra. Dado que la trayectoria de la mosca fue lanzar la nave espacial lejos de la luna, la luna se cruza en su camino original, pero esta es de corta duración ya que la Tierra, a continuación, tira de nuevo. Si la nave hubiera recogido bastante la velocidad de la maniobra orbital a estar en una parabólica o hiperbólica trayectoria, podría haber escapado de la atracción de la Tierra y ha sido enviado al espacio.
Una forma de determinar las velocidades de los test en el diseño de un vuelo para encontrar la constante de Jacobi, $C$. Para un determinado $C$, el cero curvas de velocidad se determinan. Puesto que quería llegar a la luna, $C\geq -1.6649$ que corresponde a una velocidad inicial de, al menos,$10.85762$, pero con una velocidad de $11.01$ es la velocidad de escape de la Tierra de modo que la velocidad inicial tiene que ser menos de $v_{esc}$.
Podemos derivar las ecuaciones de movimiento sólo hasta un punto. Estoy saltando las básicas de derivación de los dos cuerpo a cuerpo problema por lo que nos podemos mover en el quid de la cuestión.
\begin{alignat*}{4}
r_{12} & = \sqrt{(x_1 - x_2)^2} & \qquad & x_1
&&{}= \text{Is the } x \text{ location of } m_1\\
x_2 & = x_1 + r_{12} & & &&{}\phantom{=}
\text{relative to the center of gravity.}\\
x_1 & = \frac{-m_2}{m_1 + m_2}r_{12} & & \pi_1 &&{}=
\frac{m_1}{m_1 + m_2}\\
x_2 & = \frac{m_1}{m_1 + m_2}r_{12} & & \pi_2 &&{}=
\frac{m_2}{m_1 + m_2}\\
0 & = m_1x_1 + m_2x_2
\end{alignat*}
Podemos describir la posición de $m$
$\mathbf{r} = x\hat{\mathbf{i}} + y\hat{\mathbf{j}} + z\hat{\mathbf{k}}$ en relación con el centro de gravedad, es decir, el origen.
\begin{align}
\mathbf{r}_1
& = (x - x_1)\hat{\mathbf{i}} + y\hat{\mathbf{j}} +
z\hat{\mathbf{k}}\\
& = (x + \pi_2r_{12})\hat{\mathbf{i}} + y\hat{\mathbf{j}}
+ z\hat{\mathbf{k}}\tag{1}\\
\mathbf{r}_2
& = (x - x_2)\hat{\mathbf{i}} + y\hat{\mathbf{j}} +
z\hat{\mathbf{k}}\\
& = (x - \pi_1r_{12})\hat{\mathbf{i}} + y\hat{\mathbf{j}}
+ z\hat{\mathbf{k}}\tag{2}
\end{align}
Vamos a definir la aceleración absoluta donde $\omega$ es la inicial de angular
velocidad es constante.
A continuación,$\omega = \frac{2\pi}{T}$.
$$
\ddot{\mathbf{r}}_{\text{abs}}
= \mathbf{a}_{\text{rel}} + \mathbf{a}_{\text{CG}}
+ \mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r})
+ \dot{\mathbf{\varOmega}}\times\mathbf{r}
+ 2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}}
\etiqueta{3}
$$
donde
\begin{alignat*}{4}
\mathbf{a}_{\text{rel}}
& = \text{Rectilinear acceleration relative to the frame} & \quad &
\mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r}) &&{}=
\text{Centripetal acceleration}\\
2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}} & =
\text{Coriolis acceleration}
\end{alignat*}
Puesto que la velocidad del centro de gravedad es constante,
$\mathbf{a}_{\text{CG}} = 0$, e $\dot{\mathbf{\varOmega}} = 0$ desde el
la velocidad angular de una órbita circular es constante.
Por lo tanto, $(3)$ se convierte en:
$$
\ddot{\mathbf{r}} = \mathbf{a}_{\text{rel}}
+ \mathbf{\varOmega}\times(\mathbf{\varOmega}\times\mathbf{r})
+ 2\mathbf{\varOmega}\times\mathbf{v}_{\text{rel}}
\etiqueta{4}
$$
donde
\begin{align}
\mathbf{\varOmega} & = \varOmega\hat{\mathbf{k}}\tag{5}\\
\mathbf{r} & = x\hat{\mathbf{i}} + y\hat{\mathbf{j}}
+ z\hat{\mathbf{k}}\tag{6}\\
\dot{\mathbf{r}} & = \mathbf{v}_{\text{CG}}
+ \mathbf{\varOmega}\times\mathbf{r}
+ \mathbf{v}_{\text{rel}}\tag{7}\\
\mathbf{v}_{\text{rel}} & = \dot{x}\hat{\mathbf{i}}
+ \dot{y}\hat{\mathbf{j}}
+ \dot{z}\hat{\mathbf{k}}\tag{8}\\
\mathbf{a}_{\text{rel}} & = \ddot{x}\hat{\mathbf{i}}
+ \ddot{y}\hat{\mathbf{j}}
+ \ddot{z}\hat{\mathbf{k}}\tag{9}
\end{align}
Después de la sustitución de $(5)$, $(6)$, $(8)$, y
$(9)$ a $(3)$, obtenemos
$$
\ddot{\mathbf{r}} =
\left(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\right)\hat{\mathbf{i}}
+ \left(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\right)\hat{\mathbf{j}}
+ \ddot{z}\hat{\mathbf{k}}.
\etiqueta{10}
$$
Newton $2^{\text{nd}}$ Ley del Movimiento es
$m\mathbf{a} = \mathbf{F}_1 + \mathbf{F}_2$ donde
$\mathbf{F}_1 = -\frac{Gm_1m}{r_1^3}\mathbf{r}_1$ y
$\mathbf{F}_2 = -\frac{Gm_2m}{r_2^3}\mathbf{r}_2$.
Deje $\mu_1 = Gm_1$$\mu_2 = Gm_2$.
\begin{align}
m\mathbf{a} & = \mathbf{F}_1 + \mathbf{F}_2\\
m\mathbf{a} & = -\frac{m\mu_1}{r_1^3}\mathbf{r}_1
- \frac{m\mu_2}{r_2^3}\mathbf{r}_2\\
\mathbf{a} & = -\frac{\mu_1}{r_1^3}\mathbf{r}_1
- \frac{\mu_2}{r_2^3}\mathbf{r}_2\\
\bigl(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\bigr)\hat{\mathbf{i}}
+ \bigl(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\bigr)\hat{\mathbf{j}}
+ \ddot{z}\hat{\mathbf{k}}
& = -\frac{\mu_1}{r_1^3}\mathbf{r}_1
- \frac{\mu_2}{r_2^3}\mathbf{r}_2\\
\bigl(\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x\bigr)\hat{\mathbf{i}}
+ \bigl(\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y\bigr)\hat{\mathbf{j}}
+ \ddot{z}\hat{\mathbf{k}}
& =
\begin{aligned}
- & \frac{\mu_1}{r_1^3}\Bigl[(x +
\pi_2r_{12})\hat{\mathbf{i}} + \hat{\mathbf{j}} +
\hat{\mathbf{k}}\Bigr]\\
- & \frac{\mu_2}{r_2^3}\Bigl[(x -
\pi_1r_{12})\hat{\mathbf{i}} + \hat{\mathbf{j}} +
\hat{\mathbf{k}}\Bigr]
\end{aligned}\etiqueta{12}
\end{align}
Ahora todo lo que tienes que hacer es igualar los coeficientes.
\begin{align}
\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x
& = -\frac{\mu_1}{r_1^3}(x + \pi_2r_{12}) - \frac{\mu_2}{r_2^3}(x
- \pi_1r_{12})\\
\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y
& = -\frac{\mu_1}{r_1^3}y - \frac{\mu_2}{r_2^3}y\\
\ddot{z}
& = -\frac{\mu_1}{r_1^3}z - \frac{\mu_2}{r_2^3}z
\end{align}
Ahora tenemos el sistema no lineal de ecuaciones diferenciales ordinarias. Podemos suponer que la trayectoria es en avión y lo hacemos dejando $z = 0$, por lo que sólo tenemos dos ecuaciones restantes:
\begin{align}
\ddot{x} - 2\varOmega\dot{y} - \varOmega^2x
& = -\frac{\mu_1}{r_1^3}(x + \pi_2r_{12}) - \frac{\mu_2}{r_2^3}(x
- \pi_1r_{12})\\
\ddot{y} + 2\varOmega\dot{x} - \varOmega^2y
& = -\frac{\mu_1}{r_1^3}y - \frac{\mu_2}{r_2^3}y
\end{align}
Ahora para lograr la requieren trayectoria podemos utilizar la verdadera anomalía de $-90^{\circ}$, un ángulo de trayectoria de vuelo de $20^{\circ}$, y una velocidad inicial de $v_0 = 10.9148$ km/s. En este punto, no tenemos otra opción sino para integrar numéricamente. Aquí es un código que va a alcanzar el nivel deseado de la trama. Sin embargo, en la realidad apolo vuelo, que tenía a mediados de correcciones de rumbo y llegaron exactamente en el otro lado de la Tierra.
Python
#!/usr/bin/env ipython
import numpy as np
import scipy
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import brentq
from scipy.optimize import fsolve
import pylab
me = 5.974 * 10 ** 24 # mass of the earth
mm = 7.348 * 10 ** 22 # mass of the moon
G = 6.67259 * 10 ** -20 # gravitational parameter
re = 6378.0 # radius of the earth in km
rm = 1737.0 # radius of the moon in km
r12 = 384400.0 # distance between the CoM of the earth and moon
M = me + mm
d = 200 # distance the spacecraft is above the Earth
pi2 = mm / M
mue = 398600.0 # gravitational parameter of earth km^3/sec^2
mum = G * mm # grav param of the moon
omega = np.sqrt(mu / r12 ** 3)
vbo = 10.9148
nu = -np.pi*0.5
gamma = 20*np.pi/180.0 # angle in radians of the flight path
vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
# velocity of the bo in the x direction
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
# velocity of the bo in the y direction
xrel = (re + d)*np.cos(nu)-pi2*r12
# spacecraft x location relative to the earth
yrel = (re + 200.0) * np.sin(nu)
u0 = [xrel, yrel, 0, vx, vy, 0]
def deriv(u, dt):
return [u[3], # dotu[0] = u[3]
u[4], # dotu[1] = u[4]
u[5], # dotu[2] = u[5]
(2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
(u[0] - pi1 * r12) /
np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
# dotu[3] = that
(-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
# dotu[4] = that
0] # dotu[5] = 0
dt = np.linspace(0.0, 535000.0, 535000.0) # secs to run the simulation
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = u.T
# 3d plot
fig = pylab.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, color = 'r')
# 2d plot if you want
#fig2 = pylab.figure(2)
#ax2 = fig2.add_subplot(111)
#ax2.plot(x, y, color = 'r')
# adding the moon
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = rm * np.outer(np.cos(phi), np.sin(theta)) + r12 - pi2 * r12
ym = rm * np.outer(np.sin(phi), np.sin(theta))
zm = rm * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
# adding the earth
xe = re * np.outer(np.cos(phi), np.sin(theta)) - pi2 * r12
ye = re * np.outer(np.sin(phi), np.sin(theta))
ze = re * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xe, ye, ze, color = '#4169E1', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
ax.set_xlim3d(-20000, 385000)
ax.set_ylim3d(-20000, 80000)
ax.set_zlim3d(-50000, 50000)
pylab.show()
![enter image description here]()
Matlab
secuencia de comandos:
days = 24*3600;
G = 6.6742e-20;
rmoon = 1737;
rearth = 6378;
r12 = 384400;
m1 = 5974e21;
m2 = 7348e19;
M = m1 + m2;
pi_1 = m1/M;
pi_2 = m2/M;
mu1 = 398600;
mu2 = 4903.02;
mu = mu1 + mu2;
W = sqrt(mu/r12^3);
x1 = -pi_2*r12;
x2 = pi_1*r12;
d0 = 200;
phi = -90;
v0 = 10.9148;
gamma = 20;
t0 = 0;
tf = 6.2*days;
r0 = rearth + d0;
x = r0*cosd(phi) + x1;
y = r0*sind(phi);
vx = v0*(sind(gamma)*cosd(phi) - cosd(gamma)*sind(phi));
vy = v0*(sind(gamma)*sind(phi) + cosd(gamma)*cosd(phi));
f0 = [x; y; vx; vy];
[t,f] = rkf45(@rates, [t0 tf], f0);
x = f(:,1);
y = f(:,2);
vx = f(:,3);
vy = f(:,4);
xf = x(end);
yf = y(end);
vxf = vx(end);
vyf = vy(end);
df = norm([xf - x2, yf - 0]) - rmoon;
vf = norm([vxf, vyf]);
plot(x,y)
funciones: 2 por separado
function [tout, yout] = rkf45(ode_function, tspan, y0, tolerance)
a = [0 1/4 3/8 12/13 1 1/2];
b = [0 0 0 0 0
1/4 0 0 0 0
3/32 9/32 0 0 0
1932/2197 -7200/2197 7296/2197 0 0
439/216 -8 3680/513 -845/4104 0
-8/27 2 -3544/2565 1859/4104 -11/40];
c4 = [25/216 0 1408/2565 2197/4104 -1/5 0];
c5 = [16/135 0 6656/12825 28561/56430 -9/50 2/55];
if nargin < 4
tol = 1.e-8;
else
tol = tolerance;
end
t0 = tspan(1);
tf = tspan(2);
t = t0;
y = y0;
tout = t;
yout = y';
h = (tf - t0)/100; % Assumed initial time step.
while t < tf
hmin = 16*eps(t);
ti = t;
yi = y;
for i = 1:6
t_inner = ti + a(i)*h;
y_inner = yi;
for j = 1:i-1
y_inner = y_inner + h*b(i,j)*f(:,j);
end
f(:,i) = feval(ode_function, t_inner, y_inner);
end
te = h*f*(c4' - c5');
te_max = max(abs(te));
ymax = max(abs(y));
te_allowed = tol*max(ymax,1.0);
delta = (te_allowed/(te_max + eps))^(1/5);
if te_max <= te_allowed
h = min(h, tf-t);
t = t + h;
y = yi + h*f*c5';
tout = [tout;t];
yout = [yout;y'];
end
h = min(delta*h, 4*h);
if h < hmin
fprintf(['\n\n Warning: Step size fell below its minimum\n'...
' allowable value (%g) at time %g.\n\n'], hmin, t)
return
end
end
segunda función:
function dfdt = rates(t,f)
r12 = 384400;
m1 = 5974e21;
m2 = 7348e19;
M = m1 + m2;
pi_1 = m1/M;
pi_2 = m2/M;
mu1 = 398600;
mu2 = 4903.02;
mu = mu1 + mu2;
W = sqrt(mu/r12^3);
x1 = -pi_2*r12;
x2 = pi_1*r12;
x = f(1);
y = f(2);
vx = f(3);
vy = f(4);
r1 = norm([x + pi_2*r12, y]);
r2 = norm([x - pi_1*r12, y]);
ax = 2*W*vy + W^2*x - mu1*(x - x1)/r1^3 - mu2*(x - x2)/r2^3;
ay = -2*W*vx + W^2*y - (mu1/r1^3 + mu2/r2^3)*y;
dfdt = [vx; vy; ax; ay];
end
![enter image description here]()
