# Earth to Mars Interplanetary Travel Escape Velocity [closed]

I'm trying to create a script that simulates interplanetary travel from Earth to Mars using a Hohmann transfer, but I'm having a few difficulties getting the correct escape velocity. I'm using an x-y Heliocentric coordinate system, so I have to break down the escape velocity into its x- and y-components and then add it to the initial conditions of the spacecraft that is orbiting Earth. I've tried using some standard interplanetary escape velocity equations, but so far have come up with the spacecraft going in the right direction, but doing it with too much speed. Here's the code that solves a system of ODEs (I apologize for the horribe looking code, I'm not too sure how to add special Mathematica symbols like epsilons and sigmas to a forum post):

(*Gravitational Constant*)
G = 6.672*10^-11;
Subscript[t, max] = 100000000;
(*Masses of Sun, Earth, Mars and satellite*)
m[0] = AstronomicalData["Sun", "Mass"];
m[1] = AstronomicalData["Earth", "Mass"];
m[2] = AstronomicalData["Mars", "Mass"];
m[3] = 100;
(*Heliocentric positions of Earth, Mars and satellite*)
p[1] = AstronomicalData["Earth", {"Position", {2013, 1, 1, 0, 0, 1}}]
p[2] = AstronomicalData["Mars", {"Position", {2013, 18, 1, 0, 0, 1}}]
p[3] = AstronomicalData[
"Earth", {"Position", {2013, 1, 1, 0, 0, 1}}] + {0, -r[1] -
200000, 0};
(*Planetary radii of Sun, Earth and Mars*)
r[0] = AstronomicalData["Sun", "Radius"];
r[1] = AstronomicalData["Earth", "Radius"];
r[2] = AstronomicalData["Mars", "Radius"];
(*Heliocentric radii components of Earth and Mars*)
rex = Norm[p[1][[1]]];
rey = Norm[p[1][[2]]];
rmx = Norm[p[2][[1]]];
rmy = Norm[p[2][[2]]];
(*Heliocentric velocities of Earth, Mars and spacecraft calculated \
via central difference approximation f'(x)=(f(x+h)-f(x-h))/2h*)
v[1] = ((AstronomicalData[
"Earth", {"Position", {2013, 1, 1, 0, 0, 1.001}}] -
AstronomicalData[
"Earth", {"Position", {2013, 1, 1, 0, 0, 0.999}}]))/(2 (0.001));
v[2] = ((AstronomicalData[
"Mars", {"Position", {2013, 18, 1, 0, 0, 1.001}}] -
AstronomicalData[
"Mars", {"Position", {2013, 18, 1, 0, 0, 0.999}}]))/(2 (0.001));
(*Hyperbolic escape velocities for spacecraft travelling to Mars*)
v1x = Sqrt[(G m[0])/rex] (Sqrt[(2 rmx)/(rex + rmx)] - 1)
v1y = Sqrt[(G m[0])/rey] (Sqrt[(2 rmy)/(rey + rmy)] - 1)
(*Spacecraft initial velocity*)
v[3] = v[1] + {-v1x, -v1y, 0}
(*Softening factor*)
ϵ = 0.001;
(*Numerical solution to 4-body problem, modelling the gravitation \
interation between the Sun, Earth, Mars and an artificial satelite*)
soln = NDSolve[{
x[0]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(x[0]$$[t] - $$x[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[0]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[0]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[0]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
y[0]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(y[0]$$[t] - $$y[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[0]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[0]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[0]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
z[0]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(z[0]$$[t] - $$z[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[0]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[0]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[0]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
x[1]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(x[1]$$[t] - $$x[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[1]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[1]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[1]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
y[1]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(y[1]$$[t] - $$y[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[1]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[1]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[1]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
z[1]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(z[1]$$[t] - $$z[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[1]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[1]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[1]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
x[2]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(x[2]$$[t] - $$x[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[2]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[2]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[2]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
y[2]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(y[2]$$[t] - $$y[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[2]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[2]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[2]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
z[2]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(z[2]$$[t] - $$z[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[2]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[2]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[2]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
x[3]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(x[3]$$[t] - $$x[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[3]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[3]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[3]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
y[3]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(y[3]$$[t] - $$y[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[3]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[3]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[3]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),
z[3]''[t] == -\!$$\*UnderoverscriptBox[\(∑$$, $$i = 0$$, $$3$$]$$\((G\ m[i] \((\(z[3]$$[t] - $$z[i]$$[t])\))\)/
\*SuperscriptBox[$$( \*SuperscriptBox[\((\(x[3]$$[t] - $$x[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(y[3]$$[t] - $$y[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$(\(z[3]$$[t] - $$z[i]$$[t])\), $$2$$] +
\*SuperscriptBox[$$ϵ$$, $$2$$])\), $$3/2$$]\)\),

x[0][0] == 0, y[0][0] == 0, z[0][0] == 0, x[1][0] == p[1][[1]],
y[1][0] == p[1][[2]], z[1][0] == p[1][[3]], x[2][0] == p[2][[1]],
y[2][0] == p[2][[2]], z[2][0] == p[2][[3]], x[3][0] == p[3][[1]],
y[3][0] == p[3][[2]], z[3][0] == p[3][[3]], x[0]'[0] == 0,
y[0]'[0] == 0, z[0]'[0] == 0, x[1]'[0] == v[1][[1]],
y[1]'[0] == v[1][[2]], z[1]'[0] == v[1][[3]],
x[2]'[0] == v[2][[1]], y[2]'[0] == v[2][[2]],
z[2]'[0] == v[2][[3]], x[3]'[0] == v[3][[1]],
y[3]'[0] == v[3][[2]], z[3]'[0] == v[3][[3]]}, {x[0][t], y[0][t],
z[0][t], x[1][t], y[1][t], z[1][t], x[2][t], y[2][t], z[2][t],
x[3][t], y[3][t], z[3][t]}, {t, 0, Subscript[t, max]},
Method -> "StiffnessSwitching", MaxSteps -> 10000000]


And here's the plot code:

orbits = ParametricPlot3D[{{x[0][t], y[0][t], z[0][t]}, {x[1][t],
y[1][t], z[1][t]}, {x[2][t], y[2][t], z[2][t]}, {x[3][t],
y[3][t], z[3][t]}} /. soln, {t, 0, Subscript[t, max]},
AxesLabel -> {x, y, z}, ImageSize -> Large,
PlotRange -> {{-0.3*10^12, 0.3*10^12}, {-0.3*10^12, 0.3*10^12}, {-0.3*10^12, 0.3*10^12}}]


As you can see, the output shows the spacecraft shooting off into outer space and never returning, whereas what I would like to see is the spacecraft just getting to the orbital path of Mars (I'll worry about timing and phase angles a little later).

Does anyone know what I'm doing wrong? Any help would be appreciated, thanks guys.

-

## closed as off-topic by MarcoB, blochwave, m_goldberg, Yves Klett, halirutanSep 30 '15 at 6:30

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – MarcoB, blochwave, halirutan
If this question can be reworded to fit the rules in the help center, please edit the question.

I think your problem is more in the mathematics your using than a coding problem. Could be wrong, of course. However, I suggest you look into the concept of a Hohmann transfer orbit. – m_goldberg Sep 13 '13 at 14:14
The delta-v equation that I'm using is the delta-v1 equation that is given in the wiki article (delta-v1 = sqrt(mu/r1)(sqrt(2r2/(r1+r2))-1), but I think my problem is that I'm not putting them into their correct vector forms of the equations and don't really know how to do that, unless I just keep them in their scalar/magnitude form and then multiply by cos(theta) to get the x-component and sin(theta) to get the y-component. – user7388 Sep 13 '13 at 15:41
I'm voting to close this question as off-topic because the issue it raises is not a Mathematica issue but a mathematics or physics problem. That it is formulated in terms of Mathematica is not sufficient to make it an appropriate question for Mathematica.SE. – m_goldberg Sep 29 '15 at 23:12