File:2 body problem on a sphere.gif
Summary
| Description |
English: The 2-body problem on a plane produces nice elliptical orbits. But the 2-body problem on a sphere (meaning that the distance is computed as the great-circle distance) can easily get quite chaotic. |
| Date | |
| Source | https://twitter.com/j_bertolotti/status/1263482519606919168 |
| Author | Jacopo Bertolotti |
| Permission (Reusing this file) |
https://twitter.com/j_bertolotti/status/1030470604418428929 |
Mathematica 12.0 code
ctop[p_, r_] :=
CoordinateTransform["Cartesian" -> "Spherical",
p + {Sqrt[r^2 - 1], 0, 0}]
m1 = 1; m2 = 1; G = 40;
dist[\[Theta]1_, \[Theta]2_, \[Phi]1_, \[Phi]2_, r_] :=
r ArcTan[(\[Sqrt]((Cos[\[Pi]/2 - \[Theta]2] Sin[
RealAbs[\[Phi]1 - \[Phi]2]])^2 + (Cos[\[Pi]/
2 - \[Theta]1] Sin[\[Pi]/2 - \[Theta]2] -
Sin[\[Pi]/2 - \[Theta]1] Cos[\[Pi]/2 - \[Theta]2] Cos[
RealAbs[\[Phi]1 - \[Phi]2]])^2))/(Sin[\[Pi]/
2 - \[Theta]1] Sin[\[Pi]/2 - \[Theta]2] +
Cos[\[Pi]/2 - \[Theta]1] Cos[\[Pi]/2 - \[Theta]2] Cos[
RealAbs[\[Phi]1 - \[Phi]2]])]
rl = 1000;
point1l = ctop[{0, 0, 4}, rl] // N
point2l = ctop[{0, 0.1, -3}, rl] // N
stoc1l = CoordinateTransform[
"Spherical" -> "Cartesian", {rl, \[Theta]1[t], \[Phi]1[t]}];
stoc2l = CoordinateTransform[
"Spherical" -> "Cartesian", {rl, \[Theta]2[t], \[Phi]2[t]}];
L[t_] := m1/2 Total@(D[stoc1l, t]^2) +
m2/2 Total@(D[stoc2l, t]^2) + ((m1 + m2) G)/
dist[\[Theta]1[t], \[Theta]2[t], \[Phi]1[t], \[Phi]2[t], rl]
eq\[Theta]1 = D[D[L[t], \[Theta]1'[t] ], t] - D[L[t], \[Theta]1[t]] ;
eq\[Phi]1 = D[D[L[t], \[Phi]1'[t] ], t] - D[L[t], \[Phi]1[t]] ;
eq\[Theta]2 = D[D[L[t], \[Theta]2'[t] ], t] - D[L[t], \[Theta]2[t]] ;
eq\[Phi]2 = D[D[L[t], \[Phi]2'[t] ], t] - D[L[t], \[Phi]2[t]] ;
soll = NDSolve[{eq\[Theta]1 == 0, eq\[Phi]1 == 0, eq\[Theta]2 == 0,
eq\[Phi]2 == 0, \[Theta]1[0] == point1l[[2]], \[Phi]1[0] ==
point1l[[3]], \[Theta]2[0] == point2l[[2]], \[Phi]2[0] ==
point2l[[3]], \[Theta]1'[0] == 0., \[Phi]1'[0] == -3/
rl, \[Theta]2'[0] == 0., \[Phi]2'[0] == 3/rl}, {\[Theta]1[
t], \[Phi]1[t], \[Theta]2[t], \[Phi]2[t]}, {t, 0, 200}]
rs = 10;
point1s = ctop[{0, 0, 4}, rs] // N
point2s = ctop[{0, 0.1, -3}, rs] // N
stoc1s = CoordinateTransform[
"Spherical" -> "Cartesian", {rs, \[Theta]1[t], \[Phi]1[t]}];
stoc2s = CoordinateTransform[
"Spherical" -> "Cartesian", {rs, \[Theta]2[t], \[Phi]2[t]}];
L[t_] := m1/2 Total@(D[stoc1s, t]^2) +
m2/2 Total@(D[stoc2s, t]^2) + ((m1 + m2) G)/
dist[\[Theta]1[t], \[Theta]2[t], \[Phi]1[t], \[Phi]2[t], rs]
eq\[Theta]1 = D[D[L[t], \[Theta]1'[t] ], t] - D[L[t], \[Theta]1[t]] ;
eq\[Phi]1 = D[D[L[t], \[Phi]1'[t] ], t] - D[L[t], \[Phi]1[t]] ;
eq\[Theta]2 = D[D[L[t], \[Theta]2'[t] ], t] - D[L[t], \[Theta]2[t]] ;
eq\[Phi]2 = D[D[L[t], \[Phi]2'[t] ], t] - D[L[t], \[Phi]2[t]] ;
sols = NDSolve[{eq\[Theta]1 == 0, eq\[Phi]1 == 0, eq\[Theta]2 == 0,
eq\[Phi]2 == 0, \[Theta]1[0] == point1s[[2]], \[Phi]1[0] ==
point1s[[3]], \[Theta]2[0] == point2s[[2]], \[Phi]2[0] ==
point2s[[3]], \[Theta]1'[0] == 0., \[Phi]1'[0] == -6/
rs, \[Theta]2'[0] == 0., \[Phi]2'[0] == 6/rs}, {\[Theta]1[
t], \[Phi]1[t], \[Theta]2[t], \[Phi]2[t]}, {t, 0, 500}]
p0 = Table[
Show[
ContourPlot3D[(x + Sqrt[rl^2 - 1])^2 + (y)^2 + (z)^2 ==
rl^2, {x, -50, 50}, {y, -50, 50}, {z, -50, 50},
ContourStyle -> None, Mesh -> 5, MeshStyle -> Directive[ White],
Background -> Black, Boxed -> False, Axes -> False,
PlotPoints -> 80]
,
Graphics3D[{
Orange,
Sphere[Evaluate[{rl Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rl Cos[\[Theta]1[t]]} /. soll][[1]] /. {t -> \[Tau]}, 1],
Cyan,
Sphere[
Evaluate[{rl Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rl Cos[\[Theta]2[t]]} /. soll][[1]] /. {t -> \[Tau]}, 1]
}]
,
ParametricPlot3D[
Evaluate[{rl Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[rl^2 - 1],
rl Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rl Cos[\[Theta]1[t]]} /. soll][[1]], {t, 0.01, \[Tau]},
PlotStyle -> {Opacity[0.5], Orange}]
,
ParametricPlot3D[
Evaluate[{rl Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[rl^2 - 1],
rl Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rl Cos[\[Theta]2[t]]} /. soll][[1]], {t, 0.001, \[Tau]},
PlotStyle -> {Opacity[0.5], Cyan}]
, ViewVector -> {{40, 0, 0}, {0, 0, 0}}, ViewAngle -> 50*Degree
]
, {\[Tau], 0, 66.4, 1}];
sinstep[t_] := Sin[\[Pi]/2 t]^2;
p1 = Table[
Show[
ContourPlot3D[(x + Sqrt[rl^2 - 1])^2 + (y)^2 + (z)^2 ==
rl^2, {x, -50, 50}, {y, -50, 50}, {z, -50, 50},
ContourStyle -> None, Mesh -> 5, MeshStyle -> Directive[ White],
Background -> Black, Boxed -> False, Axes -> False,
PlotPoints -> 80]
,
Graphics3D[{
Orange,
Sphere[Evaluate[{rl Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rl Cos[\[Theta]1[t]]} /. soll][[1]] /. {t -> 0}, 1],
Cyan,
Sphere[
Evaluate[{rl Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rl Cos[\[Theta]2[t]]} /. soll][[1]] /. {t -> 0}, 1]
}]
,
ParametricPlot3D[
Evaluate[{rl Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[rl^2 - 1],
rl Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rl Cos[\[Theta]1[t]]} /. soll][[1]], {t, 0, 33.2},
PlotStyle -> {Opacity[0.5 - 0.5 sinstep[\[Tau]] ], Orange}]
,
ParametricPlot3D[
Evaluate[{rl Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[rl^2 - 1],
rl Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rl Cos[\[Theta]2[t]]} /. soll][[1]], {t, 0, 33.2},
PlotStyle -> {Opacity[0.5 - 0.5 sinstep[\[Tau]] ], Cyan}]
, ViewVector -> {sinstep[\[Tau]] ({20, -20, 20} - {40, 0,
0}) + {40, 0, 0}, {0, 0, 0}}, ViewAngle -> 50*Degree
]
, {\[Tau], 0, 1, 0.1}];
p2 = Table[
Show[
ContourPlot3D[(x +
Sqrt[(sinstep[\[Tau]] (rs - rl) + rl)^2 -
1])^2 + (y)^2 + (z)^2 == (sinstep[\[Tau]] (rs - rl) +
rl)^2, {x, -50, 50}, {y, -50, 50}, {z, -50, 50},
ContourStyle -> None, Mesh -> 5, MeshStyle -> Directive[ White],
Background -> Black, Boxed -> False, Axes -> False,
PlotPoints -> 80]
,
Graphics3D[{
Orange,
Sphere[Evaluate[{rl Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rl Cos[\[Theta]1[t]]} /. soll][[1]] /. {t -> 0}, 1],
Cyan,
Sphere[
Evaluate[{rl Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[
rl^2 - 1], rl Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rl Cos[\[Theta]2[t]]} /. soll][[1]] /. {t -> 0}, 1]
}]
, ViewVector -> {{20, -20, 20}, {0, 0, 0}},
ViewAngle -> 50*Degree
]
, {\[Tau], 0, 1, 0.05}];
p3 = Table[
Show[
ContourPlot3D[(x + Sqrt[rs^2 - 1])^2 + (y)^2 + (z)^2 ==
rs^2, {x, -50, 50}, {y, -50, 50}, {z, -50, 50},
ContourStyle -> None, Mesh -> 5, MeshStyle -> Directive[ White],
Background -> Black, Boxed -> False, Axes -> False,
PlotPoints -> 30]
,
Graphics3D[{
Orange,
Sphere[Evaluate[{rs Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[
rs^2 - 1], rs Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rs Cos[\[Theta]1[t]]} /. sols][[1]] /. {t -> \[Tau]}, 1],
Cyan,
Sphere[
Evaluate[{rs Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[
rs^2 - 1], rs Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rs Cos[\[Theta]2[t]]} /. sols][[1]] /. {t -> \[Tau]}, 1]
}]
,
ParametricPlot3D[
Evaluate[{rs Cos[\[Phi]1[t]] Sin[\[Theta]1[t]] - Sqrt[rs^2 - 1],
rs Sin[\[Theta]1[t]] Sin[\[Phi]1[t]],
rs Cos[\[Theta]1[t]]} /. sols][[1]], {t, 0.01, \[Tau]},
PlotStyle -> {Opacity[0.5], Orange}]
,
ParametricPlot3D[
Evaluate[{rs Cos[\[Phi]2[t]] Sin[\[Theta]2[t]] - Sqrt[rs^2 - 1],
rs Sin[\[Theta]2[t]] Sin[\[Phi]2[t]],
rs Cos[\[Theta]2[t]]} /. sols][[1]], {t, 0.01, \[Tau]},
PlotStyle -> {Opacity[0.5], Cyan}]
, ViewVector -> {{20, -20, 20}, {0, 0, 0}},
ViewAngle -> 50*Degree
]
, {\[Tau], 0, 100, 1}];
ListAnimate[Join[p0, p1, p2, p3]]
Licensing
I, the copyright holder of this work, hereby publish it under the following license:
  |
This file is made available under the Creative Commons CC0 1.0 Universal Public Domain Dedication. |
| The person who associated a work with this deed has dedicated the work to the public domain by waiving all of their rights to the work worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission.
|
