search

JuliaSpace / SatelliteToolboxCelestialBodies.jl

General functions related to celestial bodies for the SatelliteToolbox.jl ecosystem.
MIT License
2 stars 2 forks source link

Compared to the position of the sun in STK, there is a significant difference. Why does this situation occur? #3

Open blackbutton opened 2 months ago

blackbutton commented 2 months ago 
 jd = date_to_jd(2024,8,6,4,0,0)
s_mod = sun_position_mod(jd)
r_mod_j2000 = r_eci_to_eci(MOD(),J2000(),jd)
r_mod_j2000 * s_mod
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -1.0512053177080336e11
  1.0039274359027356e11
  4.351950513418951e10

STK: -105129703.817743 100395407.099930 43520964.349647

ronisbr commented 2 months ago

Hi @blackbutton !

The algorithm we currently have is a "low-precision" implemented from the Astronomical Almanac (and described on Vallado's book). We still lack a high precision Sun position algorithm. Nevertheless, there is a minor problem in your code because sun_position_mod expects the time in TBD (Barycentric Dynamical Time) instead of UTC. The difference between the Terrestrial time and TBD is period by Earth's orbit. Hence, we can assume that both are the same here:

julia> jd = date_to_jd(2024, 8, 6, 4, 0, 0)
2.4605286666666665e6

julia> jd_utc = date_to_jd(2024, 8, 6, 4, 0, 0)
2.4605286666666665e6

julia> jd_tt = jd_utc_to_tt(jd_utc)
2.4605286674674074e6

julia> s_mod = sun_position_mod(jd_tt)
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -1.05776297878576e11
  9.981141741230861e10
  4.326692615766075e10

julia> s_j2000 = r_eci_to_eci(MOD(), J2000(), jd) * s_mod
3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -1.0512198446655511e11
  1.003914406336404e11
  4.351894032716577e10

Now, if we take the angular difference between the STK answer and the one we just obtained, we get:

julia> s_j2000_stk = [-105129703.817743, 100395407.099930, 43520964.349647]
3-element Vector{Float64}:
 -1.05129703817743e8
  1.0039540709993e8
  4.3520964349647e7

julia> using LinearAlgebra

julia> acosd(s_j2000 / norm(s_j2000) ⋅ s_j2000_stk / norm(s_j2000_stk)) * 3600
3.4026739070005374

The result indicates that we differ from STK by less than 3.5 arcs, which is more that enough for most satellite applications. If you want a really accurate algorithm for satellite purposes, we can implement it here.