Last active
September 6, 2017 00:07
-
-
Save rebcabin/180b54ca8d7dafc25806c6d80baada39 to your computer and use it in GitHub Desktop.
Kepler Problem, Lagrangian form
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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