I've been making a Mathematica script whereby a spacecraft is transferred from Earth to Mars using a Hohmann interplanetary transfer by means of a patched conic approximation.

As far as I can tell all of the equations and numbers should be correct, yet when running the code the spacecraft doesn't make it to Mars. Without getting into any detail of the code (that will have to be for the Mathematica SX forums), I wanted to ask how accurate the patched conic approximation really is. If all of the numbers are correct should I be able to run the code and the spacecraft will hit Mars, or is it supposed to be used to find a rough trajectory to another planet? I understand the inherent inaccuracies with the approximation since it considers everything in terms of individual 2-body problems, and so thought I might be a little too hopeful when I thought it would get me to Mars.

enter image description here

EDIT: Okay, I think my initial post used an incorrect angle for applying the delta-v impulse which I have since fixed (I think!). I made a second program that looked at just the Earth and spacecraft upon departure to confirm whether my angle was correct and this is what I got (the spacecraft doesn't actually pass through the green circle (Earth). The green circle is just the position of Earth at the start of the simulation):

enter image description here

As you can see, the spacecraft's escape velocity is parallel to the Earth's velocity which is a good start. I got the correct angle by finding the angle of the Earth's velocity vector using phi = arctan(v_y/v_x), the "impulse angle" (not too sure what to call it) by using beta = arccos(1/e) and then using theta = phi + beta + pi when finding the x- and y-components for delta-v = v_p - v_c (where v_x and v_y are the x- and y-components of the Earth's velocity respectively, e is the escape hyperbola's eccentricity and v_p and v_c are the spacecraft's periapse speed and circular orbit speed about Earth respectively).

Unfortunately, when running the full program again in the Sun's frame of reference I found the following:

1) Roughly half-period of spacecraft orbit about Sun enter image description here

2) Full period of spacecraft orbit about Sun enter image description here

P.S. The equations I'm using are from Orbital Mechanics for Engineering Students by Howard Curtis.

UPDATE: Solution to Earth-Mars interplanetary transfer using Lambert solver (some much needed coding help was given my the Mathematica SX forums).

enter image description here

  • 1
    $\begingroup$ What does your plot signify? In what manner is the craft not making it to Mars? $\endgroup$ – Colin McFaul Nov 8 '13 at 23:15
  • $\begingroup$ The plot signifies a top down view of the orbits of Earth (inner circle) and Mars (outer circle) about the Sun (located at (0,0)), as well as the path the spacecraft makes from Earth to Mars, which was made by numerically integrating Newton's Law of Gravitation for 4 bodies. I was hoping to see the spacecraft either get pulled into an orbit about Mars, or have its orbital path highly perturbed by coming into Mars' sphere of influence, but as can be seen it unfortunately passed harmlessly past Mars without even getting close (relative to astronomical distances it missed by a long shot!). $\endgroup$ – user7388 Nov 10 '13 at 0:08
  • $\begingroup$ What are the times involved? When the spacecraft makes its closest approach to Mars's orbit, is Mars also at that point in its orbit? If so, what is the distance of closet approach? I don't know this area, but these might help someone else figure what's going on. In addition, you might try flagging this for migration to Space Exploration or Astronomy, as those communities might be able to give you a better answer than this one. $\endgroup$ – Colin McFaul Nov 10 '13 at 0:52
  • $\begingroup$ I set the date of departure as the 30th of March 2001 which I found in a table online (it was said to be within the window of opportunity). I made an animation of the above plot to see how close the spacecraft got to Mars and it looks like the craft arrived too late (the angle made between their radius vectors is about 10 degrees once the spacecraft reaches aphelion) so maybe the departure date was wrong or my equations were wrong. Can anyone confirm if the above date that I used would have been within a correct launch window? $\endgroup$ – user7388 Nov 11 '13 at 0:05
  • 1
    $\begingroup$ This is the website I used to find the launch windows, for those who are interested: clowder.net/hop/railroad/EMa.htm $\endgroup$ – user7388 Nov 11 '13 at 15:29

It's pretty hard to tell what's going on with just a plot, especially since you show an entire spacecraft orbit instead of just the half that you intended to be the transfer. In any case, one thing that I'm seeing that doesn't look right is that the spacecraft trajectory goes inside Earth's orbit. That shouldn't happen (at least not to the degree shown). That indicates that your departure velocity from Earth ($V_\infty$) has the wrong direction.

Are you integrating from a low Earth orbit? It takes some care to time that correctly so you get the right $V_\infty$ direction.

While the patched conic is just an approximation, if done correctly it should get you much closer to the orbit of Mars than what you show. You should not be able to see a gap at the resolution of that plot.

When you say in the comments "I was hoping to see the spacecraft either get pulled into an orbit about Mars ...", that's not possible. Without a propulsive maneuver at Mars, the spacecraft will fly by, as indicated in the rest of your sentence "... or have its orbital path highly perturbed by coming into Mars' sphere of influence". That will happen if you get the timing right for when you depart Earth.


Ok, that's an improvement.

Now you need to account for the eccentricity of the orbit of Mars. The spacecraft should have an aphelion that corresponds to the distance of Mars from the Sun at the planned time of the encounter. It will help to also take into account the eccentricity of the Earth's orbit, in particular the distance of the Earth from the Sun at the time of departure.

Side note:

Beware the perils of the numerical derivative. Take care to not use too small of a $\Delta T$. Here is an example in Mathematica when trying to determine the velocity vector of Earth using Earth position data. $\Delta T$ is on the x-axis in seconds.

f[dt_?NumericQ] := AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1 + dt}}] - AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1}}]

LogLinearPlot[f[dt][[1]]/dt, {dt, 0.001, 100}]

numerical derivative noise plot

| improve this answer | |
  • $\begingroup$ When you say I must take into account the eccentricity of Mars and Earth, I'm guessing you mean I must change the equation for v_infinity which is given by v_infinity = sqrt(mu_sun/R_earth(sqrt(2R_mars/(R_mars+R_earth))-1)? Would I set R_earth to Earth's heliocentric radius at spacecraft departure, and R_mars to Mars' heliocentric radius at spacecraft arrival? Currently they are both set to their heliocentric radii at departure. P.S. I apologise for the horrible looking equation, it's tough to make them look nice in a little box like this. $\endgroup$ – InquisitiveInquirer Nov 17 '13 at 12:07
  • $\begingroup$ I can't quite make out your equation, but yes, use Earth's heliocentric radius at departure, and Mars' heliocentric radius at arrival. $\endgroup$ – Mark Adler Nov 17 '13 at 17:20
  • $\begingroup$ By the way, you can use $\TeX$ equations in comments, e.g. $\sqrt{\mu_S\over r_1}\left(\sqrt{2r_2\over r_1+r_2}-1\right)$. I think that's the equation you meant, but it looks like you're missing a parenthesis. $\endgroup$ – Mark Adler Nov 17 '13 at 17:31
  • 1
    $\begingroup$ Perfect example. 0.002 seconds is way too small of a $\Delta T$ for the data you are using. Yes, you want a small $\Delta T$, but it is possible to pick too small of a $\Delta T$ that samples noise in the lower significant digits of your data. Run this in Mathematica to see what I'm talking about: f[dt_] := AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1 + dt}}] - AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1}}] and LogLinearPlot[f[dt][[1]]/dt, {dt, 0.001, 100}]. From that plot, it looks like around 10 seconds is a good $\Delta T$. $\endgroup$ – Mark Adler Nov 21 '13 at 17:57
  • 2
    $\begingroup$ For that purpose, I recommend the SPICE toolkit, along with a recent ephemeris file such as DE430. You'll also need a recent leap seconds kernel. $\endgroup$ – Mark Adler Feb 8 '14 at 22:22

If departure and destination points are 180 degrees apart (as in a Hohmann transfer), Lambert iterations will give a polar heliocentric transfer orbit. That is an orbit with inclination 90º. Obviously such a transfer orbit would take huge plane change expense at departure and destination.

This can be mitigated by using a broken plane transfer orbit. See Deboning The Porkchop Plot.

If you stay in the ecliptic plane, you can miss Mars by a good margin. The only places where Mars comes close to the ecliptic plane is around the ascending and descending nodes.

| improve this answer | |
  • $\begingroup$ The other solution, which is the one actually used, is to simply not force the departure and destination to be 180° apart. You pick your departure and destination dates to be around the minima on either side of the ridge on the porkchop plot, called a Type 1 and Type 2 trajectory, for less than or more than 180° respectively. No broken plane maneuver is needed. $\endgroup$ – Mark Adler Mar 15 '14 at 22:21
  • $\begingroup$ There are two types of direction change to consider here. Plane change as well as flight path angle change. For coplanar orbits you would want the transfer orbit to be tangent to the departure and destination orbit. Since Mars has a noticeably elliptical orbit, the tangent transfer orbit can be more or less than 180 degrees, but usually in that neighborhood. $\endgroup$ – HopDavid Mar 15 '14 at 23:27
  • $\begingroup$ To avoid plane change expense you want launch and or departure to be near a descending or ascending node. Sometimes this can result in decidedly un-Hohmann like orbits where the transfer orbit only covers 100º from departure to destination. This would result in a major change for heliocentric flight path angle. $\endgroup$ – HopDavid Mar 15 '14 at 23:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.