Why the track is not drifting?

[ZMAlt]

Hi. I used an example from the SatelliteToolbox.jl‘s document.

orbp = init_orbit_propagator(Val(:twobody), 
    KeplerianElements(0.0, 7130982.0, 0.001111, 98.405*pi/180, pi/2, 0.0, 0.0) )
r,_ = propagate!(orbp, collect(0 : 5 : 5 * 24 * 3600))

I convert r to geodetic coordinates. And the coordinates are drawn using GMT. Why the track is not drifting?

[helgee]

Because you are using the Two Body propagator which assumes a spherical Earth. The precession of the orbital plane is caused by Earth's oblateness which means you need to consider the J2 term and use the J2 propagator.

E.g.

orbp = init_orbit_propagator(Val(:J2), 
    KeplerianElements(0.0, 7130982.0, 0.001111, 98.405*pi/180, pi/2, 0.0, 0.0) )
r,_ = propagate!(orbp, collect(0 : 5 : 5 * 24 * 3600))

[ZMAlt]

Hi, @helgee I have modified the code as you suggested. And I got this.

However, I want to get the trajectory in the figure below. I used the same parameters to get this on HOMA.

[helgee]

Can you post your whole code or an MWE (minimum working example) so that we have an easier time debugging?

Do you know what kind of force model HOMA uses? I could not find information about that on the website.

[ZMAlt]

Can you post your whole code or an MWE (minimum working example) so that we have an easier time debugging?

Julia code

using DelimitedFiles
using SatelliteToolbox
orbp = init_orbit_propagator(Val(:J2), 
    KeplerianElements(0.0, 7130982.0, 0.001111, 98.405*pi/180, pi/2, 0.0, 0.0) )
r,_ = propagate!(orbp, collect(0 : 5 : 5 * 24 * 3600))

lat = 0.0; lon = 0.0
r_ell = zeros(length(r), 2)
@inbounds for (i,pnt) in enumerate(r)
    lat, lon = ecef_to_geodetic(pnt)
    r_ell[i,:] = [rad2deg(lon),rad2deg(lat)]
end
writedlm("nadir.txt",r_ell)

GMT script

gmt begin map pdf
    gmt coast -G220 -Rg -JQ8i -B
    gmt plot -W0.1p nadir.txt
gmt end show

---

Do you know what kind of force model HOMA uses? I could not find information about that on the website.

I have no idea.

[ronisbr]

Hi @ZMAlt !

I think the problem is because propagate! will return the state vector r, v represented in the same reference frame of the orbit elements. In this example (Amazonia 1 orbit), it is represented in the True-of-Date, which is an ECI frame. Before converting to latitude and longitude, you need to convert to ECEF. This transformation can be done by:

r_ecef = r_eci_to_ecef(TOD(), PEF(), jd) * r

In which jd is the Julian day at each propagation step. In your case, the Keplerian elements is being initialized with 0. I suggest adding some random day so that you can use to obtain this transformation matrix.

[ronisbr]

This seems to be working:

using DelimitedFiles
using SatelliteToolbox

jd0 = date_to_jd(2020, 1, 1, 0, 0, 0)
orbp = init_orbit_propagator(
    Val(:J2),
    KeplerianElements(jd0, 7130982.0, 0.001111, 98.405*pi/180, pi/2, 0.0, 0.0)
)

t = collect(0 : 5 : 5 * 24 * 3600)
r_tod, _ = propagate!(orbp, t)
r_ecef = map((t, r)->r_eci_to_ecef(TOD(), PEF(), jd0 + t / 86400) * r, t, r_tod)

lat = 0.0; lon = 0.0
r_ell = zeros(length(r_tod), 2)
@inbounds for (i,pnt) in enumerate(r_ecef)
    lat, lon = ecef_to_geodetic(pnt)
    r_ell[i,:] = [rad2deg(lon),rad2deg(lat)]
end

[ZMAlt]

It works! ! Thank you. @ronisbr @helgee