# Finding Minimium Periapse Distance Between Spacecraft and Planet

I have the following simple model that simulates a spacecraft orbiting Mars. What I'm hoping to do is be able to find the periapse radius of the spacecraft as I did below, except instead of using tables I'll use the interpolating functions that are given as output from NDSolve (would using interpolating functions to find the periapse radius be faster than using tables?):

G = 6.672*10^-11;
m = 6.4185*10^23;
m = 100;
p = {0, 0};
p = {1000000, 1000000};
v = {0, 2500};
v = {0, 0};
tmax = 1000;

soln = NDSolve[{
x''[t] == -(G m (x[t] - x[t]))/((x[t] - x[t])^2 + (y[t] - y[t])^2)^(3/2),
y''[t] == -(G m (y[t] - y[t]))/((x[t] - x[t])^2 + (y[t] - y[t])^2)^(3/2),
x''[t] == -(G m (x[t] - x[t]))/((x[t] - x[t])^2 + (y[t] - y[t])^2)^(3/2),
y''[t] == -(G m (y[t] - y[t]))/((x[t] - x[t])^2 + (y[t] - y[t])^2)^(3/2),
x == p[], y == p[], x == p[],
y == p[], x' == v[],
y' == v[], x' == v[],
y' == v[]}, {x[t], y[t], x[t], y[t],
x[t], y[t], x[t], y[t]}, {t, 0, tmax},
Method -> "StiffnessSwitching", AccuracyGoal -> 18,
PrecisionGoal -> 18, MaxSteps -> 10000000];

ParametricPlot[{{{x[t], y[t], {x[t], y[t]}}}} /. soln, {t,
0, tmax}, Prolog -> {Red, Disk[{p[], p[]}, 50000]},
AxesLabel -> {x, y}, ImageSize -> Large]

Animate[ParametricPlot[{{x[t], y[t]}, {x[t], y[t]}} /.
soln /. t -> a, {t, Max[0, a - 1000], a},
Prolog -> {Red, Disk[{p[], p[]}, 50000]},
AxesLabel -> {x, y}, Axes -> False, ImageSize -> Large], {a, 0,
tmax}, AnimationRate -> 10]

dt = 1;
MarsPosition = Table[{x[t], y[t]} /. soln, {t, 0, tmax, dt}] ;
SpaceCraftPosition = Table[{x[t], y[t]} /. soln, {t, 0, tmax, dt}];
dxy = Sqrt[(MarsPosition - SpaceCraftPosition)^2];
dr = Table[Norm[dxy[[i]]], {i, 1, Length[dxy]}];
mindr = Min[dr] (*Find closest approach of spacecraft to Mars*)
mindrindex = Flatten[Position[dr, mindr]][] (*Finds index position of mindr in dr*)


I've already had a look at a few topics that used FindMinValue, FindMinimum and NMinimize (Finding the Minimum value of an interpolating function, Find minimum value of function which calls NDsolve, Find maximum value of interpolation function - obviously wrong result), but could not get them to work as I would always get syntax errors that I could not solve. Would using these functions be the correct way of achieving what I'd like to do, and if so, could anyone help me figure out how to get past the syntax errors? Any help would be brilliant. • I think your system would benefit greatly from a judicious rescaling of your variables. Jun 13, 2015 at 13:13

G = 6672*10^-14;
m = 64185*10^19;
m = 100;
p = {0, 0};
p = {10^6, 10^6};
v = {0, 2500};
v = {0, 0};
tmax = 1000;

soln = NDSolve[{
x''[t] == -(G m (x[t] - x[t]))/
Norm[{x[t], y[t]} - {x[t], y[t]}]^3,
y''[t] == -(G m (y[t] - y[t]))/
Norm[{x[t], y[t]} - {x[t], y[t]}]^3,
x''[t] == -(G m (x[t] - x[t]))/
Norm[{x[t], y[t]} - {x[t], y[t]}]^3,
y''[t] == -(G m (y[t] - y[t]))/
Norm[{x[t], y[t]} - {x[t], y[t]}]^3,
x == p[], y == p[],
x == p[], y == p[],
x' == v[], y' == v[],
x' == v[], y' == v[]},
{x[t], y[t], x[t], y[t],
x[t], y[t], x[t], y[t]},
{t, 0, tmax},
Method -> "StiffnessSwitching",
AccuracyGoal -> 18,
PrecisionGoal -> 18,
MaxSteps -> 10000000];

mindr = NMinimize[Join[
Norm[{x[t], y[t]} - {x[t], y[t]}] /. soln,
{0 <= t <= tmax}], t]


{76703.6, {t -> 909.167}}

ParametricPlot[
{{x[t], y[t]}, {x[t], y[t]}} /. soln,
{t, 0, tmax},
Prolog -> {
Red, Disk[{p[], p[]}, 50000],
Magenta, AbsolutePointSize,
Point[{x[t], y[t]} /. soln /. mindr[[-1]]]},
AxesLabel -> {x, y},
ImageSize -> Large] • Brilliant, thank you Bob, worked like a charm. I'd also like to be able to get the x and y velocity components at periapse, but how would I isolate 909.167 from {t -> 909.167}? I tried mindr[] which produced {t -> 909.167}, so x'[t] /. soln /. t ->mindr[] unfortunately did not work. Jun 13, 2015 at 18:12
• Nevermind, using Last[Last[mindr[]]] did the trick. Jun 13, 2015 at 20:34

What you want to do is called "event location" and is realized with NDSolve using WhenEvent. In principle you give it a predicate that is true when the spaceship is at periapsis and NDSolve uses a root finding method to figure out exactly when this happens.

G = 6.672*10^-11;
m = 6.4185*10^23;
m = 100;
p = {0, 0};
p = {1000000, 1000000};
v = {0, 2500};
v = {0, 0};
tmax = 650;

system = {
{x''[t], y''[t]} == -G m {x[t] - x[t], y[t] - y[t]}/Norm[{x[t] - x[t], y[t] - y[t]}]^3,
{x''[t], y''[t]} == -G m {x[t] - x[t], y[t] - y[t]}/Norm[{x[t] - x[t], y[t] - y[t]}]^3,
{x, y} == p,
{x, y} == p,
{x', y'} == v,
{x', y'} == v,
WhenEvent[
{x[t] - x[t], y[t] - y[t]}.{x'[t] - x'[t], y'[t] - y'[t]} == 0,
Sow[t]
]
};

initial = {
x[t],
y[t],
x[t],
y[t],
x[t],
y[t],
x[t],
y[t]
};

{soln, periapsis} = Reap@NDSolve[
system,
initial,
{t, 0, tmax},
Method -> "StiffnessSwitching",
AccuracyGoal -> 18,
PrecisionGoal -> 18,
MaxSteps -> 10000000
];

pts = Table[Flatten[{x[t], y[t]} /. soln], {t, Flatten@periapsis}]

ParametricPlot[
{{x[t], y[t], {x[t], y[t]}}} /. soln,
{t, 0, tmax},
Prolog -> {
Red,
Disk[p, 50000],
Green,
PointSize[Large],
Point[pts]
},
AxesLabel -> {x, y}
] The predicate that I'm using is when the dot product between the velocity vector and the position vector changes sign. This produces two candidates; the periapsis and the apoapsis, as can be seen in the figure.

Note: I don't know which root finding algorithm NDSolve uses but conceptually it will be similar to the secant method. Therefore if you request a high precision is might tell you that it couldn't find the exact place at the precision that you requested because it required too many iterations. But it will tell you how well it can approximate it.

• The built-in root-finder is in fact Brent's method, which is effectively a hybrid of bisection, the secant method, and inverse quadratic interpolation. Jun 13, 2015 at 15:01