Use Case: Custom propagated orbit

This example was generated on 2024-11-18T22:29:37.076.

Once the general structure of the FrameSystem is understood, we can pass to a use case in which we want to build and exploit our frame system to perform computations inserting a custom orbit. This becomes especially crucial in complex cases like Trajectory Optimization and Navigation Analysis.

In such scenarios, the trajectory is under design, and the trajectory information might not be completely available. Moreover, in these cases derivatives could be required for various quantities such as time, states, and parameters.

In this context, we shall remember that FrameTransformations is able to perform operations, including AD, on the frames independent variable, e.g. only time.

A proper orbit representation is essential to avoid perturbation confusion and ensure proper custom orbit handling.

Frame system setup

First of all, a new FrameSystem shall be created.

using FrameTransformations

G = FrameSystem{4,Float64}()

add_axes!(G, :ICRF, 1)
add_point!(G, :Earth, 399, 1)

Note

As the orbit considered in this example is custom, there is no requirement to load ephemeris.

Custom orbit model

A custom orbit model is then required. In this case a simple two-body problem is considered but more complex models can be used. We interface with SciML ODE solvers to compute the trajectory solution.

using OrdinaryDiffEq
using DiffEqBase

Now we create the problem with a given initial condition, time-span and gravitational parameter.

function two_body!(du, u, p, t)
    μ = p[1]
    r = sqrt(u[1]^2 + u[2]^2 + u[3]^2)
    r3 = r^3
    du[1] = u[4]
    du[2] = u[5]
    du[3] = u[6]
    du[4] = -μ * u[1] / r3
    du[5] = -μ * u[2] / r3
    du[5] = -μ * u[2] / r3
    nothing
end

u0 = [7000.0, 0.0, 0.0, 0.0, 8.0, 0.1]
p = [398600.0,]
tspan = (0, 10 * 86400.0)

prob = ODEProblem(two_body!, u0, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 864000.0)
u0: 6-element Vector{Float64}:
 7000.0
    0.0
    0.0
    0.0
    8.0
    0.1

It is possible to create a new node using directly the solution of the ODE problem or before the problem is solved, using the integrator:

integrator = DiffEqBase.__init(prob, Vern9(), abstol=1e-14, reltol=1e-12)
add_point_dynamical!(G, :SC, -1, 399, 1, t -> @views(integrator.sol(t)[1:3]), t -> integrator.sol(t))

Now we can compute the trajectory:

solve!(integrator)

retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 767-element Vector{Float64}:
      0.0
      0.002100410678471327
      0.0047599230055715425
      0.008955696574592512
      0.014405191722968853
      0.02256127192121507
      0.03376908947224431
      0.049589055530138385
      0.0698982062387515
      0.09881641196938389
      ⋮
 806236.5951903301
 813554.5837279536
 821095.7757502404
 828936.367091165
 837208.97355692
 845655.5468850702
 853693.9091403671
 862493.6137650727
 864000.0
u: 767-element Vector{Vector{Float64}}:
 [7000.0, 0.0, 0.0, 0.0, 8.0, 0.1]
 [6999.999999982056, 0.01680328542775626, 0.0002100410678471327, -1.7086197886483454e-5, 7.999999999979493, 0.1]
 [6999.999999907846, 0.03807938404440524, 0.0004759923005571542, -3.8720516530803896e-5, 7.9999999998946825, 0.1]
 [6999.99999967378, 0.07164557259562712, 0.000895569657459251, -7.285185009298929e-5, 7.999999999627178, 0.1]
 [6999.999999155986, 0.11524153377911914, 0.0014405191722968852, -0.00011718182490732019, 7.999999999035414, 0.1]
 [6999.999997929675, 0.1804901753519266, 0.002256127192121507, -0.00018352904054244176, 7.999999997633916, 0.1]
 [6999.999995361795, 0.2701527157182867, 0.003376908947224431, -0.000274701205297076, 7.999999994699195, 0.1]
 [6999.999989998091, 0.3967124440521604, 0.004958905553013839, -0.0004033917861508917, 7.999999988569248, 0.1]
 [6999.999980127922, 0.5591856493808606, 0.00698982062387515, -0.0005686005096037746, 7.999999977289055, 0.1]
 [6999.999960283645, 0.7905312942599749, 0.00988164119693839, -0.0008038412593625989, 7.999999954609882, 0.1]
 ⋮
 [27780.683963325602, 42359.01446895586, 80623.65951903291, -0.8025971884101836, 0.7920170040605534, 0.1]
 [21593.404005364973, 47617.06581630905, 81355.45837279527, -0.8843405807263419, 0.643265710217164, 0.1]
 [14685.747014284721, 51862.035932975115, 82109.57757502396, -0.9434414826071104, 0.48150113994807625, 0.1]
 [7133.1289576091995, 54962.451067493384, 82893.63670911641, -0.9787751475642218, 0.30899747628535834, 0.1]
 [-1024.0952950747117, 56763.854488823155, 83720.89735569191, -0.9888316844294007, 0.12684156398017055, 0.1]
 [-9325.026533480774, 57062.54147480205, 84565.55468850693, -0.9723707350002112, -0.055126342355960874, 0.1]
 [-17000.28781860142, 55946.09286142069, 85369.39091403662, -0.9337249254505049, -0.22127499568585487, 0.1]
 [-24945.126440130833, 53237.49061776463, 86249.36137650716, -0.8681287897346082, -0.3921808825049796, 0.1]
 [-26242.78437442029, 52625.65098316583, 86399.9999999999, -0.8546412436047252, -0.4200754020972222, 0.1]

From now on we can use the frame system with the :SC custom node:

vector6(G, :Earth, :SC, :ICRF, 123456.0)

6-element StaticArraysCore.SVector{6, Float64} with indices SOneTo(6):
 -12000.809592690548
  10523.77670590901
  12345.600000000011
     -2.7588011391226344
     -2.247097800159169
      0.1

This page was generated using Literate.jl.