4
$\begingroup$

I have tried implementing equations from two different textbooks in Python to find if a satellite is in eclipse, but both fail a significant percentage of the time when tested against data from industry tools like FreeFlyer. I have the data in CSV form as the positions of spacecraft and the Sun in TEME coordinate vectors at individual timesteps, along with whether the spacecraft is in eclipse at that timestep.

The first equation is algorithm 34 from David Vallado's Fundamentals of Astrodynamics and Applications:

import numpy as np

earth_radius_km = 6378.1366
sun_radius_km = 695700

def norm(x: np.ndarray) -> float:
    return np.linalg.norm(x, axis=0)

def unit_vector(x: np.ndarray) -> np.ndarray:
    return x / norm(x)

def angle_between_vectors_rad(a: np.ndarray, b: np.ndarray) -> float:
    return np.arccos(np.sum(unit_vector(a) * unit_vector(b), axis=0))

def check_object_shadows_vallado_algorithm_34(earth_object_position_vector_eci: np.ndarray, earth_sun_position_vector_eci: np.ndarray) -> tuple[bool, bool]:
    earth_sun_distance = norm(earth_sun_position_vector_eci)
    umbra_angle = np.arctan((sun_radius_km-earth_radius_km)/earth_sun_distance)
    penumbra_angle = np.arctan((sun_radius_km+earth_radius_km)/earth_sun_distance)

    in_umbra = False
    in_penumbra = False
    
    if np.dot(earth_object_position_vector_eci, earth_sun_position_vector_eci) < 0.0:
        angle = angle_between_vectors_rad(-earth_sun_position_vector_eci, earth_object_position_vector_eci)
        sathoriz = norm(earth_object_position_vector_eci)*np.cos(angle)
        satvert  = norm(earth_object_position_vector_eci)*np.sin(angle)
        x = earth_radius_km/np.sin(penumbra_angle)
        penvert = np.tan(penumbra_angle)*(x + sathoriz)
        if satvert <= penvert:
            in_penumbra = True
            y = earth_radius_km/np.sin(umbra_angle)
            umbvert = np.tan(umbra_angle)*(y-sathoriz)
            if satvert <= umbvert:
                in_umbra = True

    return in_penumbra, in_umbra

I check the above with the following:

def test_vallado_alg_34():
    validation_data = pd.read_csv("validation_data.csv")

    total_checks = 0
    correct_checks = 0
    incorrect_checks_actual_sunlit = 0
    incorrect_checks_actual_eclipse = 0
    incorrect_checks_actual_sunlit_penumbra = 0
    incorrect_checks_actual_sunlit_umbra = 0
    correct_check_indices = []
    incorrect_check_indices = []

    # Validate the results
    for index, row in validation_data.iterrows():
        total_checks += 1
        # extract args for function
        earth_sun_position_vector_teme = np.array(
            [row["Sun_Position_TEME_X_km"], row["Sun_Position_TEME_Y_km"], row["Sun_Position_TEME_Z_km"]]
        )
        earth_object_position_vector_teme = np.array(
            [row["Cartesian_TEME_X_km"], row["Cartesian_TEME_Y_km"], row["Cartesian_TEME_Z_km"]]
        )

        # get actual value
        in_penumbra, in_umbra = check_object_shadows_vallado_algorithm_34(earth_object_position_vector_teme, earth_sun_position_vector_teme)
        is_in_eclipse_actual = in_penumbra or in_umbra

        # get expected value
        is_in_eclipse_expected = bool(row["Sat_InShadow_bool"]) # represented as 0 or 1 in csv

        # assert
        if is_in_eclipse_actual == is_in_eclipse_expected:
            correct_checks += 1
            correct_check_indices.append(index)
        else:
            incorrect_check_indices.append(index)
            if is_in_eclipse_actual is True and is_in_eclipse_expected is False:
                incorrect_checks_actual_sunlit += 1
                if in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra += 1
                if in_umbra:
                    incorrect_checks_actual_sunlit_umbra += 1
            elif is_in_eclipse_actual is False and is_in_eclipse_expected is True:
                incorrect_checks_actual_eclipse += 1

    print(f"# Total Checks     : {total_checks}")
    print(f"# Correct Checks   : {correct_checks}")
    print(f"# Incorrect Checks : {total_checks - correct_checks}")
    print(f"% Correct          : {correct_checks / total_checks}")
    print(f"# of times we thought we were in eclipse but actually were sunlit : {incorrect_checks_actual_sunlit}")
    print(f"# of times we thought we were sunlit but actually were in eclipse : {incorrect_checks_actual_eclipse}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : {incorrect_checks_actual_sunlit_umbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : {incorrect_checks_actual_sunlit_penumbra}")

Results when running tests for 14 different Earth-orbiting objects from VLEO to GEO across a 24-hour period for 30 second time steps:

# Total Checks     : 43200
# Correct Checks   : 27021
# Incorrect Checks : 16179
% Correct          : 0.6254861111111111
# of times we thought we were in eclipse but actually were sunlit : 8031
# of times we thought we were sunlit but actually were in eclipse : 8148
# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : 7923
# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : 8031

To understand the above code, here is algorithm 34, illustrated in his MATLAB code:

Part 1 of Algorithm 34, from Fundamentals of Astrodynamics and Applications (Vallado 2013) Part 2 of Algorithm 34, from Fundamentals of Astrodynamics and Applications (Vallado 2013)


The second equation is from Alfonso Gonzalez's Astrodynamics with Python and adapted from his example code - see check_eclipse, check_umbra, and check_penumbra

import math
import numpy as np

earth_radius_km = 6378.1366
sun_radius_km = 695700

def norm(x: np.ndarray) -> float:
    return np.linalg.norm(x, axis=0)

def check_umbra( delta_ps, Dp, proj_scalar, rej_norm):
    sun_diameter = sun_radius_km * 2
    Xu     = ( Dp * delta_ps ) / ( sun_diameter - Dp )
    alphau = math.asin( Dp / ( 2 * Xu ) )
    zeta   = ( Xu - proj_scalar ) * math.tan( alphau )
    return rej_norm <= zeta

def check_penumbra( delta_ps, Dp, proj_scalar, rej_norm):
    sun_diameter = sun_radius_km * 2
    Xp     = ( Dp * delta_ps ) / ( sun_diameter + Dp )
    alphap = math.asin( Dp / ( 2 * Xp ) )
    kappa  = ( Xp + proj_scalar ) * math.tan( alphap )
    return rej_norm <= kappa

def check_eclipse( sun_earth_position_vector, earth_spacecraft_position_vector ):
    in_penumbra = False
    in_umbra = False
    
    r_sun2body  = sun_earth_position_vector
    delta_ps    = norm( r_sun2body )
    s_hat       = r_sun2body / delta_ps
    proj_scalar = np.dot( earth_spacecraft_position_vector, s_hat )

    if proj_scalar <= 0.0:
        return in_penumbra, in_umbra

    proj     = proj_scalar * s_hat
    rej_norm = norm( earth_spacecraft_position_vector - proj )
    
    earth_diameter_km = earth_radius_km * 2

    in_umbra = check_umbra( delta_ps, earth_diameter_km, proj_scalar, rej_norm )
    in_penumbra = check_penumbra( delta_ps, earth_diameter_km, proj_scalar, rej_norm )

    return in_penumbra, in_umbra

I check the above with the following:

def test_gonzalez_awp():
    validation_data = pd.read_csv("validation_data.csv")

    total_checks = 0
    correct_checks = 0
    incorrect_checks_actual_sunlit = 0
    incorrect_checks_actual_eclipse = 0
    incorrect_checks_actual_sunlit_penumbra = 0
    incorrect_checks_actual_sunlit_umbra = 0
    incorrect_checks_actual_sunlit_penumbra_and_umbra = 0
    correct_check_indices = []
    incorrect_check_indices = []

    # Validate the results
    for index, row in validation_data.iterrows():
        total_checks += 1
        # extract args for function
        earth_sun_position_vector_teme = np.array(
            [row["Sun_Position_TEME_X_km"], row["Sun_Position_TEME_Y_km"], row["Sun_Position_TEME_Z_km"]]
        )
        earth_object_position_vector_teme = np.array(
            [row["Cartesian_TEME_X_km"], row["Cartesian_TEME_Y_km"], row["Cartesian_TEME_Z_km"]]
        )

        # get actual value
        in_penumbra, in_umbra = check_eclipse(-earth_sun_position_vector_teme, earth_object_position_vector_teme) # using negative earth-sun vector because equation expects sun-earth vectpr
        is_in_eclipse_actual = in_penumbra or in_umbra

        # get expected value
        is_in_eclipse_expected = bool(row["Sat_InShadow_bool"]) # represented as 0 or 1 in csv

        # assert
        if is_in_eclipse_actual == is_in_eclipse_expected:
            correct_checks += 1
            correct_check_indices.append(index)
        else:
            incorrect_check_indices.append(index)
            if is_in_eclipse_actual is True and is_in_eclipse_expected is False:
                incorrect_checks_actual_sunlit += 1
                if in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra += 1
                if in_umbra:
                    incorrect_checks_actual_sunlit_umbra += 1
                if in_umbra and in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra_and_umbra += 1
            elif is_in_eclipse_actual is False and is_in_eclipse_expected is True:
                incorrect_checks_actual_eclipse += 1

    print(f"# Total Checks     : {total_checks}")
    print(f"# Correct Checks   : {correct_checks}")
    print(f"# Incorrect Checks : {total_checks - correct_checks}")
    print(f"% Correct          : {correct_checks / total_checks}")
    print(f"# of times we thought we were in eclipse but actually were sunlit : {incorrect_checks_actual_sunlit}")
    print(f"# of times we thought we were sunlit but actually were in eclipse : {incorrect_checks_actual_eclipse}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : {incorrect_checks_actual_sunlit_umbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : {incorrect_checks_actual_sunlit_penumbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra+penumbra) : {incorrect_checks_actual_sunlit_penumbra}")

Results when running tests for 14 different Earth-orbiting objects from VLEO to GEO across a 24-hour period for 30 second time steps:

# Total Checks     : 43200
# Correct Checks   : 27021
# Incorrect Checks : 16179
% Correct          : 0.6254861111111111
# of times we thought we were in eclipse but actually were sunlit : 0
# of times we thought we were sunlit but actually were in eclipse : 8148
# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : 0
# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : 0

To aid in understanding Gonzalez/my adaptation of Gonzalez, see my notes from his video:

Notes from Gonzalez's video on calculating eclipse, equations only. Notes from Gonzalez's video on calculating eclipse, diagrams/visuals alongside equations for umbra. Notes from Gonzalez's video on calculating eclipse, diagrams/visuals alongside equations for penumbra.


My questions are:

  1. Have I translated these equations/their code into Python incorrectly somehow, leading to the discrepancies?

  2. Am I possibly making mistakes in my test approach?

  3. Do I potentially just have fundamental flaws in my helper functions like norm()?

EDIT: I have uploaded files to a gist for context.

eclipse_checking_validation_data.csv is the data I use for checking, consisting of a satellite position vector, a sun position vector, and whether or not the satellite is in eclipse.

gonzalez_validation_data_with_results.csv contains the provided values, expected values, and calculated values for every failed test case using the Gonzalez method.

vallado_validation_data_with_results.csv contains the provided values, expected values, and calculated values for every failed test case using the Vallado method.

$\endgroup$
4
  • $\begingroup$ What exactly and where are the discrepancies? $\endgroup$ Dec 12, 2023 at 21:13
  • 1
    $\begingroup$ The code for Vallado's book is available in a few different languages, you should compare yours to his: celestrak.org/software/vallado-sw.php $\endgroup$ Dec 12, 2023 at 21:59
  • $\begingroup$ @planetmaker I have created a gist with files of relevant data used for testing as well as information associated with failed test cases. $\endgroup$ Dec 12, 2023 at 22:52
  • $\begingroup$ @GregMiller I have followed Vallado's code as closely as I can given my use case, but his code is not available in Python, so I have had to write it myself. $\endgroup$ Dec 12, 2023 at 22:52

1 Answer 1

8
+100
$\begingroup$

The issues in each approach stemmed from bad assumptions, lack of source specificity, and badly labeled data.

Both approaches were essentially doing the same math, just with different labeling/organization. They failed and succeeded on the exact same cases.

As for badly labeled data: in my test data, the Sun_Position_<ref_frame>_* columns of the CSV should have been named Sun_Earth_Position_<ref_frame>_*. That is, the vectors were from the Sun to the Earth, as opposed to from the Earth to the Sun as I had assumed.

Vallado Approach

Given the data labeling issues above, my issue was not multiplying what was actually the Sun-Earth vector by -1 to get the Earth-Sun vector. Once I did that to the vector before feeding it into Vallado's algorithm, it worked extremely well*.

I should have:

  • Relabeled the data - Sun_Position_TEME_X_km to Sun_Earth_Position_TEME_X_km and so forth

  • Used this test function:

def test_vallado_alg_34():
    validation_data = pd.read_csv("validation_data.csv")

    total_checks = 0
    correct_checks = 0
    incorrect_checks_actual_sunlit = 0
    incorrect_checks_actual_eclipse = 0
    incorrect_checks_actual_sunlit_penumbra = 0
    incorrect_checks_actual_sunlit_umbra = 0
    correct_check_indices = []
    incorrect_check_indices = []

    # Validate the results
    for index, row in validation_data.iterrows():
        total_checks += 1
        # extract args for function
        earth_object_position_vector_teme = np.array(
            [row["Cartesian_TEME_X_km"], row["Cartesian_TEME_Y_km"], row["Cartesian_TEME_Z_km"]]
        )
        sun_earth_position_vector_teme = np.array(
            [row["Sun_Earth_Position_TEME_X_km"], row["Sun_Earth_Position_TEME_Y_km"], row["Sun_Earth_Position_TEME_Z_km"]]
        )
        earth_sun_position_vector_teme = -sun_earth_position_vector_teme

        # get actual value
        in_penumbra, in_umbra = check_object_shadows_vallado_algorithm_34(earth_object_position_vector_teme, earth_sun_position_vector_teme)
        is_in_eclipse_actual = in_penumbra or in_umbra

        # get expected value
        is_in_eclipse_expected = bool(row["Sat_InShadow_bool"]) # represented as 0 or 1 in csv

        # assert
        if is_in_eclipse_actual == is_in_eclipse_expected:
            correct_checks += 1
            correct_check_indices.append(index)
        else:
            incorrect_check_indices.append(index)
            if is_in_eclipse_actual is True and is_in_eclipse_expected is False:
                incorrect_checks_actual_sunlit += 1
                if in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra += 1
                if in_umbra:
                    incorrect_checks_actual_sunlit_umbra += 1
            elif is_in_eclipse_actual is False and is_in_eclipse_expected is True:
                incorrect_checks_actual_eclipse += 1

    print(f"# Total Checks     : {total_checks}")
    print(f"# Correct Checks   : {correct_checks}")
    print(f"# Incorrect Checks : {total_checks - correct_checks}")
    print(f"% Correct          : {correct_checks / total_checks}")
    print(f"# of times we thought we were in eclipse but actually were sunlit : {incorrect_checks_actual_sunlit}")
    print(f"# of times we thought we were sunlit but actually were in eclipse : {incorrect_checks_actual_eclipse}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : {incorrect_checks_actual_sunlit_umbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : {incorrect_checks_actual_sunlit_penumbra}")

Gonzalez Approach

Gonzalez's approach explicitly called for using Sun-Earth vectors, not Earth-Sun vectors. Because of my wrong assumptions, I was multiplying what was actually the Sun-Earth vector by -1, getting the Earth-Sun vector and throwing off the calculations. Once I removed the multiplication by -1 before feeding it into Gonzalez's algorithm, it worked extremely well*.

I should have:

  • Relabeled the data - Sun_Position_TEME_X_km to Sun_Earth_Position_TEME_X_km and so forth

  • Used this test function:

def test_gonzalez_awp():
    validation_data = pd.read_csv("validation_data.csv")

    total_checks = 0
    correct_checks = 0
    incorrect_checks_actual_sunlit = 0
    incorrect_checks_actual_eclipse = 0
    incorrect_checks_actual_sunlit_penumbra = 0
    incorrect_checks_actual_sunlit_umbra = 0
    incorrect_checks_actual_sunlit_penumbra_and_umbra = 0
    correct_check_indices = []
    incorrect_check_indices = []

    # Validate the results
    for index, row in validation_data.iterrows():
        total_checks += 1
        # extract args for function
        sun_earth_position_vector_teme = np.array(
            [row["Sun_Earth_Position_TEME_X_km"], row["Sun_Earth_Position_TEME_Y_km"], row["Sun_Earth_Position_TEME_Z_km"]]
        )
        earth_object_position_vector_teme = np.array(
            [row["Cartesian_TEME_X_km"], row["Cartesian_TEME_Y_km"], row["Cartesian_TEME_Z_km"]]
        )

        # get actual value
        in_penumbra, in_umbra = check_eclipse(sun_earth_position_vector_teme, earth_object_position_vector_teme)
        is_in_eclipse_actual = in_penumbra or in_umbra

        # get expected value
        is_in_eclipse_expected = bool(row["Sat_InShadow_bool"]) # represented as 0 or 1 in csv

        # assert
        if is_in_eclipse_actual == is_in_eclipse_expected:
            correct_checks += 1
            correct_check_indices.append(index)
        else:
            incorrect_check_indices.append(index)
            if is_in_eclipse_actual is True and is_in_eclipse_expected is False:
                incorrect_checks_actual_sunlit += 1
                if in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra += 1
                if in_umbra:
                    incorrect_checks_actual_sunlit_umbra += 1
                if in_umbra and in_penumbra:
                    incorrect_checks_actual_sunlit_penumbra_and_umbra += 1
            elif is_in_eclipse_actual is False and is_in_eclipse_expected is True:
                incorrect_checks_actual_eclipse += 1

    print(f"# Total Checks     : {total_checks}")
    print(f"# Correct Checks   : {correct_checks}")
    print(f"# Incorrect Checks : {total_checks - correct_checks}")
    print(f"% Correct          : {correct_checks / total_checks}")
    print(f"# of times we thought we were in eclipse but actually were sunlit : {incorrect_checks_actual_sunlit}")
    print(f"# of times we thought we were sunlit but actually were in eclipse : {incorrect_checks_actual_eclipse}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra)          : {incorrect_checks_actual_sunlit_umbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc penumbra)       : {incorrect_checks_actual_sunlit_penumbra}")
    print(f"# of times we thought we were in eclipse but actually were sunlit (calc umbra+penumbra) : {incorrect_checks_actual_sunlit_penumbra}")

Essentially, both are correct, carrying out the same mathematical process with different assumptions (sun-earth vector vs earth-sun vector). I'm going to go with the Vallado version for a few reasons:

  • It comes from a widely-used textbook and can be related to existing code, making maintenance easier down the line
  • A bit more succinct than Gonzalez' approach
  • In my opinion, more readable (even in example code, Vallado's variables are much better named)

* using all test cases, these approaches yielded 45-57 failures out of ~40000-43000 test cases. All inconsistencies are happening around the penumbra <-> sunlight transition. I am unable to chalk this up to the industry tool test data taking into account earth and/or sun oblateness, as the tool's eclipse-checking function explicitly states it assumes spherical Earth and Sun. I do not currently know the root cause of the differences.

$\endgroup$
1
  • 7
    $\begingroup$ Congrats, well done! I've often found that carefully formulating and composing a Stack Exchange question helps me to eventually find the solution myself, and it's great when people take the time to come back and post the solution as an answer post for future readers. $\endgroup$
    – uhoh
    Dec 14, 2023 at 23:50

You must log in to answer this question.

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