Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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*)
    \[Epsilon] = 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[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    y[0]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    z[0]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    x[1]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    y[1]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    z[1]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    x[2]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    y[2]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    z[2]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    x[3]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    y[3]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(2\)])\), \(3/2\)]\)\),
    z[3]''[t] == -\!\(
    \*UnderoverscriptBox[\(\[Sum]\), \(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[\(\[Epsilon]\), \(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).

enter image description here

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

share|improve this 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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.