# Julia now has a native SGP4/SDP4 implementation

**URL:** <https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376>\
**Category:** Astro/Space\
**Created:** [June 2, 2018, 9:13pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376 "2018-06-02T21:13:23Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 2, 2018, 9:13pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/1 "2018-06-02T21:13:23Z")

</div>

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](https://github.com/SatelliteToolbox/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 😃

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:

```fortran
    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

```julia
    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 🙂

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

---

<div class="post-metadata">

**Author:** ![crbinz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/crbinz/32/1919_2.png) [@crbinz](https://discourse.julialang.org/u/crbinz)\
**Post date:** [June 3, 2018, 12:12am UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/2 "2018-06-03T00:12:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 12:37am UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/3 "2018-06-03T00:37:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![helgee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgee/32/2022_2.png) [@helgee](https://discourse.julialang.org/u/helgee)\
**Post date:** [June 3, 2018, 11:49am UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/4 "2018-06-03T11:49:40Z")

</div>

> [@Ronis\_BR](#):
>
> help needed!

Count me in 👍

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 3, 2018, 12:07pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/5 "2018-06-03T12:07:47Z")

</div>

> [@Ronis\_BR](#):
>
> the first version of the SGP4/SDP4 implementation 100% written in Julia

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).

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 1:20pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/6 "2018-06-03T13:20:48Z")

</div>

> [@Tamas\_Papp](#):
>
> 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).

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.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 1:38pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/7 "2018-06-03T13:38:17Z")

</div>

Well, at least against Matlab it is pretty fast:

# Julia

```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

```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

```

---

<div class="post-metadata">

**Author:** ![stephancb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephancb/32/14243_2.png) [@stephancb](https://discourse.julialang.org/u/stephancb)\
**Post date:** [June 3, 2018, 8:13pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/8 "2018-06-03T20:13:29Z")

</div>

Very nice. I have used the sgdp4 from Julia, but `ccall`ing [the C version](https://github.com/cbassa/strf).

A detail comment about the Julia version in [SatelliteToolbox.jl](https://github.com/SatelliteToolbox/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](https://github.com/SatelliteToolbox/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](https://ccmc.gsfc.nasa.gov/modelweb/magnetos/tsygan.html) using internally the IGRF). And a number of high resolution models for low satellite orbits and the ground, e.g. [here](https://geomag.us/models/index.html) and [here](http://www.space.dtu.dk/english/Research/Scientific_data_and_models/Magnetic_Field_Models). 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.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 9:43pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/9 "2018-06-03T21:43:57Z")

</div>

Thanks @stephancb for the very good suggestions! 🙂

> [@stephancb](#):
>
> 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).

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.

> [@stephancb](#):
>
> 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.

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.

> [@stephancb](#):
>
> The IGRF in [SatelliteToolbox.jl](https://github.com/SatelliteToolbox/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](https://ccmc.gsfc.nasa.gov/modelweb/magnetos/tsygan.html) using internally the IGRF). And a number of high resolution models for low satellite orbits and the ground, e.g. [here](https://geomag.us/models/index.html) and [here](http://www.space.dtu.dk/english/Research/Scientific_data_and_models/Magnetic_Field_Models). 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.

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

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 3, 2018, 9:56pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/10 "2018-06-03T21:56:29Z")

</div>

> [@Ronis\_BR](#):
>
> 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](https://gmt.soest.hawaii.edu/doc/5.4.3/supplements/mgd77/mgd77magref.html) and accessible via [GMT.jl](https://github.com/GenericMappingTools/GMT.jl)

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 10:23pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/11 "2018-06-03T22:23:27Z")

</div>

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…)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 3, 2018, 10:32pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/12 "2018-06-03T22:32:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 3, 2018, 10:43pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/13 "2018-06-03T22:43:48Z")

</div>

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

# Without preallocation

```julia
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
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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![crbinz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/crbinz/32/1919_2.png) [@crbinz](https://discourse.julialang.org/u/crbinz)\
**Post date:** [June 4, 2018, 1:53am UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/14 "2018-06-04T01:53:29Z")

</div>

I ran this same benchmark using [`SGP4.jl`](https://github.com/crbinz/SGP4.jl), which is just a wrapper of the python library [`sgp4`](https://github.com/brandon-rhodes/python-sgp4):

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 4, 2018, 2:45am UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/15 "2018-06-04T02:45:51Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![stephancb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephancb/32/14243_2.png) [@stephancb](https://discourse.julialang.org/u/stephancb)\
**Post date:** [June 4, 2018, 1:43pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/16 "2018-06-04T13:43:38Z")

</div>

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).

```julia
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!`.

---

<div class="post-metadata">

**Author:** ![benelsen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/benelsen/32/1620_2.png) [@benelsen](https://discourse.julialang.org/u/benelsen)\
**Post date:** [June 4, 2018, 3:44pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/17 "2018-06-04T15:44:51Z")

</div>

> [@stephancb](#):
>
> 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.

Doesn’t `sgp4!` calculate the osculating element (`?_k`) [from the mean elements](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4.jl#L373) and then [stores them](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4.jl#L618) in the `SGP4_Structure`, while `propagate!` only [updates](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4_api.jl#L54) the `Orbit` with current osculating elements for the step?

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 4, 2018, 4:41pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/18 "2018-06-04T16:41:17Z")

</div>

> [@benelsen](#):
>
> Doesn’t `sgp4!` calculate the osculating element ( `?_k` ) [from the mean elements](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4.jl#L373) and then [stores them](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4.jl#L618) in the `SGP4_Structure` , while `propagate!` only [updates](https://github.com/SatelliteToolbox/SatelliteToolbox.jl/blob/master/src/orbit/propagators/sgp4_api.jl#L54) the `Orbit` with current osculating elements for the step?

> [@stephancb](#):
>
> 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.

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!`.

> [@stephancb](#):
>
> 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).
> 
> ```julia
> 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` ).

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

> [@stephancb](#):
>
> 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!` .

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
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 😃

---

<div class="post-metadata">

**Author:** ![stephancb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephancb/32/14243_2.png) [@stephancb](https://discourse.julialang.org/u/stephancb)\
**Post date:** [June 4, 2018, 5:04pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/19 "2018-06-04T17:04:56Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 4, 2018, 9:04pm UTC](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376/20 "2018-06-04T21:04:41Z")

</div>

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`.

[Next page](https://discourse.julialang.org/t/julia-now-has-a-native-sgp4-sdp4-implementation/11376.md?page=2)