4
$\begingroup$

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

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.

$\endgroup$
7
  • 2
    $\begingroup$ Can't really tell from the code, but a common mistake is passing UTC to the ephemeris function expecting barycentric dynamic time. Some tools will convert automatically, some don't. $\endgroup$ Dec 30, 2023 at 14:49
  • 1
    $\begingroup$ @GregMiller thankfully, the ts.from_datetime(datetime_object) is designed to ingest a datetime object of any timezone and automatically makes adjustments for it, as detailed in the documentation $\endgroup$ Jan 2 at 22:27
  • 1
    $\begingroup$ Two additional facts would help answer your question: first, could you identify the exact date and time of the greatest error measured between your two data sources? That would let users focus on the exact vectors for that moment in time, instead of the whole slate of test vectors. Second: once you have identified that moment, you could check the GCRS vectors at least against NASA HORIZONS, so you would have a "tiebreaker vote" as to whether we're looking at a Skyfield or a FreeFlyer problem, or both. $\endgroup$ Jan 4 at 4:45
  • 1
    $\begingroup$ @BrandonRhodes the only roadblock is what type and format of ephemeris to generate with NASA HORIZONS, given that you can't get SPK ephemeris files for the Sun, and I don't know if Skyfield can ingest ephemeris from the .txt format that the web tool uses. I am of course open to guidance. $\endgroup$ Jan 4 at 22:46
  • 1
    $\begingroup$ @BrandonRhodes the exact date and time of greatest error measured between the two data sources when comparing the GCRS method is 2023-12-06T14:59:30.00000. Interestingly, this is at the maximum index of the test data. The minimum index of the test data has the least error. $\endgroup$ Jan 5 at 0:10

2 Answers 2

5
$\begingroup$

Let’s try narrowing this down to a smaller test case, to make the problem tractable, and ask why the vector isn’t looking right for one particular time. I often find that much easier than trying to keep up with the effects of a code change on whole tables.

from skyfield.api import load
ts = load.timescale()
t = ts.utc(2023, 12, 5, 15)
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
p = earth.at(t).observe(sun).apparent()
print(*p.xyz.km)
print('But FreeFlyer says:')
print('43414053.86821108 129268344.84962216 56035574.242660224')

This example is also easy to work with because the goal vector x,y,z values aren't over in a data file, but right here in the Stack Exchange text where we (and other readers in the future) can see them. The output is:

-43428254.79156197 -129264324.83888216 -56033831.191084124
But FreeFlyer says:
43414053.86821108 129268344.84962216 56035574.242660224

And we see that indeed the results disagree by quite a bit, as you pointed out.

Well, what effects might FreeFlyer be failing to apply? Maybe they don't account for aberration and deflection, and so aren't producing an apparent position at all. Let's try re-running the script, but leaving off .apparent():

-43414058.085826546 -129268339.04763176 -56035571.81497113
But FreeFlyer says:
43414053.86821108 129268344.84962216 56035574.242660224

That's better! Now the positions agree to within 5 meters, which is pretty good for the Sun-Earth distance. You can try fiddling with the time maybe if you want to see if you can bring it even closer (do they handle leap seconds the same?), but hopefully your main question about the difference is now answered!

Oh, I suppose my script uses a different ephemeris, DE421. Feel free to try it with DE430; the difference should be under half a meter or so.

$\endgroup$
1
  • 1
    $\begingroup$ you taught us how to pronounce ephemerides ten years ago and it seems you just can't stop teaching! Great answer, great package, great advice. $\endgroup$
    – uhoh
    Jan 6 at 6:46
4
$\begingroup$

angle_between_vectors_deg reminds me of an issue I encountered a few years ago. θ = arccos(x) loses precision where x ≈ ±1 (θ ≈ 0° or 180°). This Math SE answer gives a more reliable formula based on Kahan 2006 pp. 46-47. I suggest something like:

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))

However, to make the arccos formula fail at the 20-arcsecond level, I had to specify dtype=np.float32 for my test vectors. If yours are float64, then the problem may be elsewhere.

$\endgroup$
1
  • 1
    $\begingroup$ While this is very informative, I am indeed using float64 for the vector dtype, and using this better angle function yielded the same results. $\endgroup$ Jan 2 at 17:57

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .