<p><code>angle_between_vectors_deg</code> reminds me of an issue I encountered a few years ago.
θ = arccos(x) loses precision where x ≈ ±1 (θ ≈ 0° or 180°).
<a href="https://math.stackexchange.com/a/1782769">This Math SE answer</a> gives a more reliable formula based on <a href="https://people.eecs.berkeley.edu/%7Ewkahan/Mindless.pdf" rel="nofollow noreferrer">Kahan 2006</a> pp. 46-47.
I suggest something like:</p>
<pre><code>def angle_between_vectors_deg(a: np.ndarray, b: np.ndarray) -&gt; 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))
</code></pre>
<p>However, to make the arccos formula fail at the 20-arcsecond level,
I had to specify <code>dtype=np.float32</code> for my test vectors.
If yours are float64, then the problem may be elsewhere.</p>

CC BY-SA 4.0
2024-01-02T17:57:28.090Z
While this is very informative, I am indeed using float64 for the vector dtype, and using this better angle function yielded the same results.
Michael Bonnet
CC BY-SA 4.0
2023-12-30T22:03:16.030Z
2023-12-31T02:44:29.577Z
2023-12-31T02:44:29.577Z
Mike G
<p>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.</p>
<pre><code>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')
</code></pre>
<p>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:</p>
<pre><code>-43428254.79156197 -129264324.83888216 -56033831.191084124
But FreeFlyer says:
43414053.86821108 129268344.84962216 56035574.242660224
</code></pre>
<p>And we see that indeed the results disagree by quite a bit, as you pointed out.</p>
<p>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 <code>.apparent()</code>:</p>
<pre><code>-43414058.085826546 -129268339.04763176 -56035571.81497113
But FreeFlyer says:
43414053.86821108 129268344.84962216 56035574.242660224
</code></pre>
<p>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!</p>
<p>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.</p>

CC BY-SA 4.0
2024-01-06T06:46:19.980Z
you taught us how to pronounce [ephemerides ten years ago](https://english.stackexchange.com/q/370655/217285) and it seems you just can't stop teaching! Great answer, great package, great advice.
uhoh
CC BY-SA 4.0
2024-01-06T05:59:34.007Z
2024-01-06T05:59:34.007Z
Brandon Rhodes
<p>I have two functions that get TEME and GCRS Earth-Sun position vectors using Skyfield:</p>
<p><code>earth_sun_vectors.py</code></p>
<pre><code>&quot;&quot;&quot;Utilities for getting Earth-Sun vectors in a reference frame at a UTC epoch.&quot;&quot;&quot;

# 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) -&gt; np.ndarray:
    &quot;&quot;&quot;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.

    &quot;&quot;&quot;
    if not isinstance(epoch, datetime):
        raise TypeError(f&quot;arg epoch must be of type datetime, not {str(type(epoch))}&quot;)
    ts = load.timescale()
    epoch = epoch.replace(tzinfo=timezone.utc)
    t = ts.from_datetime(epoch)
    bodies = load(&quot;de430_1850-2150.bsp&quot;)
    with closing(bodies):
        sun = bodies[&quot;sun&quot;]
        earth = bodies[&quot;earth&quot;]
        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) -&gt; np.ndarray:
    &quot;&quot;&quot;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.

    &quot;&quot;&quot;
    if not isinstance(epoch, datetime):
        raise TypeError(f&quot;arg epoch must be of type datetime, not {str(type(epoch))}&quot;)
    ts = load.timescale()
    epoch = epoch.replace(tzinfo=timezone.utc)
    t = ts.from_datetime(epoch)
    bodies = load(&quot;de430_1850-2150.bsp&quot;)
    with closing(bodies):
        sun = bodies[&quot;sun&quot;]
        earth = bodies[&quot;earth&quot;]
        apparent = earth.at(t).observe(sun).apparent()
        vector = apparent.frame_xyz(TEME).km
        return np.array(vector)
</code></pre>
<p>To validate the functionality, I have generated test data with a common industry tool, FreeFlyer. The FreeFlyer <a href="https://ai-solutions.com/_help_Files/planetary_forces_files.htm" rel="nofollow noreferrer">documentation</a> clearly states that <code>FreeFlyer will use the De430.dat file by default</code>, hence why I use <code>de430_1850-2150.bsp</code> as my ephemeris file.</p>
<p>Below are my test functions, written using pytest:</p>
<p><code>test_earth_sun_vectors.py</code></p>
<pre><code>&quot;&quot;&quot;Unit tests for Earth-Sun vector functions.&quot;&quot;&quot;

# 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) -&gt; np.ndarray:
    &quot;&quot;&quot;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.

    &quot;&quot;&quot;
    return x / np.linalg.norm(x, axis=0)


def angle_between_vectors_deg(a: np.ndarray, b: np.ndarray) -&gt; float:
    &quot;&quot;&quot;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.

    &quot;&quot;&quot;
    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&quot;get_earth_sun_vector_gcrs_at_epoch_{datetime.now()}&quot;)

    test_data = pd.read_csv(
        &quot;earth_sun_vector_test_data.csv&quot;
    )

    test_count = 0

    # Validate the results
    for _, row in test_data.iterrows():
        test_count += 1
        
        # extract arg for function
        epoch = datetime.fromisoformat(row[&quot;Epoch_ISO8601_UTC&quot;][: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[&quot;Sun_Earth_MJ2000_X_km&quot;],  # MJ2000 is supposedly accurate to GCRS within milliarcseconds
                row[&quot;Sun_Earth_MJ2000_Y_km&quot;],
                row[&quot;Sun_Earth_MJ2000_Z_km&quot;],
            ]
        )
        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&quot;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}&quot;)

        # accurate to 1/10th of 1%
        np.testing.assert_allclose(test_position_vector, true_position_vector, rtol=0.001)
        assert np.abs(angle_between) &lt; 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&quot;get_earth_sun_vector_teme_at_epoch_{datetime.now()}&quot;)

    test_data = pd.read_csv(
        &quot;earth_sun_vector_test_data.csv&quot;
    )

    test_count = 0

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

        # extract arg for function
        epoch = datetime.fromisoformat(row[&quot;Epoch_ISO8601_UTC&quot;][: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[&quot;Sun_Earth_TEME_X_km&quot;],
                row[&quot;Sun_Earth_TEME_Y_km&quot;],
                row[&quot;Sun_Earth_TEME_Z_km&quot;],
            ]
        )
        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&quot;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}&quot;)

        # 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) &lt; 0.01

</code></pre>
<p>I have created a <a href="https://gist.github.com/MichaelBonnet/dc6f0d78f1a03629f12aeeb4bd01bbd6" rel="nofollow noreferrer">gist</a> with these and other files, including the test data csv and Docker-related files to minimize possible differences between machines.</p>
<p>An rtol of 0.001 is only 1/10th of 1% which seems not <em>that</em> 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.</p>
<p>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?</p>
<p>Thanks for any and all help.</p>

CC BY-SA 4.0
2023-12-30T14:49:40.030Z
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.
Greg Miller
CC BY-SA 4.0
2024-01-02T22:27:53.427Z
@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](https://rhodesmill.org/skyfield/time.html)
Michael Bonnet
CC BY-SA 4.0
2024-01-04T04:45:45.313Z
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.
Brandon Rhodes
CC BY-SA 4.0
2024-01-04T22:46:32.647Z
@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.
Michael Bonnet
CC BY-SA 4.0
2024-01-05T00:10:35.253Z
@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.
Michael Bonnet
CC BY-SA 4.0
2024-01-05T00:21:37.340Z
However, it's worth noting how small the magnitude of the difference is between the max error and min error: `0.005775612799308253 degrees` minimum, `0.005776897219398649 degrees` maximum.
Michael Bonnet
CC BY-SA 4.0
2024-01-05T01:09:26.667Z
I managed to generate ephemeris from NASA HORIZONS. Checking against that ephemeris in CSV-of-position-vectors-and-epoch form gives an average angle between predicted and actual sun vectors in GCRS of ~0.00081256 degrees - that is, an order of magnitude more accurate.
Michael Bonnet
<the-sun><python><positional-astronomy><ephemeris><skyfield>
Resolving variance between Skyfield-generated sun vectors and test data when both use same ephemeris
Michael Bonnet
# Resolving variance between Skyfield-generated sun vectors and test data when both use same ephemeris

(Asked by: username_0 on 2023-12-29)

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](https://ai-solutions.com/_help_Files/planetary_forces_files.htm) 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](https://gist.github.com/username_0/dc6f0d78f1a03629f12aeeb4bd01bbd6) 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.

**Tags:** the-sun, python, positional-astronomy, ephemeris, skyfield

> **Comments:**
> **username_1**: 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.  
> **username_0**: @username_1 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](https://rhodesmill.org/skyfield/time.html)  
> **username_2**: 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.  
> **username_0**: @username_2 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.  
> **username_0**: @username_2 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.  
> **username_0**: However, it's worth noting how small the magnitude of the difference is between the max error and min error: `0.005775612799308253 degrees` minimum, `0.005776897219398649 degrees` maximum.  
> **username_0**: I managed to generate ephemeris from NASA HORIZONS. Checking against that ephemeris in CSV-of-position-vectors-and-epoch form gives an average angle between predicted and actual sun vectors in GCRS of ~0.00081256 degrees - that is, an order of magnitude more accurate.  

---

## 2 answers

### Answer by username_2 on 2024-01-06. Score: 5 | Accepted answer

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.

> **Comments:**
> **username_4**: you taught us how to pronounce [ephemerides ten years ago](https://english.stackexchange.com/q/370655/217285) and it seems you just can't stop teaching! Great answer, great package, great advice.  

---

### Answer by username_3 on 2023-12-30. Score: 4

`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](https://math.stackexchange.com/a/1782769) gives a more reliable formula based on [Kahan 2006](https://people.eecs.berkeley.edu/%7Ewkahan/Mindless.pdf) 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.

> **Comments:**
> **username_0**: While this is very informative, I am indeed using float64 for the vector dtype, and using this better angle function yielded the same results.  

---

