Fixes compat issues

JuliaAstro · Jul 13, 2024 · c2f5def · c2f5def
1 parent ca91416
commit c2f5def
Show file tree

Hide file tree

Showing 7 changed files with 1,521 additions and 148 deletions.
diff --git a/lib/AstrodynamicalCalculations/test/runtests.jl b/lib/AstrodynamicalCalculations/test/runtests.jl
diff --git a/lib/AstrodynamicalModels/Project.toml b/lib/AstrodynamicalModels/Project.toml
@@ -4,7 +4,6 @@ authors = ["Joseph D Carpinelli <joey@loopy.codes>"]
 version = "3.8.0"
 
 [deps]
-AstrodynamicalCalculations = "c0cf9fb7-f496-4999-a425-c50785d1b88b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
@@ -22,7 +21,7 @@ AstrodynamicalCalculationsExt = "AstrodynamicalCalculations"
 SPICEBodiesExt = "SPICEBodies"
 
 [compat]
-AstrodynamicalCalculations = "1"
+AstrodynamicalCalculations = "0.5, 1"
 DocStringExtensions = "0.9"
 Memoize = "0.4"
 ModelingToolkit = "9.3"

diff --git a/lib/AstrodynamicalModels/src/Attitude.jl b/lib/AstrodynamicalModels/src/Attitude.jl
@@ -201,9 +201,11 @@ model = Attitude()
     eqs =
         vcat(δ.(q) .~ (1 // 2) * (A * q), δ.(ω) .~ (-inv(J) * ωx * J * ω + inv(J) * L + f))
 
+    u = vcat(q, ω)
+
     if stm
         @variables (Φ(t))[1:7, 1:7] [description = "state transition matrix estimate"]
-        A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v))
+        A = Symbolics.jacobian(map(el -> el.rhs, eqs), u)
 
         Φ = Symbolics.scalarize(Φ)
 
@@ -224,7 +226,7 @@ model = Attitude()
         return ODESystem(
             eqs,
             t,
-            vcat(q, ω, vec(Φ)),
+            vcat(u, vec(Φ)),
             vcat(vec(J), L, f);
             name = name,
             defaults = defaults,
@@ -234,7 +236,7 @@ model = Attitude()
         return ODESystem(
             eqs,
             t,
-            vcat(q, ω),
+            u,
             vcat(vec(J), L, f);
             name = name,
             defaults = defaults,

diff --git a/lib/AstrodynamicalModels/src/NBP.jl b/lib/AstrodynamicalModels/src/NBP.jl
@@ -65,6 +65,7 @@ model = NBSystem(9)
     )
 
     eqs = vcat(poseqs, veleqs)
+    u = vcat(r..., v...)
 
     if stm
         if N ≥ 3
@@ -76,8 +77,10 @@ model = NBSystem(9)
             """
         end
 
-        @variables (Φ(t))[1:6, 1:6] [description = "state transition matrix estimate"]
-        A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v))
+        @variables (Φ(t))[1:length(eqs), 1:length(eqs)] [
+            description = "state transition matrix estimate",
+        ]
+        A = Symbolics.jacobian(map(el -> el.rhs, eqs), u)
 
         Φ = Symbolics.scalarize(Φ)
 
@@ -98,7 +101,7 @@ model = NBSystem(9)
         return ODESystem(
             eqs,
             t,
-            vcat(r..., v..., Φ...),
+            vcat(u, Φ...),
             vcat(G, m...);
             name = modelname,
             defaults = defaults,
@@ -108,7 +111,7 @@ model = NBSystem(9)
         return ODESystem(
             eqs,
             t,
-            vcat(r..., v...),
+            u,
             vcat(G, m...);
             name = modelname,
             defaults = defaults,

diff --git a/lib/AstrodynamicalSolvers/src/CR3BSolvers.jl b/lib/AstrodynamicalSolvers/src/CR3BSolvers.jl
@@ -11,7 +11,7 @@ $(IMPORTS)
 """
 module CR3BSolvers
 
-export halo, lyapunov, monodromy
+export halo, lyapunov
 
 using LinearAlgebra
 using StaticArrays
@@ -435,82 +435,4 @@ function halo(
 
 end
 
-"""
-Solve for the monodromy matrix of the periodic orbit.
-"""
-function monodromy(
-    u::AbstractVector,
-    μ,
-    T;
-    algorithm = Vern9(),
-    reltol = 1e-12,
-    abstol = 1e-12,
-    save_everystep = false,
-    kwargs...,
-)
-    problem = ODEProblem(
-        CR3BFunction(stm = true),
-        MVector{42}(
-            u[begin],
-            u[begin+1],
-            u[begin+2],
-            u[begin+3],
-            u[begin+4],
-            u[begin+5],
-            1,
-            0,
-            0,
-            0,
-            0,
-            0,
-            0,
-            1,
-            0,
-            0,
-            0,
-            0,
-            0,
-            0,
-            1,
-            0,
-            0,
-            0,
-            0,
-            0,
-            0,
-            1,
-            0,
-            0,
-            0,
-            0,
-            0,
-            0,
-            1,
-            0,
-            0,
-            0,
-            0,
-            0,
-            0,
-            1,
-        ),
-        (zero(T), T),
-        (μ,),
-    )
-    solution = solve(
-        problem,
-        algorithm;
-        reltol = reltol,
-        abstol = abstol,
-        save_everystep = save_everystep,
-        kwargs...,
-    )
-
-    if solution.u[begin][begin:begin+5] ≉ solution.u[end][begin:begin+5]
-        @warn "The orbit does not appear to be periodic!"
-    end
-
-    return reshape((solution.u[end][begin+6:end]), 6, 6)
-end
-
 end # module
diff --git a/lib/AstrodynamicalSolvers/src/Propagation.jl b/lib/AstrodynamicalSolvers/src/Propagation.jl
@@ -132,4 +132,84 @@ function monodromy(
     return AstrodynamicalModels.CartesianSTM(solution.u[end][begin+N:end])
 end
 
+
+"""
+Solve for the monodromy matrix of the periodic orbit.
+"""
+function monodromy(
+    u::AbstractVector,
+    μ,
+    T,
+    f;
+    algorithm = Vern9(),
+    reltol = 1e-12,
+    abstol = 1e-12,
+    save_everystep = false,
+    kwargs...,
+)
+    problem = ODEProblem(
+        f,
+        MVector{42}(
+            u[begin],
+            u[begin+1],
+            u[begin+2],
+            u[begin+3],
+            u[begin+4],
+            u[begin+5],
+            1,
+            0,
+            0,
+            0,
+            0,
+            0,
+            0,
+            1,
+            0,
+            0,
+            0,
+            0,
+            0,
+            0,
+            1,
+            0,
+            0,
+            0,
+            0,
+            0,
+            0,
+            1,
+            0,
+            0,
+            0,
+            0,
+            0,
+            0,
+            1,
+            0,
+            0,
+            0,
+            0,
+            0,
+            0,
+            1,
+        ),
+        (zero(T), T),
+        (μ,),
+    )
+    solution = solve(
+        problem,
+        algorithm;
+        reltol = reltol,
+        abstol = abstol,
+        save_everystep = save_everystep,
+        kwargs...,
+    )
+
+    if solution.u[begin][begin:begin+5] ≉ solution.u[end][begin:begin+5]
+        @warn "The orbit does not appear to be periodic!"
+    end
+
+    return reshape((solution.u[end][begin+6:end]), 6, 6)
+end
+
 end