Skyfield coordinates differ from r_BN_N

[Kreolis]

Hi,
I am running my simulation to compare it to some of my own implementations of sgp4 etc.
I have tried to verify the output of sgp4 to skyfield and it aligns well. However when I plug the same TLE and also start time to basilisk the output from r_BN_N is significantly different. I have set the inital position like so:

import skyfield.api as skyfield
ts = skyfield.load.timescale()
satellite = skyfield.EarthSatellite(self.tle1, self.tle2)

# Parse the string into a datetime object
parsed_datetime = datetime.strptime(self.date_and_time, '%Y %B %d %H:%M:%S.%f (UTC)')

# Extract the year, month, day, hour, and minute
year = parsed_datetime.year
month = parsed_datetime.month
day = parsed_datetime.day
hour = parsed_datetime.hour
minute = parsed_datetime.minute
second = parsed_datetime.second

# Parse UTC datetime string
t = ts.utc(year, month, day, hour, minute, second)
geocentric = satellite.at(t)


# Note that these vectors are given here relative to the Earth frame (zeroBase = 'Earth').  When we set the spacecraft()
# initial position and velocity vectors through before initialization
self.spacecraft_object.hub.r_CN_NInit = geocentric.position.m  # Initial position

Also note that zeroBase = 'Earth' is set.
So what am I doing wrong here? Is the ouput of basilisk not in ECI?

[Kreolis]

This should be simple... But also here it drifts quite much after only a couple of orbits.

import os
from copy import copy

import matplotlib.pyplot as plt
import numpy as np

import skyfield.api as skyfield
from skyfield.framelib import galactic_frame
from datetime import datetime
ts = skyfield.load.timescale()

# TLE from SUOMI-100
tle1 = "1 43804U 18099AY  25324.11800123  .00012160  00000-0  52655-3 0  9994"
tle2 = "2 43804  97.4522  20.7660 0009972 300.6855  59.3398 15.22511975381541"
satellite = skyfield.EarthSatellite(tle1, tle2)

date_and_time = '2025 November 20 11:13:0.0 (UTC)'

# Parse the string into a datetime object
parsed_datetime = datetime.strptime(date_and_time, '%Y %B %d %H:%M:%S.%f (UTC)')

# Extract the year, month, day, hour, and minute
year = parsed_datetime.year
month = parsed_datetime.month
day = parsed_datetime.day
hour = parsed_datetime.hour
minute = parsed_datetime.minute
second = parsed_datetime.second

from Basilisk import __path__

bskPath = __path__[0]
fileName = os.path.basename(os.path.splitext(__file__)[0])

from Basilisk.simulation import spacecraft

from Basilisk.utilities import (
    SimulationBaseClass,
    macros,
    orbitalMotion,
    simIncludeGravBody,
    unitTestSupport,
    vizSupport,
)


from Basilisk.utilities.supportDataTools.dataFetcher import get_path, DataFile


def run(show_plots, orbitCase, useSphericalHarmonics, planetCase):

    simTaskName = "simTask"
    simProcessName = "simProcess"


    scSim = SimulationBaseClass.SimBaseClass()

    scSim.SetProgressBar(True)


    dynProcess = scSim.CreateNewProcess(simProcessName)


    simulationTimeStep = macros.sec2nano(10.0)
    dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))

    scObject = spacecraft.Spacecraft()
    scObject.ModelTag = "bsk-Sat"


    scSim.AddModelToTask(simTaskName, scObject)

    gravFactory = simIncludeGravBody.gravBodyFactory()


    planet = gravFactory.createEarth()
    planet.isCentralBody = True  # ensure this is the central gravitational body
    if useSphericalHarmonics:
        # If extra customization is required, see the createEarth() macro to change additional values.
        # For example, the spherical harmonics are turned off by default.  To engage them, the following code
        # is used
        ggm03s_j2_only_path = get_path(DataFile.LocalGravData.GGM03S_J2_only)
        planet.useSphericalHarmonicsGravityModel(str(ggm03s_j2_only_path), 2)

        # The value 2 indicates that the first two harmonics, excluding the 0th order harmonic,
        # are included.  This harmonics data file only includes a zeroth order and J2 term.
    mu = planet.mu


    gravFactory.addBodiesTo(scObject)
    gravFactory.createSpiceInterface(time=date_and_time, epochInMsg=True)


    # Parse UTC datetime string
    t = ts.utc(year, month, day, hour, minute, second)
    geocentric = satellite.at(t)

    # To set the spacecraft initial conditions, the following initial position and velocity variables are set:
    scObject.hub.r_CN_NInit = geocentric.position.m # m   - r_BN_N
    scObject.hub.v_CN_NInit = geocentric.velocity.m_per_s # m/s - v_BN_N

    n = np.sqrt(mu / oe.a / oe.a / oe.a)
    P = 2.0 * np.pi / n
    if useSphericalHarmonics:
        simulationTime = macros.sec2nano(50.0 * P)
    else:
        simulationTime = macros.sec2nano(50.0 * P)

    # Setup data logging before the simulation is initialized
    if useSphericalHarmonics:
        numDataPoints = 400
    else:
        numDataPoints = 100
    samplingTime = unitTestSupport.samplingTime(
        simulationTime, simulationTimeStep, numDataPoints
    )

    dataRec = scObject.scStateOutMsg.recorder(samplingTime)
    scSim.AddModelToTask(simTaskName, dataRec)

    if vizSupport.vizFound:
        viz = vizSupport.enableUnityVisualization(scSim, simTaskName, scObject,
                                                  # saveFile=__file__
                                                  # liveStream=True
                                                  )

   
    scSim.InitializeSimulation()


    scSim.ConfigureStopTime(simulationTime)
    scSim.ExecuteSimulation()
 
    posData = dataRec.r_BN_N
    velData = dataRec.v_BN_N

    np.set_printoptions(precision=16)

    figureList, finalDiff = plotOrbits(
        dataRec.times(),
        posData,
        velData,
        oe,
        mu,
        P,
        orbitCase,
        useSphericalHarmonics,
        planetCase,
        planet,
    )

    if show_plots:
        plt.show()


    plt.close("all")

    return finalDiff, figureList


def plotOrbits(
    timeAxis,
    posData,
    velData,
    oe,
    mu,
    P,
    orbitCase,
    useSphericalHarmonics,
    planetCase,
    planet,
):

    plt.close("all")  # clears out plots from earlier test runs
    plt.figure(1)
    fig = plt.gcf()
    ax = fig.gca()
    ax.ticklabel_format(useOffset=False, style="plain")
    finalDiff = 0.0

    for idx in range(3):
        plt.plot(
            timeAxis * macros.NANO2SEC / P,
            posData[:, idx] / 1000.0,
            color=unitTestSupport.getLineColor(idx, 3),
            label="$r_{BN," + str(idx) + "}$",
        )
    plt.legend(loc="lower right")
    plt.xlabel("Time [orbits]")
    plt.ylabel("Inertial Position [km]")
    figureList = {}
    pltName = fileName + "1" + orbitCase + str(int(useSphericalHarmonics)) + planetCase
    figureList[pltName] = plt.figure(1)

    if useSphericalHarmonics is False:
        # draw orbit in perifocal frame
        b = oe.a * np.sqrt(1 - oe.e * oe.e)
        p = oe.a * (1 - oe.e * oe.e)
        plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
        plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
        # draw the planet
        fig = plt.gcf()
        ax = fig.gca()
        if planetCase == "Mars":
            planetColor = "#884400"
        else:
            planetColor = "#008800"
        planetRadius = planet.radEquator / 1000
        ax.add_artist(plt.Circle((0, 0), planetRadius, color=planetColor))

        rData = []
        fData = []
        for idx in range(0, len(posData)):
            oeData = orbitalMotion.rv2elem(mu, posData[idx], velData[idx])
            rData.append(oeData.rmag)
            fData.append(oeData.f + oeData.omega - oe.omega)
        plt.plot(
            rData * np.cos(fData) / 1000,
            rData * np.sin(fData) / 1000,
            color="#aa0000",
            linewidth=3.0,
        )

        fData = np.linspace(0, 2 * np.pi, 100)
        rData = []
        for idx in range(0, len(fData)):
            rData.append(p / (1 + oe.e * np.cos(fData[idx])))
        plt.plot(
            rData * np.cos(fData) / 1000,
            rData * np.sin(fData) / 1000,
            "--",
            color="#555555",
        )
        plt.xlabel("$i_e$ Cord. [km]")
        plt.ylabel("$i_p$ Cord. [km]")
        plt.grid()

        plt.figure(3)
        fig = plt.gcf()
        ax = fig.gca()
        ax.ticklabel_format(useOffset=False, style="plain")
        Deltar = np.empty((0, 3))
        E0 = orbitalMotion.f2E(oe.f, oe.e)
        M0 = orbitalMotion.E2M(E0, oe.e)
        n = np.sqrt(mu / (oe.a * oe.a * oe.a))
        oe2 = copy(oe)
        for idx in range(0, len(posData)):
            M = M0 + n * timeAxis[idx] * macros.NANO2SEC
            Et = orbitalMotion.M2E(M, oe.e)
            oe2.f = orbitalMotion.E2f(Et, oe.e)
            rv, vv = orbitalMotion.elem2rv(mu, oe2)
            Deltar = np.append(Deltar, [posData[idx] - rv], axis=0)

        skyfield_pos = []
        for idx in range(0, len(timeAxis)):
            t = ts.utc(year, month, day, hour, minute, second + timeAxis[idx] * macros.NANO2SEC)
            geocentric = satellite.at(t)
            skyfield_pos.append(geocentric.position.m)
        skyfield_pos = np.array(skyfield_pos)
        DeltarSkyfield = posData - skyfield_pos

        for idx in range(3):
            plt.plot(
                timeAxis * macros.NANO2SEC / P,
                DeltarSkyfield[:, idx],
                "--",
                color=unitTestSupport.getLineColor(idx, 3),
                label=r"$\Delta r_{Skyfield," + str(idx) + "}$",
            )
        plt.legend(loc="lower right")
        plt.xlabel("Time [orbits]")
        plt.ylabel("Trajectory Differences [m]")
        pltName = (
            fileName + "3" + orbitCase + str(int(useSphericalHarmonics)) + planetCase
        )
        figureList[pltName] = plt.figure(3)

        finalDiff = np.linalg.norm(Deltar[-1])

    else:
        plt.figure(2)
        fig = plt.gcf()
        ax = fig.gca()
        ax.ticklabel_format(useOffset=False, style="plain")
        smaData = []
        for idx in range(0, len(posData)):
            oeData = orbitalMotion.rv2elem(mu, posData[idx], velData[idx])
            smaData.append(oeData.a / 1000.0)
        plt.plot(
            timeAxis * macros.NANO2SEC / P,
            smaData,
            color="#aa0000",
        )
        plt.xlabel("Time [orbits]")
        plt.ylabel("SMA [km]")

        plt.figure(3)
        fig = plt.gcf()
        ax = fig.gca()
        ax.ticklabel_format(useOffset=False, style="plain")
        Deltar = np.empty((0, 3))
        E0 = orbitalMotion.f2E(oe.f, oe.e)
        M0 = orbitalMotion.E2M(E0, oe.e)
        n = np.sqrt(mu / (oe.a * oe.a * oe.a))
        oe2 = copy(oe)
        for idx in range(0, len(posData)):
            M = M0 + n * timeAxis[idx] * macros.NANO2SEC
            Et = orbitalMotion.M2E(M, oe.e)
            oe2.f = orbitalMotion.E2f(Et, oe.e)
            rv, vv = orbitalMotion.elem2rv(mu, oe2)
            Deltar = np.append(Deltar, [posData[idx] - rv], axis=0)

        skyfield_pos = []
        for idx in range(0, len(timeAxis)):
            t = ts.utc(year, month, day, hour, minute, second + timeAxis[idx] * macros.NANO2SEC)
            geocentric = satellite.at(t)
            skyfield_pos.append(geocentric.position.m)
        skyfield_pos = np.array(skyfield_pos)
        DeltarSkyfield = posData - skyfield_pos

        for idx in range(3):
            plt.plot(
                timeAxis * macros.NANO2SEC / P,
                DeltarSkyfield[:, idx],
                "--",
                color=unitTestSupport.getLineColor(idx, 3),
                label=r"$\Delta r_{Skyfield," + str(idx) + "}$",
            )
        plt.legend(loc="lower right")
        plt.xlabel("Time [orbits]")
        plt.ylabel("Trajectory Differences [m]")
        pltName = (
            fileName + "3" + orbitCase + str(int(useSphericalHarmonics)) + planetCase
        )
        figureList[pltName] = plt.figure(3)

    pltName = fileName + "2" + orbitCase + str(int(useSphericalHarmonics)) + planetCase
    figureList[pltName] = plt.figure(2)
    return figureList, finalDiff


if __name__ == "__main__":
    run(
        True,  # show_plots
        "LEO",  # orbit Case (LEO, GTO, GEO)
        True,  # useSphericalHarmonics
        "Earth",  # planetCase (Earth, Mars)
    )

[Kreolis]

fixed by moving the initialization of spacecraft init behind environmental (gravfactory) initialization