Skip to content

Instantly share code, notes, and snippets.

@rebcabin
Last active September 6, 2017 00:07
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rebcabin/180b54ca8d7dafc25806c6d80baada39 to your computer and use it in GitHub Desktop.
Save rebcabin/180b54ca8d7dafc25806c6d80baada39 to your computer and use it in GitHub Desktop.
Kepler Problem, Lagrangian form
ClearAll[toSI];
toSI[magn_, unit_] :=
UnitConvert[Quantity[magn, unit]] // QuantityMagnitude;
<< NumericalCalculus`
ClearAll[r, \[Theta], \[Tau], realizer, G, mE, mS, M, \[Mu],
dES, \[Theta]dot00, day];
realizer = {G -> toSI[6.67408 10^-11, "m^3/kg/s^2"],
mE -> QuantityMagnitude[
UnitConvert[
Entity["Planet", "Earth"][EntityProperty["Planet", "Mass"]]]],
mS -> QuantityMagnitude[
UnitConvert[
Entity["Star", "Sun"][EntityProperty["Star", "Mass"]]]],
dES ->
QuantityMagnitude[UnitConvert[Quantity[1, "AstronomicalUnit"]]],
\[Theta]dot00 -> toSI[360/365.25, "angular degree/day"],
day -> toSI[1, "day"],
M -> mE + mS,
\[Mu] -> (mE mS)/(mE + mS)};
Module[{coordinates, velocities, nonConservativeForces,
L, equations, numericalEquations,
initialConditions, numericalInitialConditions,
tlim = 400 day //. realizer,
soln},
coordinates = {r[\[Tau]], \[Theta][\[Tau]]};
velocities = D[coordinates, \[Tau]];
L = 1/2 \[Mu] r'[\[Tau]]^2 +
1/2 \[Mu] r[\[Tau]]^2 \[Theta]'[\[Tau]]^2 + G M \[Mu]/r[\[Tau]];
nonConservativeForces = -{0, 0}*velocities;
equations =
MapThread[{q, v, h} \[Function] D[D[L, v], \[Tau]] - D[L, q] == h,
{coordinates, velocities, nonConservativeForces}];
numericalEquations = equations //. realizer;
initialConditions = {
r[0] == dES, r'[0] == 0,
\[Theta][0] == 0, \[Theta]'[0] == \[Theta]dot00};
numericalInitialConditions = initialConditions //. realizer;
soln = NDSolve[Join[numericalEquations, numericalInitialConditions],
coordinates, {\[Tau], 0, tlim}, Method -> "BDF"];
Grid[{
{"Gravitational Force", G M \[Mu]/dES^2 //. realizer},
{"Centripetal Force", \[Mu] dES \[Theta]dot00^2 //. realizer},
{"Lagrangian", L},
{"Equations of Motion", equations},
{Plot[{r[\[Tau]] /. soln[[1]]}, {\[Tau], 0, tlim},
ImageSize -> Medium, Frame -> True,
FrameLabel -> {{"r [meter]", ""}, {"t [sec]",
"Radius: Earth-Sun Kepler Problem"}},
LabelStyle -> Directive[Red, Bold, Medium]],
Plot[{\[Theta][\[Tau]]/\[Degree] /. soln[[1]]}, {\[Tau], 0,
tlim},
ImageSize -> Medium, Frame -> True,
FrameLabel -> {{"\[Theta] [degree]", ""}, {"t [sec]",
"Angle: Earth-Sun Kepler Problem"}},
LabelStyle -> Directive[Red, Bold, Medium]]
},
{With[{\[Theta]dot =
Quiet@ND[\[Theta][\[Tau]] /. soln[[1]], \[Tau], t]},
Plot[{\[Mu] r[\[Tau]]^2 \[Theta]dot/\[Degree] /. soln[[1]] //.
realizer /. {\[Tau] -> t}}, {t, 0, tlim},
ImageSize -> Medium, Frame -> True,
FrameLabel -> {{"\[ScriptCapitalL] [joule-sec]",
""}, {"t [sec]",
"Angular Momentum: Earth-Sun Kepler Problem"}},
LabelStyle -> Directive[Red, Bold, Medium]]],
With[{\[Theta]dot =
Quiet@ND[\[Theta][\[Tau]] /. soln[[1]], \[Tau], t],
rdot = Quiet@ND[\[Theta][\[Tau]] /. soln[[1]], \[Tau], t]},
Plot[{1/2 \[Mu] rdot^2 + 1/2 \[Mu] r[\[Tau]]^2 \[Theta]dot^2 -
G M \[Mu]/r[\[Tau]] /. soln[[1]] //.
realizer /. {\[Tau] -> t}}, {t, 0, tlim},
ImageSize -> Medium, Frame -> True,
FrameLabel -> {{"\[ScriptCapitalT]+\[ScriptCapitalV] [joule]",
""}, {"t [sec]", "Total Energy: Earth-Sun Kepler Problem"}},
LabelStyle -> Directive[Red, Bold, Medium]]]}
}]]
ClearAll[r, \[Theta], \[Tau], realizer, G, mE, mS, M, \[Mu],
dES, \[Theta]dot00, day];
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment