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've written an interplanetary trajectory solver/plotter that plots the path taken by a spacecraft on an Earth-Mars mission, but have run into a little trouble when the spacecraft actually reaches Mars. I'm hoping to be able to get the spacecraft into an orbit around Mars upon arrival, but before I can do that I need to find its closest approach in order to apply the correct delta-v at the right time (using the wonders of NDSolve's WhenEvent function). The trajectory plot relies on the output of an interpolating function from NDSolve, and it is this interpolating function that I am trying to work with. So far I know when the closest approach occurs to within a day, and it is between days 254 and 255 of the trip. I have put the Heliocentric x-y positions of Mars and the spacecraft (between days 254 and 255) into separate tables using the following:

MarsPosition = Table[{x[2][t], y[2][t]} /. Soln, {t, 254*86400, 255*86400, 86400/100}] 
SpaceCraftPosition = Table[{x[3][t], y[3][t]} /. Soln, {t, 254*86400, 255*86400, 86400/100}]

And have then calculated their relative x-y positions using the following:

dxy = Sqrt[(MarsPosition - SpaceCraftPosition)^2]

But I am now trying to calculate the minimum approach radius (i.e. combining the x-y positions into a magnitude) from the above bits of code. My futile attempt of trying

dr = Min[Sqrt[dxy[1]^2 + dxy[2]^2]]

and other similar attempts have unfortunately ended in failure. Some sample output looks like this:

  dxy = {{{2.00946*10^8, 9.69241*10^7}}, {{2.03432*10^8, 
       9.83663*10^7}}, {{2.05917*10^8, 9.98081*10^7}},...}

Where, for example, the first entry shows the x and y relative position as {2.00946*10^8, 9.69241*10^7} and I'm trying to combine those two components into a radius. Combining the x-y position pairs into a radius using Sqrt[x^2 + y^2] for the first 3 position pairs given above I would get

 dr = {{2.231*10^8},{2.25965*10^8},{2.2883*10^8},...}

But this was done manually for each entry and is not practical for hundreds of entries. From this list of radius positions, I can then use Min[] to find the closes approach point for the spacecraft relative to Mars.

As can be seen below, the spacecraft gets very close to Mars (where the inner semi-circle is Earth's orbit, the outer semi-circle is Mars' orbit and the line joining the two is the spacecraft's trajectory):enter image description here

And if propagated further we can see that it gets close enough for Mars' gravity to affect its trajectory:enter image description here

So it looks pretty close, but exactly how close it what I'm hoping to find out. Any help would be greatly appreciated.

share|improve this question
You need to give us Soln first or all your flybys are belong to us :D PS: Look up Norm. –  Yves Klett Mar 14 '14 at 12:33
I was hoping someone would know a general way to work with tables where each entry gives an x-y position pair. I can give the code for Soln if needed, but it is a ton of code what will take some time to edit into a "stackexchange friendly" format. –  user7388 Mar 14 '14 at 12:36
yeah, but some data to work with would be great, otherwise (how) are we supposed to make that up? Did you already look up Norm or Nearest? –  Yves Klett Mar 14 '14 at 12:38
I would think you would be able to use Soln and FindMinimum. You could probably find the closest point with WhenEvent, too, if Soln is produced by NDSolve. Hard to tell. –  Michael E2 Mar 14 '14 at 13:11
Very true, my apologies, I've edited the original post and given some example output for the first 3 x-y position pairs from the dxy table as well as what I hope to get as a table for dr. Doing that now I realise that I should have just give a general example in the first place instead of writing out what I'm using it for. –  user7388 Mar 14 '14 at 13:13

1 Answer 1

up vote 1 down vote accepted

I guess I'll answer my own question in case someone else needs to use it for future reference. The following code was used to find the minimum approach radius:

    MarsPosition = Table[{x[2][t], y[2][t]} /. Soln, {t, 253*86400, 256*86400, 60}] ;
    SpaceCraftPosition = Table[{x[3][t], y[3][t]} /. Soln, {t, 253*86400, 256*86400, 60}] ;
    (*xy positions for spacecraft relative to Mars*)
    dxy = Sqrt[(MarsPosition - SpaceCraftPosition)^2];
    (*Table used to hold approach radii*)
    dr = Table[Norm[dxy[[i]]], {i, 1, Length[dxy]}]
    (*Minimum approach radius*)
    mindr = Min[dr]
    (*Index position at which mindr occurred*)
    mindrindex = Position[dr, mindr]
    (*Time at which mindr occurred, t = (253 days + 60n seconds), where n is the index of mindr*)
    mindrtime = (253*86300) + (60*mindrindex)

Thanks for trying to help Yves and Michael, I appreciate it.

share|improve this answer
You are welcome! Only without incomplete code to this question will only be partially helpful. You could add in some test data to make it that much better (and get some more rep). –  Yves Klett Mar 14 '14 at 14:20
Yup, I agree it was very vague. I have edited it to include test data for the dxy table. –  user7388 Mar 14 '14 at 14:25

Your Answer


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

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