diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index 7bba809e..c8fcedfa 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -1,25 +1,26 @@ // For format details, see https://aka.ms/devcontainer.json. For config options, see the -// README at: https://github.com/devcontainers/templates/tree/main/src/universal +// README at: https://github.com/JuliaLang/devcontainer-templates/tree/main/src/julia { - "name": "Default Linux Universal", + "name": "Julia", // Or use a Dockerfile or Docker Compose file. More info: https://containers.dev/guide/dockerfile - "image": "mcr.microsoft.com/devcontainers/universal:2-linux", + "image": "mcr.microsoft.com/devcontainers/base:ubuntu", + // Features to add to the dev container. More info: https://containers.dev/features. "features": { - "ghcr.io/julialang/devcontainer-features/julia:1": {} + // A Feature to install Julia via juliaup. More info: https://github.com/JuliaLang/devcontainer-features/tree/main/src/julia. + "ghcr.io/julialang/devcontainer-features/julia:1": {}, + "ghcr.io/devcontainers/features/python:1": {}, + "ghcr.io/devcontainers/features/github-cli": {}, + "ghcr.io/devcontainers-contrib/features/starship": {} + }, + "customizations": { + "vscode": { + "extensions": [ + "GitHub.codespaces", + "github.vscode-github-actions", + "julialang.language-julia", + "ms-toolsai.jupyter", + "valentjn.vscode-ltex" + ] + } } - - // Features to add to the dev container. More info: https://containers.dev/features. - // "features": {}, - - // Use 'forwardPorts' to make a list of ports inside the container available locally. - // "forwardPorts": [], - - // Use 'postCreateCommand' to run commands after the container is created. - // "postCreateCommand": "uname -a", - - // Configure tool-specific properties. - // "customizations": {}, - - // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root. - // "remoteUser": "root" -} +} \ No newline at end of file diff --git a/Project.toml b/Project.toml index b3716b1b..a15c4ecd 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ julia = "1.7" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [targets] -test = ["Test"] +test = ["Test", "Pkg"] diff --git a/lib/AstrodynamicalModels/Project.toml b/lib/AstrodynamicalModels/Project.toml index 70b305df..7a989b21 100644 --- a/lib/AstrodynamicalModels/Project.toml +++ b/lib/AstrodynamicalModels/Project.toml @@ -25,7 +25,7 @@ SPICEBodiesExt = "SPICEBodies" AstrodynamicalCalculations = "1" DocStringExtensions = "0.9" Memoize = "0.4" -ModelingToolkit = "<9.15" +ModelingToolkit = "9.3" SciMLBase = "2" StaticArrays = "1" Symbolics = "5" diff --git a/lib/AstrodynamicalModels/src/Attitude.jl b/lib/AstrodynamicalModels/src/Attitude.jl index 08174aa2..61105e34 100644 --- a/lib/AstrodynamicalModels/src/Attitude.jl +++ b/lib/AstrodynamicalModels/src/Attitude.jl @@ -14,8 +14,27 @@ Base.@kwdef mutable struct AttitudeState{F} <: AstrodynamicalState{F,7} AttitudeState{F}(::UndefInitializer) where {F} = new{F}() AttitudeState(::UndefInitializer) = AttitudeState{Float64}(undef) - AttitudeState{F}(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) where {F} = new{F}(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) - AttitudeState(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) = new{promote_type(typeof(q₁), typeof(q₂), typeof(q₃), typeof(q₄), typeof(ω₁), typeof(ω₂), typeof(ω₃))}(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) + AttitudeState{F}(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) where {F} = + new{F}(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) + AttitudeState(q₁, q₂, q₃, q₄, ω₁, ω₂, ω₃) = new{ + promote_type( + typeof(q₁), + typeof(q₂), + typeof(q₃), + typeof(q₄), + typeof(ω₁), + typeof(ω₂), + typeof(ω₃), + ), + }( + q₁, + q₂, + q₃, + q₄, + ω₁, + ω₂, + ω₃, + ) end """ @@ -38,8 +57,75 @@ Base.@kwdef struct AttitudeParameters{F} <: AstrodynamicalParameters{F,15} f₂::F = 0.0 f₃::F = 0.0 - AttitudeParameters{F}(J₁₁, J₂₁, J₃₁, J₁₂, J₂₂, J₃₂, J₁₃, J₂₃, J₃₃, L₁, L₂, L₃, f₁, f₂, f₃) where {F} = new{F}(J₁₁, J₂₁, J₃₁, J₁₂, J₂₂, J₃₂, J₁₃, J₂₃, J₃₃, L₁, L₂, L₃, f₁, f₂, f₃) - AttitudeParameters(J₁₁, J₂₁, J₃₁, J₁₂, J₂₂, J₃₂, J₁₃, J₂₃, J₃₃, L₁, L₂, L₃, f₁, f₂, f₃) = new{promote_type(typeof(J₁₁), typeof(J₂₁), typeof(J₃₁), typeof(J₁₂), typeof(J₂₂), typeof(J₃₂), typeof(J₁₃), typeof(J₂₃), typeof(J₃₃), typeof(L₁), typeof(L₂), typeof(L₃), typeof(f₁), typeof(f₂), typeof(f₃))}(J₁₁, J₂₁, J₃₁, J₁₂, J₂₂, J₃₂, J₁₃, J₂₃, J₃₃, L₁, L₂, L₃, f₁, f₂, f₂) + AttitudeParameters{F}( + J₁₁, + J₂₁, + J₃₁, + J₁₂, + J₂₂, + J₃₂, + J₁₃, + J₂₃, + J₃₃, + L₁, + L₂, + L₃, + f₁, + f₂, + f₃, + ) where {F} = + new{F}(J₁₁, J₂₁, J₃₁, J₁₂, J₂₂, J₃₂, J₁₃, J₂₃, J₃₃, L₁, L₂, L₃, f₁, f₂, f₃) + AttitudeParameters( + J₁₁, + J₂₁, + J₃₁, + J₁₂, + J₂₂, + J₃₂, + J₁₃, + J₂₃, + J₃₃, + L₁, + L₂, + L₃, + f₁, + f₂, + f₃, + ) = new{ + promote_type( + typeof(J₁₁), + typeof(J₂₁), + typeof(J₃₁), + typeof(J₁₂), + typeof(J₂₂), + typeof(J₃₂), + typeof(J₁₃), + typeof(J₂₃), + typeof(J₃₃), + typeof(L₁), + typeof(L₂), + typeof(L₃), + typeof(f₁), + typeof(f₂), + typeof(f₃), + ), + }( + J₁₁, + J₂₁, + J₃₁, + J₁₂, + J₂₂, + J₃₂, + J₁₃, + J₂₃, + J₃₃, + L₁, + L₂, + L₃, + f₁, + f₂, + f₂, + ) end system(::AttitudeParameters, args...; kwargs...) = AttitudeSystem(args...; kwargs...) @@ -78,7 +164,12 @@ spherical planet. model = Attitude() ``` """ -@memoize function AttitudeSystem(; stm=false, name=:Attitude, defaults=Pair{ModelingToolkit.Num,<:Number}[], kwargs...) +@memoize function AttitudeSystem(; + stm = false, + name = :Attitude, + defaults = Pair{ModelingToolkit.Num,<:Number}[], + kwargs..., +) @variables t @variables (q(t))[1:4] [description = "scalar-last attitude quaternion"] @@ -107,20 +198,17 @@ model = Attitude() -ω[2] ω[1] 0 ] - eqs = vcat( - δ.(q) .~ (1 // 2) * (A * q), - δ.(ω) .~ (-inv(J) * ωx * J * ω + inv(J) * L + f) - ) + eqs = + vcat(δ.(q) .~ (1 // 2) * (A * q), δ.(ω) .~ (-inv(J) * ωx * J * ω + inv(J) * L + f)) if stm - @variables (Φ(t))[1:7, 1:7] - Φ = Symbolics.scalarize(Φ) - M = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(q, ω)) + @variables (Φ(t))[1:7, 1:7] [description = "state transition matrix estimate"] + A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v)) - LHS = map(δ, Φ) - RHS = M * Φ + LHS = δ.(Φ) + RHS = A * Φ - eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in eachindex(LHS)]) + eqs = vcat(eqs, vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])) end if string(name) == "Attitude" && stm @@ -131,9 +219,25 @@ model = Attitude() if stm append!(defaults, vec(Φ .=> I(7))) - return ODESystem(eqs, t, vcat(q, ω, vec(Φ)), vcat(vec(J), L, f); name=name, defaults=defaults, kwargs...) + return ODESystem( + eqs, + t, + vcat(q, ω, vec(Φ)), + vcat(vec(J), L, f); + name = name, + defaults = defaults, + kwargs..., + ) else - return ODESystem(eqs, t, vcat(q, ω), vcat(vec(J), L, f); name=name, defaults=defaults, kwargs...) + return ODESystem( + eqs, + t, + vcat(q, ω), + vcat(vec(J), L, f); + name = name, + defaults = defaults, + kwargs..., + ) end end @@ -155,12 +259,14 @@ let u = randn(7), p = randn(15), t = NaN # time invariant end ``` """ -@memoize function AttitudeFunction(; stm=false, name=:Attitude, kwargs...) - defaults = (; jac=true) +@memoize function AttitudeFunction(; stm = false, name = :Attitude, kwargs...) + defaults = (; jac = true) options = merge(defaults, kwargs) - sys = complete(AttitudeSystem(; stm=stm, name=name); split=false) + sys = complete(AttitudeSystem(; stm = stm, name = name); split = false) return ODEFunction{true,SciMLBase.FullSpecialize}( - sys, ModelingToolkit.unknowns(sys), ModelingToolkit.parameters(sys); - options... + sys, + ModelingToolkit.unknowns(sys), + ModelingToolkit.parameters(sys); + options..., ) end diff --git a/lib/AstrodynamicalModels/src/CR3BP.jl b/lib/AstrodynamicalModels/src/CR3BP.jl index 2832adbb..88009f38 100644 --- a/lib/AstrodynamicalModels/src/CR3BP.jl +++ b/lib/AstrodynamicalModels/src/CR3BP.jl @@ -2,8 +2,7 @@ # Circular Restricted Three-body Problem models # -@doc CartesianState -const CR3BState = CartesianState +@doc CartesianState const CR3BState = CartesianState """ A paremeter vector for CR3BP dynamics. @@ -44,7 +43,12 @@ systems in our solar system. model = CR3BSystem(; stm=true) ``` """ -@memoize function CR3BSystem(; stm=false, name=:CR3B, defaults=Pair{ModelingToolkit.Num,<:Number}[], kwargs...) +@memoize function CR3BSystem(; + stm = false, + name = :CR3B, + defaults = Pair{ModelingToolkit.Num,<:Number}[], + kwargs..., +) @parameters t μ @variables x(t) y(t) z(t) ẋ(t) ẏ(t) ż(t) @@ -54,20 +58,31 @@ model = CR3BSystem(; stm=true) eqs = vcat( δ.(r) .~ v, - δ(ẋ) ~ x + 2ẏ - (μ * (x + μ - 1) * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3)) - ((x + μ) * (sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ)), - δ(ẏ) ~ y - (2ẋ) - (y * (μ * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) + (sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ))), - δ(ż) ~ z * (-μ * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) - ((sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ))) + δ(ẋ) ~ + x + 2ẏ - (μ * (x + μ - 1) * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3)) - + ((x + μ) * (sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ)), + δ(ẏ) ~ + y - (2ẋ) - ( + y * ( + μ * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) + + (sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ) + ) + ), + δ(ż) ~ + z * ( + -μ * (sqrt(y^2 + z^2 + (x + μ - 1)^2)^-3) - + ((sqrt(y^2 + z^2 + (x + μ)^2)^-3) * (1 - μ)) + ), ) if stm @variables (Φ(t))[1:6, 1:6] [description = "state transition matrix estimate"] - Φ = Symbolics.scalarize(Φ) A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v)) - LHS = map(δ, Φ) - RHS = map(simplify, A * Φ) + LHS = δ.(Φ) + RHS = A * Φ - eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in eachindex(LHS)]) + eqs = vcat(eqs, vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])) end if string(name) == "CR3B" && stm @@ -79,14 +94,23 @@ model = CR3BSystem(; stm=true) if stm append!(defaults, vec(Φ .=> I(6))) return ODESystem( - eqs, t, vcat(r, v, vec(Φ)), [μ]; - name=modelname, - defaults=defaults, - kwargs... + eqs, + t, + vcat(r, v, vec(Φ)), + [μ]; + name = modelname, + defaults = defaults, + kwargs..., ) else return ODESystem( - eqs, t, vcat(r, v), [μ]; name=modelname, defaults=defaults, kwargs... + eqs, + t, + vcat(r, v), + [μ]; + name = modelname, + defaults = defaults, + kwargs..., ) end end @@ -113,13 +137,15 @@ let u = randn(6), p = randn(1), t = 0 end ``` """ -@memoize function CR3BFunction(; stm=false, name=:CR3B, kwargs...) - defaults = (; jac=true) +@memoize function CR3BFunction(; stm = false, name = :CR3B, kwargs...) + defaults = (; jac = true) options = merge(defaults, kwargs) - sys = complete(CR3BSystem(; stm=stm, name=name); split=false) + sys = complete(CR3BSystem(; stm = stm, name = name); split = false) return ODEFunction{true,SciMLBase.FullSpecialize}( - sys, ModelingToolkit.unknowns(sys), ModelingToolkit.parameters(sys); - options... + sys, + ModelingToolkit.unknowns(sys), + ModelingToolkit.parameters(sys); + options..., ) end @@ -127,12 +153,22 @@ end An `Orbit` which exists within CR3BP dynamics. """ const CR3BOrbit = Orbit{<:CR3BState,<:CR3BParameters} -AstrodynamicalModels.CR3BOrbit(state::AbstractVector, parameters::AbstractVector) = Orbit(CR3BState(state), CR3BParameters(parameters)) -AstrodynamicalModels.CR3BOrbit(; state::AbstractVector, parameters::AbstractVector) = Orbit(CR3BState(state), CR3BParameters(parameters)) +AstrodynamicalModels.CR3BOrbit(state::AbstractVector, parameters::AbstractVector) = + Orbit(CR3BState(state), CR3BParameters(parameters)) +AstrodynamicalModels.CR3BOrbit(; state::AbstractVector, parameters::AbstractVector) = + Orbit(CR3BState(state), CR3BParameters(parameters)) """ Return an `ODEProblem` for the provided CR3B system. """ CR3BProblem(u0, tspan, p; kwargs...) = ODEProblem(CR3BFunction(), u0, tspan, p; kwargs...) -CR3BProblem(orbit::AstrodynamicalOrbit, tspan::Union{<:Tuple,<:AbstractArray}; kwargs...) = ODEProblem(CR3BFunction(), AstrodynamicalModels.state(orbit), tspan, AstrodynamicalModels.parameters(orbit); kwargs...) -CR3BProblem(orbit::AstrodynamicalOrbit, Δt; kwargs...) = CR3BProblem(orbit, (zero(Δt), δt); kwargs...) \ No newline at end of file +CR3BProblem(orbit::AstrodynamicalOrbit, tspan::Union{<:Tuple,<:AbstractArray}; kwargs...) = + ODEProblem( + CR3BFunction(), + AstrodynamicalModels.state(orbit), + tspan, + AstrodynamicalModels.parameters(orbit); + kwargs..., + ) +CR3BProblem(orbit::AstrodynamicalOrbit, Δt; kwargs...) = + CR3BProblem(orbit, (zero(Δt), δt); kwargs...) \ No newline at end of file diff --git a/lib/AstrodynamicalModels/src/NBP.jl b/lib/AstrodynamicalModels/src/NBP.jl index 8e32100c..c5afa2fa 100644 --- a/lib/AstrodynamicalModels/src/NBP.jl +++ b/lib/AstrodynamicalModels/src/NBP.jl @@ -35,7 +35,13 @@ That's about right for a model in a package called model = NBSystem(9) ``` """ -@memoize function NBSystem(N::Int; stm=false, name=:NBP, defaults=Pair{ModelingToolkit.Num,<:Number}[], kwargs...) +@memoize function NBSystem( + N::Int; + stm = false, + name = :NBP, + defaults = Pair{ModelingToolkit.Num,<:Number}[], + kwargs..., +) N > 0 || throw(ArgumentError("`N` must be a number greater than zero!")) T = N * 6 + (N * 6)^2 @@ -43,20 +49,20 @@ model = NBSystem(9) @variables (x(t))[1:N] (y(t))[1:N] (z(t))[1:N] (ẋ(t))[1:N] (ẏ(t))[1:N] (ż(t))[1:N] δ = Differential(t) - r = [[x[i], y[i], z[i]] for i in 1:N] - v = [[ẋ[i], ẏ[i], ż[i]] for i in 1:N] + r = [[x[i], y[i], z[i]] for i = 1:N] + v = [[ẋ[i], ẏ[i], ż[i]] for i = 1:N] - poseqs = reduce(vcat, [ - δ.(r[i]) .~ v[i] - for i ∈ 1:N - ]) + poseqs = reduce(vcat, [δ.(r[i]) .~ v[i] for i ∈ 1:N]) - veleqs = reduce(vcat, [ - δ.(v[i]) .~ sum(( - ((m[i] * m[j]) / norm(r[j] .- r[i])^3 * (r[j] .- r[i])) .* (j == N ? G / m[i] : 1) - for j ∈ 1:N if j ≢ i - )) for i ∈ 1:N - ]) + veleqs = reduce( + vcat, + [ + δ.(v[i]) .~ sum(( + ((m[i] * m[j]) / norm(r[j] .- r[i])^3 * (r[j] .- r[i])) .* + (j == N ? G / m[i] : 1) for j ∈ 1:N if j ≢ i + )) for i ∈ 1:N + ], + ) eqs = vcat(poseqs, veleqs) @@ -70,14 +76,13 @@ model = NBSystem(9) """ end - @variables (Φ(t))[1:6N, 1:6N] [description = "state transition matrix estimate"] - Φ = Symbolics.scalarize(Φ) - A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r..., v...)) + @variables (Φ(t))[1:6, 1:6] [description = "state transition matrix estimate"] + A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v)) - LHS = map(δ, Φ) - RHS = map(simplify, A * Φ) + LHS = δ.(Φ) + RHS = A * Φ - eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in eachindex(LHS)]) + eqs = vcat(eqs, vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])) end if string(name) == "NBP" && stm @@ -89,17 +94,23 @@ model = NBSystem(9) if stm append!(defaults, vec(Φ .=> I(6N))) return ODESystem( - eqs, t, vcat(r..., v..., Φ...), vcat(G, m...); - name=modelname, - defaults=defaults, - kwargs... + eqs, + t, + vcat(r..., v..., Φ...), + vcat(G, m...); + name = modelname, + defaults = defaults, + kwargs..., ) else return ODESystem( - eqs, t, vcat(r..., v...), vcat(G, m...); - name=modelname, - defaults=defaults, - kwargs... + eqs, + t, + vcat(r..., v...), + vcat(G, m...); + name = modelname, + defaults = defaults, + kwargs..., ) end end @@ -148,8 +159,8 @@ let u = randn(3*6), p = randn(1 + 3), t = 0 end ``` """ -@memoize function NBFunction(N::Int; stm=false, name=:R2BP, kwargs...) - defaults = (; jac=false) +@memoize function NBFunction(N::Int; stm = false, name = :R2BP, kwargs...) + defaults = (; jac = false) options = merge(defaults, kwargs) if N ≥ 2 && stm && options.jac @warn """ @@ -159,9 +170,11 @@ end calculations! Consider setting `jac=false`, `stm=false`, or both. """ end - sys = complete(NBSystem(N; stm=stm, name=name); split=false) + sys = complete(NBSystem(N; stm = stm, name = name); split = false) return ODEFunction{true,SciMLBase.FullSpecialize}( - sys, ModelingToolkit.unknowns(sys), ModelingToolkit.parameters(sys); - options... + sys, + ModelingToolkit.unknowns(sys), + ModelingToolkit.parameters(sys); + options..., ) end diff --git a/lib/AstrodynamicalModels/src/R2BP.jl b/lib/AstrodynamicalModels/src/R2BP.jl index 973b197a..7c60a4bd 100644 --- a/lib/AstrodynamicalModels/src/R2BP.jl +++ b/lib/AstrodynamicalModels/src/R2BP.jl @@ -2,8 +2,7 @@ # Restricted Two-body Problem models # -@doc CartesianState -const R2BState = CartesianState +@doc CartesianState const R2BState = CartesianState """ A parameter vector for R2BP dynamics. @@ -45,7 +44,12 @@ spacecraft orbiting Earth. model = R2BSystem() ``` """ -@memoize function R2BSystem(; stm=false, name=:R2B, defaults=Pair{ModelingToolkit.Num,<:Number}[], kwargs...) +@memoize function R2BSystem(; + stm = false, + name = :R2B, + defaults = Pair{ModelingToolkit.Num,<:Number}[], + kwargs..., +) @parameters t μ @variables x(t) y(t) z(t) ẋ(t) ẏ(t) ż(t) @@ -53,20 +57,16 @@ model = R2BSystem() r = [x, y, z] v = [ẋ, ẏ, ż] - eqs = vcat( - δ.(r) .~ v, - δ.(v) .~ -μ .* (r ./ norm(r)^3) - ) + eqs = vcat(δ.(r) .~ v, δ.(v) .~ -μ .* (r ./ norm(r)^3)) if stm @variables (Φ(t))[1:6, 1:6] [description = "state transition matrix estimate"] - Φ = Symbolics.scalarize(Φ) A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(r, v)) - LHS = map(δ, Φ) - RHS = map(simplify, A * Φ) + LHS = δ.(Φ) + RHS = A * Φ - eqs = vcat(eqs, [LHS[i] ~ RHS[i] for i in 1:length(LHS)]) + eqs = vcat(eqs, vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])) end if string(name) == "R2B" && stm @@ -78,14 +78,23 @@ model = R2BSystem() if stm append!(defaults, vec(Φ .=> I(6))) return ODESystem( - eqs, t, vcat(r, v, vec(Φ)), [μ]; - name=modelname, - defaults=defaults, - kwargs... + eqs, + t, + vcat(r, v, vec(Φ)), + [μ]; + name = modelname, + defaults = defaults, + kwargs..., ) else return ODESystem( - eqs, t, vcat(r, v), [μ]; name=modelname, defaults=defaults, kwargs... + eqs, + t, + vcat(r, v), + [μ]; + name = modelname, + defaults = defaults, + kwargs..., ) end end @@ -112,13 +121,15 @@ let u = randn(6), p = randn(1), t = 0 end ``` """ -@memoize function R2BFunction(; stm=false, name=:R2B, kwargs...) - defaults = (; jac=true) +@memoize function R2BFunction(; stm = false, name = :R2B, kwargs...) + defaults = (; jac = true) options = merge(defaults, kwargs) - sys = complete(R2BSystem(; stm=stm, name=name); split=false) + sys = complete(R2BSystem(; stm = stm, name = name); split = false) return ODEFunction{true,SciMLBase.FullSpecialize}( - sys, ModelingToolkit.unknowns(sys), ModelingToolkit.parameters(sys); - options... + sys, + ModelingToolkit.unknowns(sys), + ModelingToolkit.parameters(sys); + options..., ) end @@ -131,5 +142,16 @@ const R2BOrbit = Orbit{<:CartesianState,<:R2BParameters} Return an `ODEProblem` for the provided R2B system. """ R2BProblem(u0, tspan, p; kwargs...) = ODEProblem(R2BFunction(), u0, tspan, p; kwargs...) -R2BProblem(orbit::AstrodynamicalOrbit{<:CartesianState}, tspan::Union{<:Tuple,<:AbstractArray}; kwargs...) = ODEProblem(R2BFunction(), AstrodynamicalModels.state(orbit), tspan, AstrodynamicalModels.parameters(orbit); kwargs...) -R2BProblem(orbit::AstrodynamicalOrbit{<:CartesianState}, Δt; kwargs...) = R2BProblem(orbit, (zero(Δt), δt); kwargs...) +R2BProblem( + orbit::AstrodynamicalOrbit{<:CartesianState}, + tspan::Union{<:Tuple,<:AbstractArray}; + kwargs..., +) = ODEProblem( + R2BFunction(), + AstrodynamicalModels.state(orbit), + tspan, + AstrodynamicalModels.parameters(orbit); + kwargs..., +) +R2BProblem(orbit::AstrodynamicalOrbit{<:CartesianState}, Δt; kwargs...) = + R2BProblem(orbit, (zero(Δt), δt); kwargs...)