Implement the True Equator Mean Equinox (TEME) coordinate system for more accurate ra/dec calculations

[szabolcsvelkei]

Is your feature request related to a problem? Please describe. I found that with the existing frames there is a residual error when calculating the ra/dec position of known satellites based on TLE.

Additional context https://celestrak.org/publications/AIAA/2008-6770/AIAA-2008-6770.pdf

A good solution example with Skyfield:

from skyfield.api import load
from skyfield.api import EarthSatellite
from skyfield.api import wgs84
ts = load.timescale()
line1 = '1 39348U 10057N   24238.91466777  .00000306  00000-0  19116-2 0  9995'
line2 = '2 39348  20.0230 212.2863 7218258 312.9449   5.6833  2.25781763 89468'
satellite = EarthSatellite(line1, line2, 'CZ-3C DEB', ts)
t = ts.utc(2024, 8, 26, 22, 34, 20)
k88 = wgs84.latlon(+47.91748, 19.89367, 984)
difference = satellite - k88
topocentric = difference.at(t)
print(topocentric.position.km)
ra, dec, distance = topocentric.radec()  # ICRF ("J2000")
print(t)
print(ra)
print(dec)

Skyfield results:

<Time tt=2460549.44131>
22h 06m 21.76s
+11deg 51' 32.3"

And the values ​​obtained based on our telescope observation:

22h 06m 23.67s
+11deg 50' 52.9"

Beware of the time transformation! 2460549.44131 = 2024.08.26 22:35:29.2 UTC = 2024.08.26 22:34:20 +69.2 sec! ( dT: http://ytliu.epizy.com/eclipse/dynamical_time.html?i=1 ) You have to use the tdt transformation to calculate the positions for the right epoch var epoch = new Time(t, TimeFrame.UTCFrame).ToTDT();

[sylvain-guillet]

Hello @szabolcsvelkei,

Yesterday I released a new preview of the framework that support TEME frame.

You can reproduce your scenario with this code :

TLE tle = TLE.Create("CZ-3C DEB", "1 39348U 10057N   24238.91466777  .00000306  00000-0  19116-2 0  9995", "2 39348  20.0230 212.2863 7218258 312.9449   5.6833  2.25781763 89468");

TimeSystem.Time epoch = new TimeSystem.Time("2024-08-26T22:34:20.00000Z");

var stateVector = tle.ToStateVector(epoch);

Site site = new Site(14, "SiteA", TestHelpers.EarthAtJ2000, new Planetodetic(19.89367 * Astrodynamics.Constants.Deg2Rad, 47.91748 * Astrodynamics.Constants.Deg2Rad, 984));

var eq = stateVector.RelativeTo(site, Aberration.None).ToEquatorial();

double ra = eq.RightAscension * Astrodynamics.Constants.Rad2Deg; //So, 331.591° is approximately 22 hours, 6 minutes, and 21.93 seconds
double dec = eq.Declination * Astrodynamics.Constants.Rad2Deg; //So, 11.859° is approximately 11°51'32.4"

You can update your project with this command : dotnet add package IO.Astrodynamics --version 6.0.0-preview-2

I hope this will be useful, don't hesitate to give me feedback.

Sylvain

[szabolcsvelkei]

Wonderful. Thank you, and sorry for the late reply