Plotting a parametric solution from NDSolve

Im trying to animate multiple points around multiple parametric plots which depicts orbital motion of body around a planet. Initial conditons are;

μ = 3.986004418*10^14
a = {7.92597218162462*^6, 7.359004757830201*^6,
6.970300551929753*^6}
r = {7.388961739897817*^6, 7.352270990873303*^6,
6.989358622500435*^6}
v = {7589.490674834727, 7366.424498576473,7541.464522475139}
Ecc = {0.014, 0.025, 0.03}(*eccentricity*)
RAAN = {150, 160, 170}
Inc = {20, 80, 320}
ArgPer = {220, 110, 330}


The following parametric equations are

Table[{x[t] =
a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Cos[t],{y[t] =
a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Sin[t]}}, {i, 1, 3}]


which are used in NDSolve

orbsoln =
Table[NDSolve[{(x''[
t] + (μ*(x[t][[i]]))/((x[t][[i]])^2 + (y[t][[i]])^2)^(3/2) ==
0, (y''[
t] + (μ* (y[t][[i]]))/((x[t][[i]])^2 + (y[t][[i]])^2)^(3/2) ==
0, x[0] == r[[j]], y[0] == 0,
Derivative[1][x][0] == 0, Derivative[1][y][0] == v[[j]]}, {x ,
y}, {t, 0, 2 Pi}], {j, 1, 3}]


The results are then plotted using Animate and rotated into plane using an euler rotation matrix

EulerRotationMatrix[{i_, Ω_, ω_}] := RotationMatrix[
i, {0, 0, 1}].RotationMatrix[Ω, {0, 1,
0}].RotationMatrix[ω, {0, 0, 1}];

plot = Table[
ParametricPlot3D[
Evaluate[(EulerRotationMatrix[{RAAN[[m]] Degree, Inc[[m]] Degree,
ArgPer[[m]] Degree}].{x[t], y[t], 0}) /. First[orbsoln]], {t,
0, 2 Pi}, PlotRange -> 0.8*10^7], {m, 1, 3}]

Animate[Show[plot,
Graphics3D[{PointSize[0.02], Red,
Point[Table[
Evaluate[(EulerRotationMatrix[{RAAN[[o]] Degree,
Inc[[o]] Degree, ArgPer[[o]] Degree}].{x[t], y[t], 0}) /.
First[orbsoln]], {o, 1, 3}]]}]], {t, 0, 2 Pi}]


I'm aiming to show 3 points following their parametric path that has been defined by NDsolve, but all i get is this plot;

and the error message

Coordinate {{{6.978455209925671*^6, 1.4717698904511775*^6, 1.9318209890130716*^6}}, {{-1.9573444509619826]^6, -6.669739440694639^6, 2.5060436058655437*^6}}, {{-4.195765089771511*^6, 4.47488342280073*^6, 4.119204986980711*^6}}} should be a triple of numbers, or a Scaled form.

Can anyone point out where I have gone wrong?

• orbsoln doesn't work properly Mar 29, 2019 at 23:15

2 Answers

Your expression for assigning values to x[t] and y[t] look suspect.

Table[
{x[t] = a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Cos[t],
{y[t] = a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Sin[t]}},
{i, 1, 3}]


generates

but most of that is thrown away;x[t], y[t] are just

Column[{x[t], y[t]}]


I find it hard to believe this what you want.

Update

After I spent quite a bit more time studying your code I concluded that main error is in formulating orbsoln. However, I going to ignore it, because I think it might not be necessary for what you want to do.

μ = 3.986004418*10^14;
a = {7.92597218162462*^6, 7.359004757830201*^6, 6.970300551929753*^6};
r = {7.388961739897817*^6, 7.352270990873303*^6, 6.989358622500435*^6};
v = {7589.490674834727, 7366.424498576473, 7541.464522475139};
Ecc = {0.014, 0.025, 0.03};
RAAN = {150, 160, 170};
Inc = {20, 80, 320};
ArgPer = {220, 110, 330};

EulerRotationMatrix[{i_, Ω_, ω_}] :=
RotationMatrix[i, {0, 0, 1}] .
RotationMatrix[Ω, {0, 1, 0}] .
RotationMatrix[ω, {0, 0, 1}]


All you really need is the table of orbit expressions written this way.

Block[{t, orbits2D},
orbits2D =
Transpose[
{Table[
a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Cos[t], {i, 1, 3}],
Table[
a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]]*Cos[t])*Sin[t], {i, 1, 3}],
{0, 0, 0}}];
orbits3D =
{orb[1][t_], orb[2][t_], orb[3][t_]} =
Table[
EulerRotationMatrix[
{RAAN[[m]] Degree, Inc[[m]] Degree, ArgPer[[m]] Degree}].orbits2D[[m]],
{m, 1, 3}]];


Now we can make the background plot showing the orbit traces with

plot = ParametricPlot3D[Evaluate @ orbits3D, {t, 0, 2 Pi}];


and make a demonstration showing particles following those obits with Manipulate. I leave it to you to convert the Manipulate to an Animate expression if that's what you prefer. It's an easy exercise.

Manipulate[
Show[
plot,
Graphics3D[{Red, AbsolutePointSize[5], Point[Evaluate[orb[#][t] & /@ Range[3]]]}]],
{t, 0, 2 Pi}]


Note: As the problem is reformulated here, the independent variable t does not represent time. So this is a purely kinematic solution. The orbital parameters μ, r and v are not used. I think you were trying to get the true dynamics of the orbits with orbsoln, but I find that code so confusing that I don't see how to fix it.

• Thanks for your reply, im trying to use the individual column values of x[t] and y[t] in NDsolve, i thought this was working as parametric plots were generated for the respective parametric paths. Could you explain why although the paths are generated, the animated points are not? Mar 30, 2019 at 9:01

I thought it might be useful as a reference to demonstrate a side-by-side comparison of two ways to visualize Keplerian orbits from their orbital elements, either by using the exact polar equations, or by using ParametricNDSolveValue[]:

μ = QuantityMagnitude[PlanetData["Earth", "GravitationalConstantMassProduct"]];

(* semimajor axes *)
A = {7.92597218162462*^6, 7.359004757830201*^6, 6.970300551929753*^6};
(* eccentricities *)
ε = {0.014, 0.025, 0.03};

(* argument of periapsis *)
ω = {220, 110, 330} °;
(* inclination *)
ι = {20, 80, 320} °;
(* longitude of the ascending node *)
Ω = {150, 160, 170} °;

kepler[a_, e_, θ_] := a (1 - e^2)/(1 + e Cos[θ])

orbit = ParametricNDSolveValue[{r''[t] == -μ r[t]/Norm[r[t]]^3,
r[0] == {a (1 - e), 0},
r'[0] == {0, Sqrt[μ (1 + e)/(a (1 - e))]}},
r, {t, 0, 2 π a Sqrt[a/μ]}, {a, e}];

GraphicsRow[{ParametricPlot3D[MapThread[
Function[{a, e, Ω, ι, ω},
Append[kepler[a, e, θ] AngleVector[θ], 0] .
EulerMatrix[{Ω, ι, ω}, {3, 1, 3}]],
{A, ε, Ω, ι, ω}] // Evaluate, {θ, 0, 2 π},
PlotLabel -> "Exact"],
ParametricPlot3D[MapThread[Function[{a, e, Ω, ι, ω},
Append[IdentityMatrix[2], {0, 0}] .
orbit[a, e][a Sqrt[a/μ] t] .
EulerMatrix[{Ω, ι, ω}, {3, 1, 3}]],
{A, ε, Ω, ι, ω}] // Evaluate, {t, 0, 2 π},
PlotLabel -> "NDSolve"]},
PlotLabel -> "Keplerian Orbits"]