# Fixing satellite eclipse equations from textbooks that are seemingly failing

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

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

earth_sun_distance = norm(earth_sun_position_vector_eci)

in_umbra = False
in_penumbra = False

if np.dot(earth_object_position_vector_eci, earth_sun_position_vector_eci) < 0.0:
sathoriz = norm(earth_object_position_vector_eci)*np.cos(angle)
satvert  = norm(earth_object_position_vector_eci)*np.sin(angle)
penvert = np.tan(penumbra_angle)*(x + sathoriz)
if satvert <= penvert:
in_penumbra = True
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():

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

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

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

def check_umbra( delta_ps, Dp, proj_scalar, rej_norm):
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):
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 )

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

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:

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.

• What exactly and where are the discrepancies? Dec 12, 2023 at 21:13
• 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 Dec 12, 2023 at 21:59
• @planetmaker I have created a gist with files of relevant data used for testing as well as information associated with failed test cases. Dec 12, 2023 at 22:52
• @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. Dec 12, 2023 at 22:52

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.

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

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

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.

• 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.
– uhoh
Dec 14, 2023 at 23:50