 # Julia now has a native SGP4/SDP4 implementation

Hi guys!

I am happy to announce that I finished the first version of the SGP4/SDP4 implementation 100% written in Julia. The source code is part of SatelliteToolbox.jl. Notice that I did not tag yet, so you need to use `master`.

As for the correctness of the algorithm, all the test TLEs available in  has been added into the test set, and everything passes. Thus, I am not expecting any major problems. But I will be very happy if you can help me to test in other scenarios and report weird behaviors Notice that I have not yet implemented the verifications for things like the satellite crashed or numerical errors arising in veeeery long integrations (this is easy and will be done soon).

Furthermore, I decided to make a “readable” code instead of the most performant possible. Because this will be used to teach students about orbit propagators. So, for example, the FORTRAN code:

``````    SING=SIN(OMEGAO)
COSG=COS(OMEGAO)
TSI=1./(AODP-S4)
ETA=AODP*EO*TSI
ETASQ=ETA*ETA
EETA=EO*ETA
PSISQ=ABS(1.-ETASQ)
COEF=QOMS24*TSI**4
COEF1=COEF/PSISQ**3.5
C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
1         CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
C1=BSTAR*C2
``````

was translated into

``````    aux0 = abs(1-η^2)
aux1 = aux0^(-7/2)
aux2 = ξ^4*all_0*β_0^2*aux1

C2 = QOMS2T*ξ^4*nll_0*aux1*
( all_0*( 1 + (3/2)η^2 + 4e_0*η + e_0*η^3) +
3/2*(k_2*ξ)/aux0*(-1/2 + (3/2)θ²)*(8 + 24η^2 + 3η^4) )

C1 = bstar*C2
``````

(give or take some lines of code…)

I know that FORTRAN version should be faster, but the Julia one is more easily related to the equations we have in the books.

Finally, I want the code to be very documented about the theory. So that, in the far far far away future, when someone needs to rewrite in another language, he/she will not have the problems I had… Hence, if you are an export in Astrodynamics and want to add comments about the equations in the code (specially the deep space bits…), be my guest : Vallado, D. A., Crawford, P., Hujsak, R., Kelso, T. S (2006). Revisiting Spacetrack Report #3: Rev1. AIAA.

10 Likes

This is awesome, well done. I’d be really interested to see benchmarks against the variety of other implementations that are out there.

1 Like

Me too! But be aware that my coding expertise in Julia is very limited yet! I will start shortly a campaign to optimize the code (help needed!). I am pretty sure there is a lot of thing that can be done.

2 Likes

Count me in I would like to integrate your propagator with the generic propagator interface I am developing at some point.

2 Likes

I am not an astrophysicist, but still find these announcements very interesting. Short explanations of the acronyms would make it easier to understand (eg here I understood that these are commonly used perturbation models for orbital mechanics).

3 Likes

Sorry, my fault!

SGP4 (Simplified General Perturbations 4) is common orbit propagator. As input, you can get the Two Line Elements (TLE, which is basically all the 6 orbit elements + some additional info written in two lines of text) and the SGP4 will predict where those object will be in the future.

The original SGP4 accounts for the Earth gravity and atmospheric drag only, but for objects with a very high orbit (orbital period higher than 225 min.) we need to consider perturbations effects from the Deep Space, specially the Sun and the Moon. Those perturbations are coded in what we called SDP4, which is the SGP4 + the perturbations.

Modern algorithms get the input elements and automatically choose whether we should add the deep space computations or not. This is why we are doing here.

Vallado (which is the author of THE book of Astrodynamics) has in his website the SGP4/SDP4 written in some languages: FORTRAN (which is the original), C++, Pascal, Matlab, C#. It is also possible to find the algorithm in Python and Java, however I am not sure if those are native implementations or wrappers of a C library.

Now, Julia has its own SGP4/SDP4 algorithm natively written. It will be very interesting to compare the speed against the very well tested and debugged FORTRAN version.

3 Likes

Well, at least against Matlab it is pretty fast:

# Julia

``````julia> tle = read_tle("/Users/ronan.arraes/tmp/teste.tle")
julia> orbp = init_orbit_propagator(Val{:sgp4}, tle)
julia> r,v = sgp4!(orbp.sgp4d, 0)
julia> @time r,v = sgp4!(orbp.sgp4d, 5*24*60)
0.000881 seconds (121 allocations: 9.029 KiB)
([-3149.8, 24100.6, 2890.71], [-2.61957, -0.196803, -0.04949])

julia> @time r,v = sgp4!(orbp.sgp4d, 5*24*60)
0.000014 seconds (28 allocations: 1.984 KiB)
([-3149.8, 24100.6, 2890.71], [-2.61957, -0.196803, -0.04949])
``````

# Matlab

``````l1 = '1 23599U 95029B   06171.76535463  .00085586  12891-6  12956-2 0  2905';
l2 = '2 23599   6.9327   0.2849 5782022 274.4436  25.2425  4.47796565123555';
[~, ~, ~, satrec] = twoline2rv(l1,l2,'c','m','a',84);
[satrec, r, v] = sgp4(satrec, 0);
tic; [satrec, r, v] = sgp4(satrec, 5*24*60); toc
Elapsed time is 0.000893 seconds.
tic; [satrec, r, v] = sgp4(satrec, 5*24*60); toc
Elapsed time is 0.000256 seconds.
>> r,v

r =

1.0e+04 *

-0.3150    2.4101    0.2891

v =

-2.6196   -0.1968   -0.0495
``````

Very nice. I have used the sgdp4 from Julia, but `ccall`ing the C version.

A detail comment about the Julia version in SatelliteToolbox.jl:

For a large nr of propagation steps the following would probably increase the speed: `propagator!` accepts a `Vector` of times `t`. Then the Arrays for the results can be preallocated avoiding the `push!` (similar as in Matlab extending Arrays beyond the last index at every step is not efficient). For a more Julian interface provide (two) `iterate` methods, (or `start`, `next`, and `done` in v0.6). For example, if one wants to propagate the orbits until satellites are within the reach of a position on the ground, and only then save the results, such an Iterator interface would be more natural.

The IGRF in SatelliteToolbox.jl should perhaps be moved eventually into a Geomag.jl package. There are also the magnetospheric models (e.g. Tsyganenko Magnetic Field Model and GEOPACK using internally the IGRF). And a number of high resolution models for low satellite orbits and the ground, e.g. here and here. They all use the method of spherical harmonics expansion. If at least some of them get converted to Julia, they would live conveniently in the same package.

Thanks @stephancb for the very good suggestions! I tested this when the deep space bits were not available. In my test, the overhead to create a new element in array using `push!` was negligible w.r.t. the computational burden of the propagator itself. However, I have changed many things (including some type-stability fixes) and I will check again! Thanks for this.

This seems very nice! Can you please tell me more about those methods? I have used `start`, `next`, and `done` with arrays and I am not entirely sure how can I adapt to an orbit propagator. For example, `start(orbp)` is easy because it will return the orbit elements, position and velocity at epoch. `next(orbp, t)` should give me the result at time `t`, but what about `done`? `done(orbp, t)` wouldn’t be exactly the same as `next`? I think I am a little confused here.

I agree! As soon as I (or someone) implements a new geomagnetic field algorithm, I will suggest to create another package!

About this, let me tell that the CM4 model is available in GMT and accessible via GMT.jl

2 Likes

Nice! Do you plan to implement CM4 natively on Julia? This can be a good point to create a separate project that will contain all related stuff to the Geomagnetic field.

My, how can I say, will to have everything written in Julia is because, IMO, 1 or 2 years from now it will be faster then calling a library and the code will get much more integrated and easy to maintain.

I am even thinking about ask for support to build a CubeSat to embedded a Julia compiler + code to handle attitude and orbit control (ok, let me dream a little right…)

Well, not really. I did the GMT’s CM4 integration via f2c + boring_cleaning_work and my geomag processing chain is still based in GMT. The CM4 code isn’t particularly fast so maybe it would get better if directly rewritten in Julia, but that’s quite some work. Besides, there is a new CM5 model but that is not yet free in the wild. I have contacted the authors but got no answer.

2 Likes

I have change to preallocate that array and the results were:

# Without preallocation

``````julia> @benchmark propagate!(\$orbp, collect(0:100:24*60)*60)
BenchmarkTools.Trial:
memory estimate:  31.97 KiB
allocs estimate:  376
--------------
minimum time:     20.032 μs (0.00% GC)
median time:      20.727 μs (0.00% GC)
mean time:        22.437 μs (5.83% GC)
maximum time:     1.079 ms (94.01% GC)
--------------
samples:          10000
evals/sample:     1
``````

# With preallocation

``````julia> @benchmark propagate!(\$orbp, collect(0:100:24*60)*60)
BenchmarkTools.Trial:
memory estimate:  31.55 KiB
allocs estimate:  367
--------------
minimum time:     19.422 μs (0.00% GC)
median time:      20.037 μs (0.00% GC)
mean time:        21.522 μs (5.13% GC)
maximum time:     898.383 μs (93.36% GC)
--------------
samples:          10000
evals/sample:     1
``````

Maybe I did something wrong. The new code is:

``````function propagate!(orbp::OrbitPropagatorSGP4{T}, t::Vector) where T
# Auxiliary variables.
orb   = orbp.orb
sgp4d = orbp.sgp4d

# Output.
result_orb = Array{Orbit}(length(t))
result_r   = Array{Vector}(length(t))
result_v   = Array{Vector}(length(t))

for i = 1:length(t)
k = t[i]
# Propagate the orbit.
(r_teme_k, v_teme_k) = sgp4!(sgp4d, k/60)

# Convert km to m.
r_teme_k *= 1000
v_teme_k *= 1000

# Update the elements in the `orb` structure.
@update_orb!(orbp, k)

result_orb[i] = copy(orb)
result_r[i]   = r_teme_k
result_v[i]   = v_teme_k
end

(result_orb, result_r, result_v)
end
``````

Any tips?

I ran this same benchmark using `SGP4.jl`, which is just a wrapper of the python library `sgp4`:

``````julia> l1 = "1 23599U 95029B   06171.76535463  .00085586  12891-6  12956-2 0  2905"
julia> l2 = "2 23599   6.9327   0.2849 5782022 274.4436  25.2425  4.47796565123555"
julia> wgs84 = SGP4.GravityModel("wgs84")
julia> satellite = twoline2rv(l1,l2,wgs84)
julia> @btime r,v = SGP4.sgp4(satellite, 5*24*60)
89.124 μs (89 allocations: 2.83 KiB)
([-3149.8, 24100.6, 2890.71], [-2.61957, -0.196803, -0.04949])
``````

compared to the native implementation:

``````julia> @btime r,v = sgp4!(orbp.sgp4d, 5*24*60)
1.308 μs (21 allocations: 1.75 KiB)
([-3149.8, 24100.6, 2890.71], [-2.61957, -0.196803, -0.04949])
``````

Not sure what to make of the small differences in the results, but the `SatelliteToolbox` implementation is obviously way faster.

Update: Edited the SGP4.jl benchmark to use a method with less overhead.

1 Like

This is awesome! Thanks for this. I don’t think that my implementation is good (I believe there is a lot improvements to do) but that most of this time in SGP4.jl is the overhead of calling the external library. I saw this before, and that is the reason why I impelemented the EGM96 and IGRF natively in Julia.

With respect to the differences, SatelliteToolbox.jl algorithm was based on the modifications described here https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev3.pdf

This TLE test was build precesily to catch a bug in deep space integration. However, I am pretty sure the lib already correct this. Hence, my bet is the algorithm to compute the Greenwich Mean Sederal time. Is it possible to propagate both to 500 min? You may have found a bug EDIT: The lib uses basically Vallado’s algorithm, which I use too. There is an option called `afspc_mode` or something like this to use a modern approach to compute GMST. In SatelliteToolbox case, this is always `true`. Now I am even more curious to see why we have this difference.

2 Likes

Your effort is really appreciated, but I do have a concern about the `propagate!`: It seems to be used as an integrator, i.e. the orbit state seems to be updated at each given time in the vector `t`. But an integrator is not provided by the original code. The original `sgp4` returns the position and velocity of the satellite at a single time `tsince`. My understanding is that the orbit state would be updated when new observations are available, as described in your reference. If SGP4 is used to update the orbit state based only on calculations, it may work in practice, but, as far as I can see, this would be untested/undefined, and should then perhaps be documented as such.

Now the for Julians more interesting part: if a Julia implementation of `sgp4(tle::OrbitState, tsince::Real)` returning position and velocity vectors is implemented (it may just call the C++ or Fortran code), then Julia will provide a dot version `sgp4.` that works for anything that can iterate over many `tsince`: if it has a `length(tsince)` (like a Vector), then (probably) `push!` will be avoided, if not `push!` will (probably) be used, but everything is just automatic, you don’t need to worry about it.

The following code snippet might list the times when the ISS would be above 60 deg elevation at some place on the ground, for tomorrow at 1 s resolution (which is needed, if you want to try to catch with a camera, v about 7.5.km/s).

``````for t in DateTime(2018,6,5):Dates.Second(1):DateTime(2018,6,6)
(r,v) = sgp4(ISSsatrec, t - DateTime(ISStle))
lc    = localcoordinates(r, t)
if lc.elevation>60*pi/180
@show t,lc
end
end
``````

(assuming that `localcoordinates(r, t)` returns local elevation and azimuth given a position `r` in ITRF and time `t`).

So I think that all you should provide is the basic `sgp4(satrec, tsince)`, as the original code, but returning `r` and `v`, i.e. not a `sgp4!`.

1 Like

Doesn’t `sgp4!` calculate the osculating element (`?_k`) from the mean elements and then stores them in the `SGP4_Structure`, while `propagate!` only updates the `Orbit` with current osculating elements for the step?

1 Like

Actually it is more or less what @benelsen said. When you call `sgp4!` with a specified `t`, it updates the osculating elements in `SGP4_Structure`. Notice that if you call it again with another `t`, it does not necessarily propagates everything from the epoch, because somethings (specially the deep space bits) can be computed recursively. This is precisely what the routine `sgp4` in Vallado’s code do.

The `propagate!` function is just an API that do this process for every instant in the input vector. So the only integrator is the one built inside the SGP4. Notice that all the tests in  are used to validade this algorithm using the `propagate!` function, which, under the hood, calls `sgp4!`.

I think I understood! I will study this a little bit because it seems very nice.

But `sgp4` is precisely the same algorithm of `sgp4(satrec,tsince)`. It returns the position and velocity. The problem is the `satrec` (in my case `sgp4d`) is modified inside the call to update the osculating elements and some other auxiliary variables. Hence, by Julia coding standards, it must have a `!`. The `propagate!` and `step!` API helps the user to select different propagators (Two Body, J2, SGP4, etc.) by changing just one argument in the function, like:

``````julia> orbp = init_orbit_propagator(Val{:J2}, Orbit(0.0,7130982.0,0.001111,98.405*pi/180,pi/2,0.0,0.0))
julia> (o,r,v) = propagate!(orbp, collect(0:3:24)*60*60)

julia> orbp = init_orbit_propagator(Val{:sgp4}, Orbit(0.0,7130982.0,0.001111,98.405*pi/180,pi/2,0.0,0.0))
julia> (o,r,v) = propagate!(orbp, collect(0:3:24)*60*60)
``````

Finally, thanks very much you all for all the contributions! I think we can do very good things here Obviously I had not checked the code very thoroughly. But then, to make use of the Julia dot functions, wouldn’t it be better to break up the SGP_structure, the mean elements are only output of sgp4? That is something like

`(r, v, m) = sgp4j(satrec, tsince)`.

The function names in the original code should be reserved for `ccall`ing them (which does not need to be implemented in the toolbox). Functions and structs written in pure Julia and possibly adjusted from the original to make things work the Julian way would be marked, for example by a `j` at the end?

I am not sure if I understood. There is not a SGP4 library in which you can `ccall` in SatelliteToolbox.jl. At least for SGP4, there will never be, since we now have a native implementation. Hence, I am not entirely sure why `sgp4!` should be renamed to `sgp4j`.