Julia now has a native SGP4/SDP4 implementation

Specific Domains Astro/Space
1

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 [1] 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 :smiley:

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 :slight_smile:

[1]: Vallado, D. A., Crawford, P., Hujsak, R., Kelso, T. S (2006). Revisiting Spacetrack Report #3: Rev1. AIAA.

12 Likes
2

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
3

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
4

Count me in :+1:

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

2 Likes
5

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
6

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
7

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[1])
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
8

Very nice. I have used the sgdp4 from Julia, but ccalling 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.

9

Thanks @stephancb for the very good suggestions! :slight_smile:

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!

10

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

2 Likes
11

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 :smiley: 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 :smiley: (ok, let me dream a little right…)

12

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
13

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?

14

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
15

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 :smile:

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
16

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
17

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
18

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 [1] 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 :smiley:

19

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 ccalling 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?

20

Hi @stephancb

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.