Getting Earth-Sun vectors using same ephemeris gives different results with different tools

[MichaelBonnet]

I have two functions that get TEME and GCRS Earth-Sun position vectors:

earth_sun_vectors.py

"""Utilities for getting Earth-Sun vectors in a reference frame at a UTC epoch."""

# Standard library imports
from contextlib import closing
from datetime import datetime, timezone

# Third part imports
import numpy as np
from skyfield.api import load
from skyfield.sgp4lib import TEME

def get_earth_sun_vector_gcrs_at_epoch(epoch: datetime) -> np.ndarray:
    """Get the Earth-Sun position vector in GCRS at given epoch.

    Args:
        epoch (datetime): The epoch to use for getting the sun vector.

    Returns:
        earth_sun_vector_gcrs_at_epoch (np.ndarray): the position of the sun in GCRS frame at epoch.

    """
    if not isinstance(epoch, datetime):
        raise TypeError(f"arg epoch must be of type datetime, not {str(type(epoch))}")
    ts = load.timescale()
    epoch = epoch.replace(tzinfo=timezone.utc)
    t = ts.from_datetime(epoch)
    bodies = load("de430_1850-2150.bsp")
    with closing(bodies):
        sun = bodies["sun"]
        earth = bodies["earth"]
        earth_sun_vector_gcrs_at_epoch = earth.at(t).observe(sun).apparent().position.km
        return np.array(earth_sun_vector_gcrs_at_epoch)

def get_earth_sun_vector_teme_at_epoch(epoch: datetime) -> np.ndarray:
    """Get the Earth-Sun position vector in TEME at given epoch.

    Args:
        epoch (datetime): The epoch to use for getting the sun vector.

    Returns:
        earth_sun_vector_teme_at_epoch (np.ndarray): the position of the sun in TEME frame at epoch.

    """
    if not isinstance(epoch, datetime):
        raise TypeError(f"arg epoch must be of type datetime, not {str(type(epoch))}")
    ts = load.timescale()
    epoch = epoch.replace(tzinfo=timezone.utc)
    t = ts.from_datetime(epoch)
    bodies = load("de430_1850-2150.bsp")
    with closing(bodies):
        sun = bodies["sun"]
        earth = bodies["earth"]
        apparent = earth.at(t).observe(sun).apparent()
        vector = apparent.frame_xyz(TEME).km
        return np.array(vector)

To validate the functionality, I have generated test data with a common industry tool, FreeFlyer. The FreeFlyer documentation clearly states that FreeFlyer will use the De430.dat file by default, hence why I use de430_1850-2150.bsp as my ephemeris file.

Below are my test functions, written using pytest:

test_earth_sun_vectors.py

"""Unit tests for Earth-Sun vector functions."""

# Standard library imports
from datetime import datetime

# Third party imports
import numpy as np
import pandas as pd
import pytest

# Local application imports
from custom_logger import setup_logging
from earth_sun_vectors import get_earth_sun_vector_gcrs_at_epoch, get_earth_sun_vector_teme_at_epoch

def unit_vector(x: np.ndarray) -> np.ndarray:
    """Returns the unit vector of the input.

    Args:
        x (np.ndarray): the input vector.

    Returns:
        x_unit (np.ndarray): the unit vector of the input vector.

    """
    return x / np.linalg.norm(x, axis=0)

def angle_between_vectors_deg(a: np.ndarray, b: np.ndarray) -> float:
    """Returns the angle between two vectors in R^3 in degrees

    Args:
        a (np.ndarray): the first vector.
        b (np.ndarray): the second vector.

    Returns:
        angle_deg (float): the angle between the two vectors in degrees.

    """
    return np.rad2deg(np.arccos(np.sum(unit_vector(a) * unit_vector(b), axis=0)))

# get_earth_sun_vector_gcrs_at_epoch()
def test_get_earth_sun_vector_gcrs_at_epoch_with_freeflyer_data():

    # logging
    logger = setup_logging(f"get_earth_sun_vector_gcrs_at_epoch_{datetime.now()}")

    test_data = pd.read_csv(
        "earth_sun_vector_test_data.csv"
    )

    test_count = 0

    # Validate the results
    for _, row in test_data.iterrows():
        test_count += 1

        # extract arg for function
        epoch = datetime.fromisoformat(row["Epoch_ISO8601_UTC"][:26])

        # get actual sun vector
        test_position_vector = get_earth_sun_vector_gcrs_at_epoch(epoch)

        # get desired sun vector
        true_position_vector = np.array(
            [
                row["Sun_Earth_MJ2000_X_km"],  # MJ2000 is supposedly accurate to GCRS within milliarcseconds
                row["Sun_Earth_MJ2000_Y_km"],
                row["Sun_Earth_MJ2000_Z_km"],
            ]
        )
        true_position_vector = -true_position_vector  # validation data has Sun-Earth, we want Earth-Sun

        angle_between = angle_between_vectors_deg(test_position_vector, true_position_vector)

        logger.info(f"Test Case: {test_count}\tCalculated Vector: {test_position_vector}\tTrue Vector: {true_position_vector}\tAngle Diff (degs): {angle_between}\tAngle Diff (arcsec): {angle_between * 3600}")

        # accurate to 1/10th of 1%
        np.testing.assert_allclose(test_position_vector, true_position_vector, rtol=0.001)
        assert np.abs(angle_between) < 0.01

# get_earth_sun_vector_teme_at_epoch()
def test_get_earth_sun_vector_teme_at_epoch_with_freeflyer_data():

    # logging
    logger = setup_logging(f"get_earth_sun_vector_teme_at_epoch_{datetime.now()}")

    test_data = pd.read_csv(
        "earth_sun_vector_test_data.csv"
    )

    test_count = 0

    # Validate the results
    for _, row in test_data.iterrows():
        test_count += 1

        # extract arg for function
        epoch = datetime.fromisoformat(row["Epoch_ISO8601_UTC"][:26])

        # get actual sun vector
        test_position_vector = get_earth_sun_vector_teme_at_epoch(epoch)

        # get desired sun vector
        true_position_vector = np.array(
            [
                row["Sun_Earth_TEME_X_km"],
                row["Sun_Earth_TEME_Y_km"],
                row["Sun_Earth_TEME_Z_km"],
            ]
        )
        true_position_vector = -true_position_vector  # validation data has Sun-Earth, we want Earth-Sun

        angle_between = angle_between_vectors_deg(test_position_vector, true_position_vector)

        logger.info(f"Test Case: {test_count}\tCalculated Vector: {test_position_vector}\tTrue Vector: {true_position_vector}\tAngle Diff (degs): {angle_between}\tAngle Diff (arcsec): {angle_between * 3600}")

        # accurate to one-tenth of 1% on each element
        np.testing.assert_allclose(test_position_vector, true_position_vector, rtol=0.001)
        # accurate to about ~0.00578 degrees, or about 21 arcseconds
        assert np.abs(angle_between) < 0.01

I have created a gist with these and other files, including the test data csv and Docker-related files to minimize possible differences between machines.

An rtol of 0.001 is only 1/10th of 1% which seems not that accurate, especially given that they're seemingly using the same ephemeris. More concerning is the variance in the 1e1 arcseconds range, which could throw off things like eclipse calculations.

My chief concern is understanding why my methodology is resulting in different results from my test data, when both my methods and the test data generating tool are using the same ephemeris. What might be the root cause of this difference? Am I using an incorrect function?

Thanks for any and all help.

[MichaelBonnet]

For what it's worth, I switched to using a better angle calculation function that gets around potential precision issues inherent to arccos:

def angle_between_vectors_deg(a: np.ndarray, b: np.ndarray) -> float:
    n_a = np.linalg.norm(a, axis=0)
    n_b = np.linalg.norm(b, axis=0)
    y = np.linalg.norm(n_b * a - n_a * b, axis=0)
    x = np.linalg.norm(n_b * a + n_a * b, axis=0)
    return np.rad2deg(2 * np.arctan2(y, x))

but my issue still remains.

[brandon-rhodes]

I've answered more at length over on your Stack Exchange question:

https://astronomy.stackexchange.com/questions/55633/resolving-variance-between-skyfield-generated-sun-vectors-and-test-data-when-bot/55677#55677

But, in short, you are applying the affects of aberration and deflection when you call .apparent(), but FreeFlyer doesn't appear to do those calculations, so if you leave that call off you should get much closer agreement. I'm going to close this issue, as I think that answers your question, but feel free to comment further about any lingering questions.