Dear Ch'Prof,
There is an example animation at the end of this. Just evaluate the code.
I've been working on developing Mathematica content for teaching STEM subjects at MIT. I often get asked to develop content from an astrodynamics class.
I am going to cut and paste some code that I debugged for a colleague showing a Mars-Earth transfer orbit.
Unfortunately, this is only debugged and not commented. I think you may see some of the benefits of teaching via coding like this.
Imagine teaching your students to code this up, and not just watch it.
Again, this is not a good and compact example of WL programming style, but it works.
r11 = 1; r21 = 1.5237`; \[Mu]sAU = 0.000295912`; v1 = \
Sqrt[\[Mu]sAU/r11]
v2 = Sqrt[\[Mu]sAU/r21]
valpha = Sqrt[\[Mu]sAU (2/r21 - 2/(r11 + r21))]
deltav = valpha - v2
solmefunc[theta0_, alpha_] :=
NDSolve[{(r1^\[Prime]\[Prime])[t] - r1[t]*th1'[t]^2 == -\[Mu]sAU/
r1[t]^2,
r1[t]*(th1^\[Prime]\[Prime])[t] + 2*r1'[t]*th1'[t] ==
0, (r2^\[Prime]\[Prime])[t] - r2[t]*th2'[t]^2 == -\[Mu]sAU /
r2[t]^2,
r2[t]*(th2^\[Prime]\[Prime])[t] + 2*r2'[t]*th2'[t] ==
0, (r3^\[Prime]\[Prime])[t] - r3[t]*th3'[t]^2 == -\[Mu]sAU /
r3[t]^2, r3[t]*(th3^\[Prime]\[Prime])[t] + 2*r3'[t]*th3'[t] == 0,
r1[0] == r11, th1[0] == 0, Derivative[1][r1][0] == 0,
Derivative[1][th1][0] == Sqrt[\[Mu]sAU/r11], r2[0] == r21,
th2[0] == 0, Derivative[1][r2][0] == 0, Derivative[1][th2][0] ==
Sqrt[(\[Mu]sAU/r21)]/r21, r3[0] == r21, th3[0] == 0,
Derivative[1][r3][0] == alpha*deltav Sin[theta0],
Derivative[1][th3][0] == (v2 + alpha*deltav Cos[theta0])/r21}, {r1,
th1, r2, th2, r3, th3}, {t, -100, 1000}]
Show[
ParametricPlot[{rsol1[t - 99]*Cos[thsol1[t - 99]],
rsol1[t - 99]*Sin[thsol1[t - 99]]}, {t, 0, 500},
PlotStyle -> {Blue, Thick}, AspectRatio -> 1,
PlotRange -> {{-1.6, 1.6}, {-1.6, 1.6}}],
ParametricPlot[{rsol2[t]*Cos[thsol2[t]],
rsol2[t]*Sin[thsol2[t]]}, {t, 0, 500}, PlotStyle -> {Red, Thick},
PlotRange -> {{-1.6, 1.6}, {-1.6, 1.6}}],
ParametricPlot[{rsol3[t]*Cos[thsol3[t]],
rsol3[t]*Sin[thsol3[t]]}, {t, 0, 500}, PlotStyle -> {Green, Thick},
PlotRange -> {{-1.6, 1.6}, {-1.6, 1.6}}]]
Manipulate[{{rsol1, thsol1, rsol2, thsol2, rsol3, thsol3}} = {r1, th1,
r2, th2, r3, th3} /. solmefunc[theta0, alpha]; Show[
ParametricPlot[{rsol1[t]*Cos[thsol1[t]],
rsol1[t]*Sin[thsol1[t]]}, {t, 0, 800}, PlotStyle -> {Blue, Thick},
AspectRatio -> 1, PlotRange -> {{-2.2, 1.8}, {-2.0, 2.0}}],
ParametricPlot[{rsol2[t]*Cos[thsol2[t]],
rsol2[t]*Sin[thsol2[t]]}, {t, 0, 800}, PlotStyle -> {Red, Thick},
PlotRange -> {{-2.2, 1.8}, {-2.0, 2.0}}],
ParametricPlot[{rsol3[t]*Cos[thsol3[t]],
rsol3[t]*Sin[thsol3[t]]}, {t, 0, 800},
PlotStyle -> {Green, Thick},
PlotRange -> {{-2.2, 1.8}, {-2.0, 2.0}}]], {{theta0, -Pi/3}, -Pi/3,
Pi/3}, {{alpha, 1}, 1, 6}]
TimeGap = 83;
anim = Animate[
ParametricPlot[{{rsol1[t - TimeGap]*Cos[thsol1[t - TimeGap]],
rsol1[t - TimeGap]*Sin[thsol1[t - TimeGap]]}, {rsol2[t]*
Cos[thsol2[t]],
rsol2[t]*Sin[thsol2[t]]}, {rsol3[t]*Cos[thsol3[t]],
rsol3[t]*Sin[thsol3[t]]}}, {t, 0, T},
PlotStyle -> { {Blue, Thick}, {Red, Thick}, {Green, Thick}},
AspectRatio -> 1, PlotRange -> {{-1.8, 1.6}, {-1.6, 1.6}}], {T, 0,
259 }]