AstrodynamicalModels.jl

Common models within astrodynamics!

Overview

This package extends ModelingToolkit to represent common astrodynamical models. All available models are shown on the Docstrings page. Consult the Models pages for more detail about each model in this package!

Usage

If you're familiar with ModelingToolkit.jl, then you'll be able to use this package! Some AstrodynamicalModels-specific usage instructions are provided here. Please don't be shy about making Discourse posts, or filing issues on GitHub!

Installation & Setup

This package can be installed just like any other registered Julia package.

# To install wherever Julia code runs...
import Pkg
Pkg.add("AstrodynamicalModels") # or ]add AstrodynamicalModels in Julia's REPL

To load the package, simply enter using AstrodynamicalModels.

julia> using AstrodynamicalModels

Retrieving a Model

Each model within this package is implemented with a function – each function returns some AbstractSystem from ModelingToolkit.jl. Typically, this will be an ODESystem. If you're worried about overhead from calling each function every time you need a particular model, don't! Each function is implemented with @memoize, so all results are cached the first time you call a model's function with a particular function signature.

julia> R2BPModel = R2BP() # Restricted Two-body Problem dynamicsERROR: UndefVarError: `R2BP` not defined
julia> CR3BPModel = CR3BP() # Circular Restricted Three-body Problem dynamicsERROR: UndefVarError: `CR3BP` not defined
julia> CR3BPModelWithSTM = CR3BP(; stm=true) # Optionally include state transition matrix dynamicsERROR: UndefVarError: `CR3BP` not defined

Using a Model

To actually use each model, you probably also want to load ModelingToolkit (and any other SciML packages of your choice).

julia> using ModelingToolkit

Now you can use any method defined for ModelingToolkit.AbstractSystem instances. Once again, the ModelingToolkit Documentation are the best place to learn how to interact with AbstractSystem instances! Some quick examples are shown below.

Check the Equations of Motion

julia> eqs = equations(R2BSystem())6-element Vector{Equation}:
 Differential(t)(x(t)) ~ ẋ(t)
 Differential(t)(y(t)) ~ ẏ(t)
 Differential(t)(z(t)) ~ ż(t)
 Differential(t)(ẋ(t)) ~ (-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)
 Differential(t)(ẏ(t)) ~ (-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)
 Differential(t)(ż(t)) ~ (-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)

List the States and Parameters

julia> x = states(R2BSystem())6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 x(t)
 y(t)
 z(t)
 ẋ(t)
 ẏ(t)
 ż(t)
julia> p = parameters(R2BSystem())1-element Vector{SymbolicUtils.BasicSymbolic{Real}}: μ

Calculate the Jacobian

julia> J = calculate_jacobian(R2BSystem())6×6 Matrix{Num}:
                                                                                                                                                                     0  …  1  0  0
                                                                                                                                                                     0     0  1  0
                                                                                                                                                                     0     0  0  1
 (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3x(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))     0  0  0
                                                        -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*x(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))     0  0  0
                                                        -3x(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))  …  0  0  0

Generate Code to Replicate the Model

julia> print(build_function(R2BSystem()))function ()
    #= /juliateam/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:373 =#
    #= /juliateam/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:374 =#
    #= /juliateam/.julia/packages/SymbolicUtils/r1pzW/src/code.jl:375 =#
    ODESystem(0x0000000000000006, Equation[Differential(t)(x(t)) ~ ẋ(t), Differential(t)(y(t)) ~ ẏ(t), Differential(t)(z(t)) ~ ż(t), Differential(t)(ẋ(t)) ~ (-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ẏ(t)) ~ (-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ż(t)) ~ (-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), z(t), ẋ(t), ẏ(t), ż(t)], SymbolicUtils.BasicSymbolic{Real}[μ], nothing, Dict{Any, Any}(:ż => ż(t), :μ => μ, :y => y(t), :ẏ => ẏ(t), :ẋ => ẋ(t), :z => z(t), :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}((Num[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3x(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3z(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*x(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*y(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*z(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3x(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3z(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0], (false, false))), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :R2B, ODESystem[], Dict{Any, Any}(), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Equation[], Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing)
end

Generate Code which Implements the Dynamics

julia> print(R2BFunction())ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#k#503"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x054b07f6, 0xa4ef70cd, 0xd58d1cd7, 0x3cda4354, 0x9d6c242b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x052b3c68, 0x124778d4, 0xd77b848d, 0xb960d875, 0xf486362f), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#___jac#509"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x67420626, 0x1a4595aa, 0x37e5bbe7, 0xf10f7fc6, 0xf0e5b87a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a0dff05, 0xbcb247fb, 0xa2658adb, 0xedbf2c81, 0xe293a4dd), Nothing}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.var"#50963#generated_observed#513"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}, Nothing, ODESystem}(ModelingToolkit.var"#k#503"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x054b07f6, 0xa4ef70cd, 0xd58d1cd7, 0x3cda4354, 0x9d6c242b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x052b3c68, 0x124778d4, 0xd77b848d, 0xb960d875, 0xf486362f), Nothing}}(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x054b07f6, 0xa4ef70cd, 0xd58d1cd7, 0x3cda4354, 0x9d6c242b), Nothing}(nothing), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x052b3c68, 0x124778d4, 0xd77b848d, 0xb960d875, 0xf486362f), Nothing}(nothing)), LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, ModelingToolkit.var"#___jac#509"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x67420626, 0x1a4595aa, 0x37e5bbe7, 0xf10f7fc6, 0xf0e5b87a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a0dff05, 0xbcb247fb, 0xa2658adb, 0xedbf2c81, 0xe293a4dd), Nothing}}(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x67420626, 0x1a4595aa, 0x37e5bbe7, 0xf10f7fc6, 0xf0e5b87a), Nothing}(nothing), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9a0dff05, 0xbcb247fb, 0xa2658adb, 0xedbf2c81, 0xe293a4dd), Nothing}(nothing)), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, ModelingToolkit.var"#50963#generated_observed#513"{Bool, ODESystem, Dict{Any, Any}, Vector{SymbolicUtils.BasicSymbolic{Real}}}(Core.Box(nothing), Core.Box(nothing), false, Core.Box(Equation[]), ODESystem(0x0000000000000006, Equation[Differential(t)(x(t)) ~ ẋ(t), Differential(t)(y(t)) ~ ẏ(t), Differential(t)(z(t)) ~ ż(t), Differential(t)(ẋ(t)) ~ (-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ẏ(t)) ~ (-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ż(t)) ~ (-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), z(t), ẋ(t), ẏ(t), ż(t)], SymbolicUtils.BasicSymbolic{Real}[μ], nothing, Dict{Any, Any}(:ż => ż(t), :μ => μ, :y => y(t), :ẏ => ẏ(t), :ẋ => ẋ(t), :z => z(t), :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}((Num[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3x(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3z(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*x(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*y(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*z(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3x(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3z(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0], (false, false))), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :R2B, ODESystem[], Dict{Any, Any}(), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Equation[], Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing), Dict{Any, Any}(), SymbolicUtils.BasicSymbolic{Real}[μ]), nothing, ODESystem(0x0000000000000006, Equation[Differential(t)(x(t)) ~ ẋ(t), Differential(t)(y(t)) ~ ẏ(t), Differential(t)(z(t)) ~ ż(t), Differential(t)(ẋ(t)) ~ (-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ẏ(t)) ~ (-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3), Differential(t)(ż(t)) ~ (-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3)], t, SymbolicUtils.BasicSymbolic{Real}[x(t), y(t), z(t), ẋ(t), ẏ(t), ż(t)], SymbolicUtils.BasicSymbolic{Real}[μ], nothing, Dict{Any, Any}(:ż => ż(t), :μ => μ, :y => y(t), :ẏ => ẏ(t), :ẋ => ẋ(t), :z => z(t), :x => x(t)), Any[], Equation[], Base.RefValue{Vector{Num}}(Num[]), Base.RefValue{Any}((Num[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3x(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3z(t)*((-x(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*x(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*y(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3((-y(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*z(t)*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0; -3x(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) -3y(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) (-μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^3) - 3z(t)*((-z(t)*μ) / (sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t)))^6))*sqrt(abs2(y(t)) + abs2(x(t)) + abs2(z(t))) 0 0 0], (false, false))), Base.RefValue{Any}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), Base.RefValue{Matrix{Num}}(Matrix{Num}(undef, 0, 0)), :R2B, ODESystem[], Dict{Any, Any}(), nothing, nothing, nothing, ModelingToolkit.SymbolicContinuousCallback[ModelingToolkit.SymbolicContinuousCallback(Equation[], Equation[])], ModelingToolkit.SymbolicDiscreteCallback[], nothing, nothing, nothing, nothing, false, nothing, nothing, nothing, nothing))

Generate C/C++ and MATLAB Code

julia> print(build_function([eq.rhs for eq in equations(R2BSystem())], states(R2BSystem()), parameters(R2BSystem()); target=Symbolics.CTarget()))#include <math.h>
void diffeqf(double* du, const double* RHS1, const double* RHS2) {
  du[0] = RHS1[3];
  du[1] = RHS1[4];
  du[2] = RHS1[5];
  du[3] = (-1 * RHS1[0] * RHS2[0]) / pow(sqrt(abs2(RHS1[1]) + abs2(RHS1[0]) + abs2(RHS1[2])), 3);
  du[4] = (-1 * RHS1[1] * RHS2[0]) / pow(sqrt(abs2(RHS1[1]) + abs2(RHS1[0]) + abs2(RHS1[2])), 3);
  du[5] = (-1 * RHS1[2] * RHS2[0]) / pow(sqrt(abs2(RHS1[1]) + abs2(RHS1[0]) + abs2(RHS1[2])), 3);
}
julia> print(build_function([eq.rhs for eq in equations(R2BSystem())], states(R2BSystem()), parameters(R2BSystem()); target=Symbolics.MATLABTarget()))diffeqf = @(internal_var___t,internal_var___u) [ internal_var___u(4); internal_var___u(5); internal_var___u(6); (-1 * internal_var___u(1) * internal_var___p(1)) / sqrt(abs2(internal_var___u(2)) + abs2(internal_var___u(1)) + abs2(internal_var___u(3))) ^ 3; (-1 * internal_var___u(2) * internal_var___p(1)) / sqrt(abs2(internal_var___u(2)) + abs2(internal_var___u(1)) + abs2(internal_var___u(3))) ^ 3; (-1 * internal_var___u(3) * internal_var___p(1)) / sqrt(abs2(internal_var___u(2)) + abs2(internal_var___u(1)) + abs2(internal_var___u(3))) ^ 3; ];

If you're interested in learning a bit about each astrodynamical model, or you'd like more specific examples which show how to use each model, consult the Models pages!