|
1 | 1 | # Examples |
2 | 2 |
|
3 | 3 | _Usage examples._ |
4 | | - |
5 | | -## Periodic Orbits |
6 | | - |
7 | | -This package contains differential correctors, and helpful wrapper functions, for |
8 | | -finding periodic orbits within Circular Restricted Three Body Problem dynamics. |
9 | | - |
10 | | -### `Plots.jl` |
11 | | - |
12 | | -```{julia} |
13 | | -#| echo: true |
14 | | -using AstrodynamicalSolvers |
15 | | -using AstrodynamicalModels |
16 | | -using OrdinaryDiffEq |
17 | | -using Plots |
18 | | -
|
19 | | -μ = 0.012150584395829193 |
20 | | -
|
21 | | -planar = let |
22 | | - ic = halo(μ, 1) # lyapunov (planar) orbit |
23 | | - u = Vector(CartesianState(ic)) |
24 | | - problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,)) |
25 | | - solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14) |
26 | | - Plots.plot(solution, idxs=(:x,:y,:z), title = "Lyapunov Orbit", label=:none, size=(1600,900), dpi=400, aspect_ratio=1) |
27 | | -end |
28 | | -
|
29 | | -extraplanar = let |
30 | | - ic = halo(μ, 2; amplitude=0.01) # halo (non-planar) orbit |
31 | | - u = Vector(CartesianState(ic)) |
32 | | - problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,)) |
33 | | - solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14) |
34 | | - Plots.plot(solution, idxs=(:x,:y,:z), title = "Halo Orbit", label=:none, size=(1600,900), dpi=400, aspect_ratio=1) |
35 | | -end |
36 | | -
|
37 | | -Plots.plot(planar, extraplanar, layout=(1,2)) |
38 | | -``` |
39 | | - |
40 | | -### `Makie.jl` |
41 | | - |
42 | | -```{julia} |
43 | | -#| echo: true |
44 | | -using AstrodynamicalSolvers |
45 | | -using AstrodynamicalModels |
46 | | -using OrdinaryDiffEq |
47 | | -using CairoMakie |
48 | | -
|
49 | | -μ = 0.012150584395829193 |
50 | | -
|
51 | | -sol_planar = let |
52 | | - ic = halo(μ, 1) # lyapunov (planar) orbit |
53 | | - u = Vector(CartesianState(ic)) |
54 | | - problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,)) |
55 | | - solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14) |
56 | | -end |
57 | | -
|
58 | | -sol_extraplanar = let |
59 | | - ic = halo(μ, 2; amplitude=0.01) # halo (non-planar) orbit |
60 | | - u = Vector(CartesianState(ic)) |
61 | | - problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,)) |
62 | | - solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14) |
63 | | -end |
64 | | -
|
65 | | -fig = CairoMakie.Figure(size=(800, 400); fontsize=11) |
66 | | -
|
67 | | -ax_kwargs_common = (; aspect=:equal, azimuth=-π/3) |
68 | | -
|
69 | | -ax_left = CairoMakie.Axis3(fig[1, 1]; |
70 | | - title = "Lyapunov Orbit", |
71 | | - limits = (0.78, 0.90, -0.09, 0.09, -0.02, 1.04), |
72 | | - ax_kwargs_common..., |
73 | | -) |
74 | | -ax_right = CairoMakie.Axis3(fig[1, 2]; |
75 | | - title = "Halo Orbit", |
76 | | - limits = (1.05, 1.26, -0.1, 0.1, -0.02, 0.01), |
77 | | - protrusions = (30, 100, 0, 0), |
78 | | - ax_kwargs_common..., |
79 | | -) |
80 | | -
|
81 | | -CairoMakie.plot!(ax_left, sol_planar; idxs=(:x, :y, :z)) |
82 | | -CairoMakie.plot!(ax_right, sol_extraplanar; idxs=(:x, :y, :z)) |
83 | | -
|
84 | | -fig |
85 | | -``` |
86 | | - |
87 | | -## Manifold Computations |
88 | | - |
89 | | -Manifold computations, provided by `AstrodynamicalCalculations.jl`, can perturb |
90 | | -halo orbits onto their unstable or stable manifolds. |
91 | | - |
92 | | -### `Plots.jl` |
93 | | - |
94 | | -```{julia} |
95 | | -#| echo: true |
96 | | -using AstrodynamicalSolvers |
97 | | -using AstrodynamicalCalculations |
98 | | -using AstrodynamicalModels |
99 | | -using OrdinaryDiffEq |
100 | | -using LinearAlgebra |
101 | | -using Plots |
102 | | -
|
103 | | -μ = 0.012150584395829193 |
104 | | -
|
105 | | -unstable = let |
106 | | - ic = halo(μ, 1; amplitude=0.005) |
107 | | -
|
108 | | - u = CartesianState(ic) |
109 | | - Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true)) |
110 | | -
|
111 | | - ics = let |
112 | | - problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,)) |
113 | | - solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10)) |
114 | | -
|
115 | | - solution.u |
116 | | - end |
117 | | -
|
118 | | - perturbations = [ |
119 | | - diverge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=-1e-7) |
120 | | - for ic in ics |
121 | | - ] |
122 | | -
|
123 | | - problem = EnsembleProblem( |
124 | | - ODEProblem(CR3BFunction(), u, (0.0, 2 * ic.Δt), (μ,)), |
125 | | - prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]), |
126 | | - ) |
127 | | -
|
128 | | - solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14) |
129 | | -end |
130 | | -
|
131 | | -stable = let |
132 | | - ic = halo(μ, 2; amplitude=0.005) |
133 | | -
|
134 | | - u = CartesianState(ic) |
135 | | - Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true)) |
136 | | -
|
137 | | - ics = let |
138 | | - problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,)) |
139 | | - solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10)) |
140 | | -
|
141 | | - solution.u |
142 | | - end |
143 | | -
|
144 | | - perturbations = [ |
145 | | - converge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=1e-7) |
146 | | - for ic in ics |
147 | | - ] |
148 | | -
|
149 | | - problem = EnsembleProblem( |
150 | | - ODEProblem(CR3BFunction(), u, (0.0, -2.1 * ic.Δt), (μ,)), |
151 | | - prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]), |
152 | | - ) |
153 | | -
|
154 | | - solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14) |
155 | | -end |
156 | | -
|
157 | | -figure = Plots.plot(; |
158 | | - aspect_ratio = 1.0, |
159 | | - background = :transparent, |
160 | | - grid = true, |
161 | | - title = "Unstable and Stable Invariant Manifolds", |
162 | | -) |
163 | | -
|
164 | | -Plots.plot!(figure, unstable, idxs=(:x, :y), aspect_ratio=1, label=:none, palette=:blues) |
165 | | -Plots.plot!(figure, stable, idxs=(:x, :y), aspect_ratio=1, label=:none, palette=:blues) |
166 | | -Plots.scatter!(figure, [1-μ], [0], label="Moon", xlabel="X (Earth-Moon Distance)", ylabel="Y (Earth-Moon Distance)", marker=:x, color=:black, markersize=10,) |
167 | | -
|
168 | | -figure # hide |
169 | | -``` |
170 | | - |
171 | | -## `Makie.jl` |
172 | | - |
173 | | -```{julia} |
174 | | -#| echo: true |
175 | | -using AstrodynamicalSolvers |
176 | | -using AstrodynamicalCalculations |
177 | | -using AstrodynamicalModels |
178 | | -using OrdinaryDiffEq |
179 | | -using LinearAlgebra |
180 | | -using CairoMakie |
181 | | -
|
182 | | -μ = 0.012150584395829193 |
183 | | -
|
184 | | -unstable = let |
185 | | - ic = halo(μ, 1; amplitude=0.005) |
186 | | -
|
187 | | - u = CartesianState(ic) |
188 | | - Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true)) |
189 | | -
|
190 | | - ics = let |
191 | | - problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,)) |
192 | | - solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10)) |
193 | | -
|
194 | | - solution.u |
195 | | - end |
196 | | -
|
197 | | - perturbations = [ |
198 | | - diverge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=-1e-7) |
199 | | - for ic in ics |
200 | | - ] |
201 | | -
|
202 | | - problem = EnsembleProblem( |
203 | | - ODEProblem(CR3BFunction(), u, (0.0, 2 * ic.Δt), (μ,)), |
204 | | - prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]), |
205 | | - ) |
206 | | -
|
207 | | - solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14) |
208 | | -end |
209 | | -
|
210 | | -stable = let |
211 | | - ic = halo(μ, 2; amplitude=0.005) |
212 | | -
|
213 | | - u = CartesianState(ic) |
214 | | - Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true)) |
215 | | -
|
216 | | - ics = let |
217 | | - problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,)) |
218 | | - solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10)) |
219 | | -
|
220 | | - solution.u |
221 | | - end |
222 | | -
|
223 | | - perturbations = [ |
224 | | - converge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=1e-7) |
225 | | - for ic in ics |
226 | | - ] |
227 | | -
|
228 | | - problem = EnsembleProblem( |
229 | | - ODEProblem(CR3BFunction(), u, (0.0, -2.1 * ic.Δt), (μ,)), |
230 | | - prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]), |
231 | | - ) |
232 | | -
|
233 | | - solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14) |
234 | | -end |
235 | | -
|
236 | | -fig = CairoMakie.Figure(size=(800, 400), fontsize=20) |
237 | | -
|
238 | | -ax = CairoMakie.Axis(fig[1, 1]; |
239 | | - xreversed = true, |
240 | | - xticks = LinearTicks(5), |
241 | | - yticks = LinearTicks(5), |
242 | | - aspect = DataAspect(), |
243 | | - xlabel = "X (Earth-Moon Distance)", |
244 | | - ylabel = "Y (Earth-Moon Distance)", |
245 | | - title = "Unstable and Stable Invariant Manifolds", |
246 | | - titlesize = 24, |
247 | | -) |
248 | | -
|
249 | | -idxs = (:x, :y) |
250 | | -
|
251 | | -# TODO: replace this manual workaround when |
252 | | -# https://github.com/SciML/SciMLBase.jl/issues/697#issuecomment-2135801331 |
253 | | -# is addressed |
254 | | -for (traj, color) in zip(unstable, resample_cmap(:blues, length(unstable))) |
255 | | - CairoMakie.plot!(ax, traj; idxs, color) |
256 | | -end |
257 | | -
|
258 | | -for (traj, color) in zip(stable, resample_cmap(:blues, length(stable))) |
259 | | - CairoMakie.plot!(ax, traj; idxs, color) |
260 | | -end |
261 | | -
|
262 | | -CairoMakie.scatter!(ax, [1-μ], [0]; marker='⨯', color=:black, markersize=50, label="Moon") |
263 | | -
|
264 | | -fig |
265 | | -``` |
0 commit comments