[WIP] Upgrade modeling dependencies (i.e. ModelingToolkit & Symbolics)

Project.toml

@@ -53,7 +53,7 @@ GeneralizedGenerated = "0.3.2"
 KernelDensity = "0.6.2"
 LineSearches = "7.1.1"
 Memoize = "0.4.4"
-ModelingToolkit = "5.6.4"
+ModelingToolkit = "6"
 NLsolve = "4.5.1"
 OrdinaryDiffEq = "5.49.1"
 PhysicalConstants = "0.2.1"
@@ -67,7 +67,7 @@ SciMLBase = "1.8"
 SimpleTraits = "0.9.3"
 StaticArrays = "1"
 StatsPlots = "0.14.17"
-SymbolicUtils = "0.8.4, 0.9, 0.10, 0.11"
-Symbolics = "0.1.2"
+SymbolicUtils = "0.13"
+Symbolics = "3"
 Unitful = "1.5.0"
 julia = "1.6"

src/dynamics/dynamical_models.jl

@@ -35,7 +35,8 @@ end
     eqs, dvs = ode_order_lowering(eqs, independent_variable(_ode), ModelingToolkit.states(_ode)) 
     eqs = [eqs[end - num_orig_eqs + 1:end]..., eqs[1:num_orig_eqs]...]
     dvs = [dvs[end - num_orig_eqs + 1:end]..., dvs[1:num_orig_eqs]...]
-    ODESystem(simplify.(eqs), independent_variable(_ode), dvs, parameters(_ode))
+
+    ODESystem(simplify.(eqs), independent_variable(_ode), dvs, parameters(_ode); name=nameof(_ode))
 end
 
 

src/dynamics/models/bc4bp.jl

@@ -38,7 +38,7 @@ function ModelingToolkit.ODESystem(::Type{BC4BP_ODESystem})
     eqs    = @. D2(p) ~ [+2D(y) + x, -2D(x) + y, 0] + forces
 
     # Simplify and return the system
-    return ODESystem(simplify.(eqs), t, [x, y, z], [μ, μ2, a3, α0, a3, ω3])
+    return ODESystem(simplify.(expand_derivatives.(eqs)), t, [x, y, z], [μ, μ2, a3, α0, a3, ω3]; name=:BC4BP)
 end
 
 # Build the equations at pre-compile time

src/dynamics/models/cr3bp.jl

@@ -14,7 +14,7 @@ function ModelingToolkit.ODESystem(::Type{CR3BP_ODESystem})
     eqs_er3bp = ODESystem(ER3BP_ODESystem)
     (μ, e) = parameters(eqs_er3bp)
     eqs = [eq.lhs ~ simplify(substitute(eq.rhs, e => 0)) for eq in equations(eqs_er3bp)]
-    return ODESystem(eqs, independent_variable(eqs_er3bp), states(eqs_er3bp), [μ])
+    return ODESystem(eqs, independent_variable(eqs_er3bp), states(eqs_er3bp), [μ]; name=:CR3BP)
 end
 
 # Build the equations at pre-compile time

src/dynamics/models/er3bp.jl

@@ -9,7 +9,7 @@ struct ER3BP_ODESystem{S,F} <: Abstract_AstrodynamicalODESystem
     ode_f      :: F
 end
 
-function ModelingToolkit.ODESystem(::Type{ER3BP_ODESystem})
+function ModelingToolkit.ODESystem(::Type{T}) where {T<:ER3BP_ODESystem}
     @parameters  μ  # Mass fraction
     @parameters  e  # Eccentricity
     @parameters  f  # True anomaly
@@ -20,14 +20,15 @@ function ModelingToolkit.ODESystem(::Type{ER3BP_ODESystem})
     # [DeiTos2017, Eqs. 14]
     # NOTE: need to expand the RHS derivatives before calling order-lowering
     ω = elliptical_potential(μ, (x, y, z), f, e)
-    return ODESystem([
+    return ODESystem(expand_derivatives.([
             D2(x) ~ +2D(y) + Dx(ω),
             D2(y) ~ -2D(x) + Dy(ω),
             D2(z) ~        + Dz(ω)
-        ], 
+        ]),
         f,
         [x, y, z],
-        [μ, e]
+        [μ, e];
+        name=nameof(T)
     )
 end
 

src/dynamics/models/nbp.jl

@@ -69,7 +69,7 @@ end
     eqs = @. D2(pos) ~ sum(accelerations)
 
     # Build the 2nd-order ODE props
-    return ODESystem(eqs, t, pos, [])
+    return ODESystem(expand_derivatives.(eqs), t, pos, []; name=:EphemerisNBP)
 end
 
 #---------------------#