Simulate a CubeSat Power System
Model a CubeSat power system in orbit using SGP4 and real TLE data.
Introduction
30 Second Summary
Satellites don't have wall outlets. Every watt of power a spacecraft uses has to come from its solar panels, and every orbit brings 30+ minutes of darkness where those panels produce nothing.
In this project, you will build a Python simulation that models how a small satellite's power system behaves in orbit. You'll propagate a real satellite's orbit, detect when it enters Earth's shadow, and track whether the battery can survive each eclipse period.
What You'll Build
A Python-based CubeSat power system simulator that propagates a satellite orbit using real TLE data, detects eclipse periods with the cylindrical shadow model, models solar panel generation, and visualizes battery state of charge across multiple orbits.
By the end of this project, you'll have:
- A simulation engine that propagates a satellite orbit using real Two-Line Element data and the SGP4 algorithm.
- An eclipse detection algorithm using the cylindrical shadow model to determine when the satellite enters Earth's shadow.
- A 4-panel visualization showing orbital altitude, sunlight and eclipse windows, power generation, and battery state of charge over time.
- Secret Mission: Add mission phase load scheduling that varies power consumption by phase and implements autonomous battery protection.
Are there any prerequisites?
You should be comfortable with basic Python (functions, loops, and imports) and have uv installed as your Python package manager. No prior knowledge of orbital mechanics or space systems is required.
Before We Start
Before We Start...
Set Up the Project
Before you can simulate how a CubeSat behaves in orbit, you need a Python project with the right scientific libraries installed.
In this step, you'll use uv to scaffold a project and install three dependencies: sgp4 for orbit propagation, NumPy for numerical computation, and Matplotlib for visualization.
In this step, get ready to:
- Create and initialize a new Python project with uv.
- Install sgp4, numpy, and matplotlib as dependencies.
- Run a test script confirming all libraries load and sgp4 works.
Create the project directory
Every simulation project starts with a dedicated folder and a proper dependency manager. This keeps your libraries isolated from anything else on your system.
- Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
- Type Terminal and press Enter to open it.
- Move to your Desktop by running this command:
cd ~/Desktop- Create and initialize your project by running this command:
uv init cubesat-power-simWhat does this command do?
The uv init command creates a new folder called cubesat-power-sim with a pyproject.toml (your project configuration), a main.py starter file, and a few other scaffolding files like .python-version and README.md.
- Move into the new project folder by running:
cd cubesat-power-sim- Confirm the project files exist by running:
lsYou should see main.py, pyproject.toml, README.md, and .python-version in the output.
Don't see the files?
Make sure you ran uv init cubesat-power-sim (not just uv init) from your Desktop. If you see an error about uv not being found, confirm uv is installed correctly.
help me troubleshoot why uv init did not create my project files
Install the simulation dependencies
Your simulation needs three libraries. The SGP4 library implements the same orbit propagation algorithm that NORAD uses to track every satellite in space. NumPy handles the heavy math, and Matplotlib draws the plots.
- Install all three dependencies by running:
uv add sgp4 numpy matplotlibWhat does this command do?
The uv add command downloads each library, records it in pyproject.toml, and locks the exact versions in a uv.lock file. This ensures your project is reproducible on any machine.
You should see output showing each package being resolved and installed. The process may take a few seconds while uv downloads everything.
Seeing a network or resolution error?
Check your internet connection. If you see a version conflict, try running uv add sgp4 numpy matplotlib again. Transient network issues sometimes cause partial installs.
help me fix a uv add error when installing sgp4 numpy matplotlib
Verify the setup
A quick test script confirms that all three libraries load correctly and that sgp4 can create a satellite object from real orbital data.
- Open main.py inside the cubesat-power-sim folder in your code editor.
- Replace everything in the file with the following code:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
print("Libraries loaded successfully.")
# Confirm sgp4 C++ acceleration
test_sat = Satrec.twoline2rv(
"1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990",
"2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
)
print(f"sgp4 satellite object created: {test_sat}")
print("Setup complete!")What does this code do?
The imports pull in numpy, matplotlib, and two key sgp4 functions: Satrec (creates satellite objects from orbital data) and jday (converts dates to Julian format for orbital calculations).
The Satrec.twoline2rv() call uses real TLE (Two-Line Element) data from the International Space Station. If sgp4 can parse this and create a satellite object, your installation is working.
- Save main.py.
- Run the script by executing this command in your terminal:
uv run python main.pyYou should see three lines printed: Libraries loaded successfully., a line showing the sgp4 satellite object was created, and Setup complete!.
Seeing an import error?
If you see ModuleNotFoundError, make sure you run with uv run python main.py (not just python main.py). The uv run prefix ensures the script uses the project's installed dependencies.
help me fix an import error when running my cubesat power sim setup script
✔️ Awesome, I've got everything!
Great. Double check you've saved main.py and you're ready to move on.
ⓧ I'd like to double check the full code
Here is the complete main.py at this point:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
print("Libraries loaded successfully.")
# Confirm sgp4 C++ acceleration
test_sat = Satrec.twoline2rv(
"1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990",
"2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
)
print(f"sgp4 satellite object created: {test_sat}")
print("Setup complete!")Your project is scaffolded and all dependencies are ready. Next up, you'll use sgp4 to propagate a real satellite orbit and watch it trace a path around Earth.
Propagate a Satellite Orbit with SGP4
Your project environment is ready. Now it's time to make your simulation actually do something: track a real satellite as it circles Earth.
In this step, you'll feed real orbital data into the SGP4 algorithm and propagate a satellite's position over several orbits. SGP4 is the same algorithm NASA and NORAD use to track every object in orbit around Earth.
In this step, get ready to:
- Define real ISS orbital data using the Two-Line Element (TLE) format.
- Write a propagation loop that computes the satellite's position every 30 seconds for 6 orbits.
- Plot the orbital altitude to verify a sinusoidal pattern.
Define TLE data and simulation constants
A TLE (Two-Line Element set) is a standardized format that encodes a satellite's orbital parameters. Satellites like the ISS have their TLE published publicly on CelesTrak.
You'll replace the test code from Step 1 with real configuration constants and ISS orbital data that SGP4 can use to compute positions.
- In main.py, select and delete everything below the import lines (the print statements and the test_sat block).
- Replace it with the following configuration constants and TLE data:
# --- Configuration ---
EARTH_RADIUS_KM = 6371.0
SOLAR_CONSTANT = 1361.0 # W/m^2
DEG2RAD = np.pi / 180.0
# ISS TLE data (example - replace with current TLE for best accuracy)
TLE_LINE1 = "1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990"
TLE_LINE2 = "2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
# Simulation parameters
SIM_START = (2024, 2, 14, 12, 0, 0) # Year, month, day, hour, min, sec
NUM_ORBITS = 6
TIME_STEP_SEC = 30What does this code do?
EARTH_RADIUS_KM and SOLAR_CONSTANT are physical constants you'll use later for eclipse detection and power calculations.
TLE_LINE1 and TLE_LINE2 encode the ISS's orbital shape, inclination, altitude, and current position as of February 2024. SGP4 reads this format directly.
The simulation parameters tell SGP4 when to start the simulation and how finely to sample positions. Six orbits at 30-second intervals gives you roughly 9 hours of data with about 1,080 data points.
- Save main.py.
What is a TLE?
A Two-Line Element set is how the U.S. military (NORAD) distributes satellite orbit data. Every tracked object in space has a TLE that gets updated every few hours. The two lines pack all the orbital parameters into a fixed-width text format designed in the 1960s.
CelesTrak distributes current TLEs for thousands of satellites. The ISS TLE in your code is from February 2024. For a real mission, you'd pull the latest TLE to get the most accurate prediction.
Write the orbit propagation function
SGP4 takes a TLE and a timestamp, then returns the satellite's exact position and velocity at that moment. Your propagation function will call SGP4 repeatedly at 30-second intervals to trace the satellite's path through space.
- Below the simulation parameters in main.py, add the following function:
# --- Orbit Propagation ---
def propagate_orbit(tle_line1, tle_line2, start_time, num_orbits, dt_sec):
"""Propagate satellite orbit using SGP4 and return positions over time."""
satellite = Satrec.twoline2rv(tle_line1, tle_line2)
# Compute orbital period from mean motion (radians/minute)
orbital_period_min = 2 * np.pi / satellite.no_kozai
total_time_min = orbital_period_min * num_orbits
num_steps = int(total_time_min * 60 / dt_sec)
# Start time as Julian Date
jd_start, fr_start = jday(*start_time)
times_min = np.zeros(num_steps)
positions_km = np.zeros((num_steps, 3))
altitudes_km = np.zeros(num_steps)
jd_array = np.zeros(num_steps)
fr_array = np.zeros(num_steps)
for i in range(num_steps):
t_min = i * dt_sec / 60.0
times_min[i] = t_min
# Compute Julian Date for this timestep
fr = fr_start + (t_min / 1440.0) # Convert minutes to fraction of day
jd = jd_start
# Handle day overflow
while fr >= 1.0:
fr -= 1.0
jd += 1.0
jd_array[i] = jd
fr_array[i] = fr
e, r, v = satellite.sgp4(jd, fr)
if e != 0:
print(f"SGP4 error at step {i}: error code {e}")
continue
positions_km[i] = r
altitudes_km[i] = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) - EARTH_RADIUS_KM
return times_min, positions_km, altitudes_km, jd_array, fr_arrayWhat does this code do?
- Satrec.twoline2rv() creates a satellite object from the two TLE lines. This object holds all the orbital parameters SGP4 needs.
- satellite.no_kozai gives the mean motion in radians per minute. Dividing 2*pi by this gives the orbital period. For the ISS, this is about 92 minutes.
- jday() converts a calendar date into Julian Date format, which SGP4 requires. It returns two values: the integer day and the fractional part.
- satellite.sgp4(jd, fr) returns position (km) and velocity (km/s) in the TEME coordinate frame. The error code e is 0 when propagation succeeds.
- Altitude is computed by finding the distance from Earth's center (sqrt(x^2 + y^2 + z^2)) and subtracting Earth's radius (6371 km).
- Save main.py.
Seeing a syntax error?
Make sure the function definition starts at the leftmost indentation level (not indented inside anything else). The def keyword should be at column 0.
Check that the f-string on the error line uses matching quotes. The outer quotes and the inner quotes must be different types.
help me debug a syntax error in my propagate_orbit function
Plot the altitude to verify the orbit
A correctly propagated LEO orbit should show altitude oscillating between perigee (lowest point) and apogee (highest point) in a roughly sinusoidal pattern. Plotting this is how you confirm the simulation is working before building on top of it.
- Below the propagate_orbit function, add a main() function that calls the propagator and plots the result:
def main():
print("Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f"Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f"Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
# Plot altitude
plt.figure(figsize=(12, 3))
plt.plot(times_min / 60.0, altitudes_km, "b-", linewidth=0.8)
plt.xlabel("Time (hours)")
plt.ylabel("Altitude (km)")
plt.title("Orbital Altitude over Time")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
main()What does this code do?
The main() function calls your propagator with the ISS TLE and simulation parameters, then prints a summary of the results. The altitude range for the ISS should be roughly 410-430 km.
The matplotlib plot shows altitude on the Y-axis and time in hours on the X-axis. You should see 6 smooth wave-like cycles, one per orbit.
- Save main.py.
- Run the simulation by executing this command:
uv run python main.pyYou should see terminal output showing the altitude range (roughly 410-430 km for the ISS) and a matplotlib window with a sinusoidal altitude plot cycling 6 times over about 9 hours.
Plot window not appearing?
On some systems, matplotlib needs a display backend. If you're running in a headless environment, try adding plt.savefig("altitude.png") before plt.show() and open the saved image instead.
If you see an altitude range that looks wrong (negative values or thousands of km), double check that your TLE strings are copied exactly with no extra spaces or missing characters.
help me debug my SGP4 orbit propagation
✔️ Awesome, I've got everything!
Great. Double check you've saved your file before moving on.
ⓧ I'd like to double check the full code
Here is the complete main.py at this point:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
# --- Configuration ---
EARTH_RADIUS_KM = 6371.0
SOLAR_CONSTANT = 1361.0 # W/m^2
DEG2RAD = np.pi / 180.0
# ISS TLE data (example - replace with current TLE for best accuracy)
TLE_LINE1 = "1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990"
TLE_LINE2 = "2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
# Simulation parameters
SIM_START = (2024, 2, 14, 12, 0, 0) # Year, month, day, hour, min, sec
NUM_ORBITS = 6
TIME_STEP_SEC = 30
# --- Orbit Propagation ---
def propagate_orbit(tle_line1, tle_line2, start_time, num_orbits, dt_sec):
"""Propagate satellite orbit using SGP4 and return positions over time."""
satellite = Satrec.twoline2rv(tle_line1, tle_line2)
# Compute orbital period from mean motion (radians/minute)
orbital_period_min = 2 * np.pi / satellite.no_kozai
total_time_min = orbital_period_min * num_orbits
num_steps = int(total_time_min * 60 / dt_sec)
# Start time as Julian Date
jd_start, fr_start = jday(*start_time)
times_min = np.zeros(num_steps)
positions_km = np.zeros((num_steps, 3))
altitudes_km = np.zeros(num_steps)
jd_array = np.zeros(num_steps)
fr_array = np.zeros(num_steps)
for i in range(num_steps):
t_min = i * dt_sec / 60.0
times_min[i] = t_min
# Compute Julian Date for this timestep
fr = fr_start + (t_min / 1440.0) # Convert minutes to fraction of day
jd = jd_start
# Handle day overflow
while fr >= 1.0:
fr -= 1.0
jd += 1.0
jd_array[i] = jd
fr_array[i] = fr
e, r, v = satellite.sgp4(jd, fr)
if e != 0:
print(f"SGP4 error at step {i}: error code {e}")
continue
positions_km[i] = r
altitudes_km[i] = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) - EARTH_RADIUS_KM
return times_min, positions_km, altitudes_km, jd_array, fr_array
def main():
print("Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f"Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f"Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
# Plot altitude
plt.figure(figsize=(12, 3))
plt.plot(times_min / 60.0, altitudes_km, "b-", linewidth=0.8)
plt.xlabel("Time (hours)")
plt.ylabel("Altitude (km)")
plt.title("Orbital Altitude over Time")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
main()Your orbit propagator is working. You now have real satellite positions over time. Next, you'll use these positions to figure out when the satellite passes through Earth's shadow.
Detect Eclipse Periods
Your orbit propagator is working. You can see the satellite's altitude varying over time as it traces its path around Earth.
Now you need to answer the question every power engineer asks first: when is the satellite in sunlight, and when is it in Earth's shadow? A satellite generates zero power during eclipse, so knowing exactly when those dark periods happen is the foundation of the entire power budget.
In this step, get ready to:
- Implement a sun position approximation using Julian Date math.
- Build an eclipse detector using the cylindrical shadow model.
- Verify that your ISS orbit spends roughly 35% of its time in shadow.
Compute the sun's direction in space
To detect when a satellite is in Earth's shadow, you first need to know where the sun is. The sun's position relative to Earth changes throughout the year as Earth orbits it.
- In main.py, find the end of your propagate_orbit() function.
- Add the following function below it:
# --- Sun Position (Simplified) ---
def compute_sun_direction(jd, fr):
"""Compute approximate sun direction vector in ECI frame."""
# Days since J2000.0 (2000 Jan 1, 12:00 TT)
T = (jd + fr) - 2451545.0
# Mean longitude and anomaly of the sun (degrees)
mean_longitude = (280.460 + 0.9856474 * T) % 360.0
mean_anomaly = (357.528 + 0.9856003 * T) % 360.0
mean_anomaly_rad = mean_anomaly * DEG2RAD
# Ecliptic longitude with correction terms
ecliptic_lon = mean_longitude + 1.915 * np.sin(mean_anomaly_rad) + 0.020 * np.sin(2 * mean_anomaly_rad)
ecliptic_lon_rad = ecliptic_lon * DEG2RAD
# Earth's axial tilt
obliquity_rad = 23.439 * DEG2RAD
# Sun direction unit vector in ECI (equatorial) frame
sun_x = np.cos(ecliptic_lon_rad)
sun_y = np.cos(obliquity_rad) * np.sin(ecliptic_lon_rad)
sun_z = np.sin(obliquity_rad) * np.sin(ecliptic_lon_rad)
magnitude = np.sqrt(sun_x**2 + sun_y**2 + sun_z**2)
return np.array([sun_x, sun_y, sun_z]) / magnitudeWhat does this code do?
This function computes a unit vector pointing from Earth toward the sun at any given moment in time.
- T counts days since the J2000 epoch (January 1, 2000). This is the standard time reference in orbital mechanics.
- The mean longitude and anomaly formulas model how fast the sun appears to move along the ecliptic (about 0.986 degrees per day, since 360/365.25 = 0.986).
- The obliquity (23.439 degrees) accounts for Earth's axial tilt when converting from ecliptic to equatorial coordinates.
- The output is a normalized direction vector in the ECI frame, which is the same coordinate system your satellite positions use.
What is the ECI frame?
ECI stands for Earth-Centered Inertial. It is a coordinate system where the origin is at Earth's center, the X-axis points toward the vernal equinox, and the Z-axis points toward the North Pole. Both SGP4 positions and this sun vector use this same frame, so you can directly compare them.
- Save main.py.
You will not see a visible change yet. The next function uses this sun direction to detect eclipses.
Implement the cylindrical shadow model
A satellite enters eclipse when Earth blocks the sun from reaching it. The cylindrical shadow model checks two conditions: is the satellite behind Earth relative to the sun, and is it within the cylinder of Earth's shadow?
- Below your compute_sun_direction() function, add the eclipse detection function:
# --- Eclipse Detection ---
def detect_eclipses(positions_km, jd_array, fr_array):
"""Detect eclipse periods using the cylindrical shadow model."""
num_steps = len(positions_km)
in_eclipse = np.zeros(num_steps, dtype=bool)
for i in range(num_steps):
sat_pos = positions_km[i]
sun_dir = compute_sun_direction(jd_array[i], fr_array[i])
# Project satellite position onto sun direction
projection = np.dot(sat_pos, sun_dir)
# Positive projection means satellite is on the sun-side of Earth
if projection >= 0:
in_eclipse[i] = False
continue
# Compute perpendicular distance from Earth-Sun line
perp_vector = sat_pos - projection * sun_dir
perp_distance = np.sqrt(np.dot(perp_vector, perp_vector))
# In eclipse if within Earth's shadow cylinder
in_eclipse[i] = perp_distance < EARTH_RADIUS_KM
return in_eclipseHow does the cylindrical shadow model work?
Imagine a cylinder extending from Earth's night side, stretching away from the sun. Any satellite inside this cylinder cannot see the sun.
- np.dot(sat_pos, sun_dir) projects the satellite's position onto the sun direction. A negative result means the satellite is on the opposite side of Earth from the sun.
- The perpendicular distance check measures how far the satellite is from the Earth-Sun line. If that distance is less than Earth's radius (6371 km), the satellite falls within the shadow cylinder.
- Both conditions must be true for eclipse: behind Earth AND within the cylinder's radius.
Why the cylindrical model instead of the conical model?
The full conical shadow model distinguishes between penumbra (partial shadow) and umbra (full shadow). For LEO power budget analysis, the cylindrical model produces results within a few percent of the conical model. Most CubeSat teams use this simpler approach in early mission planning.
- Save main.py.
You will not see a change yet. The next substep connects this function to your main simulation loop.
Run eclipse detection and print a summary
Now you will wire up the eclipse detector to your orbit data and print statistics that tell you whether the results make physical sense.
- In main.py, scroll down to your main() function.
- Find the altitude plotting code that starts with this line:
# Plot altitude
plt.figure(figsize=(12, 3))- Delete everything from that # Plot altitude comment down to and including plt.show().
- Replace it with the following eclipse detection code:
# Detect eclipse periods using cylindrical shadow model
print("\nDetecting eclipse periods...")
in_eclipse = detect_eclipses(positions_km, jd_array, fr_array)
eclipse_fraction = np.sum(in_eclipse) / len(in_eclipse)
print(f"Eclipse fraction: {eclipse_fraction*100:.1f}%")
print(f"Sunlight fraction: {(1-eclipse_fraction)*100:.1f}%")
# Count individual eclipse periods and compute average duration
eclipse_transitions = np.diff(in_eclipse.astype(int))
num_eclipses = np.sum(eclipse_transitions == 1)
print(f"Number of eclipse periods: {num_eclipses}")
if num_eclipses > 0:
avg_eclipse_min = (np.sum(in_eclipse) * TIME_STEP_SEC / 60) / num_eclipses
print(f"Average eclipse duration: {avg_eclipse_min:.1f} min")What does this code do?
- detect_eclipses() runs the shadow model at every timestep and returns a boolean array. True means eclipse, False means sunlight.
- np.diff(in_eclipse.astype(int)) finds transitions between sunlight and eclipse. A value of 1 in the diff means the satellite just entered eclipse.
- The average duration divides total eclipse time by the number of eclipse periods detected.
- Save main.py.
- Run the simulation by executing this command:
uv run python main.pyYou should see output showing an eclipse fraction of roughly 35% and 6 eclipse periods (one per orbit).
Why roughly 35%?
A typical LEO satellite at ISS altitude (~420 km) spends about 35% of each orbit behind Earth. This matches real-world data. If your number is significantly different (below 25% or above 45%), something is off in the eclipse geometry.
Seeing 0% eclipse or unexpected numbers?
- Make sure compute_sun_direction() is placed above detect_eclipses() in your file. Python reads top to bottom, so the function must be defined before it is called.
- Check that you are using EARTH_RADIUS_KM (not a hardcoded number) in the perpendicular distance comparison.
- Verify the projection check uses >= 0 to skip sun-side positions. A flipped sign would mark the wrong half of the orbit as eclipse.
Why is my eclipse detection returning unexpected results?
✔️ Awesome, I've got everything!
Great. Double check you have saved main.py before moving on.
ⓧ I'd like to double check the full code
Here is the complete main.py as it should look at the end of this step:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
# --- Configuration ---
EARTH_RADIUS_KM = 6371.0
SOLAR_CONSTANT = 1361.0 # W/m^2
DEG2RAD = np.pi / 180.0
# ISS TLE data (example - replace with current TLE for best accuracy)
TLE_LINE1 = "1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990"
TLE_LINE2 = "2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
# Simulation parameters
SIM_START = (2024, 2, 14, 12, 0, 0) # Year, month, day, hour, min, sec
NUM_ORBITS = 6
TIME_STEP_SEC = 30
# --- Orbit Propagation ---
def propagate_orbit(tle_line1, tle_line2, start_time, num_orbits, dt_sec):
"""Propagate satellite orbit using SGP4 and return positions over time."""
satellite = Satrec.twoline2rv(tle_line1, tle_line2)
orbital_period_min = 2 * np.pi / satellite.no_kozai
total_time_min = orbital_period_min * num_orbits
num_steps = int(total_time_min * 60 / dt_sec)
jd_start, fr_start = jday(*start_time)
times_min = np.zeros(num_steps)
positions_km = np.zeros((num_steps, 3))
altitudes_km = np.zeros(num_steps)
jd_array = np.zeros(num_steps)
fr_array = np.zeros(num_steps)
for i in range(num_steps):
t_min = i * dt_sec / 60.0
times_min[i] = t_min
fr = fr_start + (t_min / 1440.0)
jd = jd_start
while fr >= 1.0:
fr -= 1.0
jd += 1.0
jd_array[i] = jd
fr_array[i] = fr
e, r, v = satellite.sgp4(jd, fr)
if e != 0:
print(f"SGP4 error at step {i}: error code {e}")
continue
positions_km[i] = r
altitudes_km[i] = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) - EARTH_RADIUS_KM
return times_min, positions_km, altitudes_km, jd_array, fr_array
# --- Sun Position (Simplified) ---
def compute_sun_direction(jd, fr):
"""Compute approximate sun direction vector in ECI frame."""
T = (jd + fr) - 2451545.0
mean_longitude = (280.460 + 0.9856474 * T) % 360.0
mean_anomaly = (357.528 + 0.9856003 * T) % 360.0
mean_anomaly_rad = mean_anomaly * DEG2RAD
ecliptic_lon = mean_longitude + 1.915 * np.sin(mean_anomaly_rad) + 0.020 * np.sin(2 * mean_anomaly_rad)
ecliptic_lon_rad = ecliptic_lon * DEG2RAD
obliquity_rad = 23.439 * DEG2RAD
sun_x = np.cos(ecliptic_lon_rad)
sun_y = np.cos(obliquity_rad) * np.sin(ecliptic_lon_rad)
sun_z = np.sin(obliquity_rad) * np.sin(ecliptic_lon_rad)
magnitude = np.sqrt(sun_x**2 + sun_y**2 + sun_z**2)
return np.array([sun_x, sun_y, sun_z]) / magnitude
# --- Eclipse Detection ---
def detect_eclipses(positions_km, jd_array, fr_array):
"""Detect eclipse periods using the cylindrical shadow model."""
num_steps = len(positions_km)
in_eclipse = np.zeros(num_steps, dtype=bool)
for i in range(num_steps):
sat_pos = positions_km[i]
sun_dir = compute_sun_direction(jd_array[i], fr_array[i])
projection = np.dot(sat_pos, sun_dir)
if projection >= 0:
in_eclipse[i] = False
continue
perp_vector = sat_pos - projection * sun_dir
perp_distance = np.sqrt(np.dot(perp_vector, perp_vector))
in_eclipse[i] = perp_distance < EARTH_RADIUS_KM
return in_eclipse
def main():
print("Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f"Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f"Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
print("\nDetecting eclipse periods...")
in_eclipse = detect_eclipses(positions_km, jd_array, fr_array)
eclipse_fraction = np.sum(in_eclipse) / len(in_eclipse)
print(f"Eclipse fraction: {eclipse_fraction*100:.1f}%")
print(f"Sunlight fraction: {(1-eclipse_fraction)*100:.1f}%")
eclipse_transitions = np.diff(in_eclipse.astype(int))
num_eclipses = np.sum(eclipse_transitions == 1)
print(f"Number of eclipse periods: {num_eclipses}")
if num_eclipses > 0:
avg_eclipse_min = (np.sum(in_eclipse) * TIME_STEP_SEC / 60) / num_eclipses
print(f"Average eclipse duration: {avg_eclipse_min:.1f} min")
if __name__ == "__main__":
main()You now know exactly when your satellite is in sunlight and when it is hiding behind Earth. Next up, you will use this eclipse data to model how the power subsystem charges and discharges across each orbit.
Model the Power Subsystem
You now know exactly when your satellite enters and exits Earth's shadow. The next question is: what does that mean for the hardware?
In this step, you will define the CubeSat's power system parameters and simulate how solar panels charge the battery during sunlight, while the battery drains during eclipse. The result tells you whether your satellite can survive the dark periods without running out of power.
In this step, get ready to:
- Define the CubeSat power system configuration (panel area, battery capacity, load).
- Write a power simulation function that computes generation, consumption, and battery state of charge.
- Run the simulation and check whether the battery survives eclipse periods.
Define the CubeSat power system parameters
A CubeSat's power budget depends on a handful of physical parameters: how large the solar panels are, how efficient the cells are, how much energy the battery stores, and how much power the subsystems consume. These numbers come directly from the hardware spec sheet.
- In main.py, find the line DEG2RAD = np.pi / 180.0 near the top of the file.
- Add the following configuration dictionary below it (before the TLE data):
# CubeSat power system parameters
CUBESAT_CONFIG = {
"panel_area_m2": 0.06, # 3U body-mounted panels
"cell_efficiency": 0.30, # GaAs triple-junction
"battery_capacity_wh": 40.0, # Watt-hours
"initial_soc": 0.80, # 80% initial charge
"load_watts": 2.5, # Average subsystem load
"sun_incidence_deg": 23.0, # Average sun incidence angle
}What do these parameters mean?
These values model a typical 3U CubeSat with body-mounted solar panels.
- panel_area_m2: 0.06 m² is the exposed panel surface area on a 3U CubeSat body.
- cell_efficiency: 30% is the conversion efficiency of GaAs triple-junction solar cells, the standard for space applications.
- battery_capacity_wh: 40 Wh is a typical lithium-ion battery pack for a 3U CubeSat.
- load_watts: 2.5 W represents the average power draw from all onboard subsystems (computer, radio, sensors).
- sun_incidence_deg: 23° is the average angle between the panel normal and the incoming sunlight. Real missions compute this from attitude data, but the average works for first-pass analysis.
- Save main.py.
Write the power generation and battery simulation function
The power generation formula for solar panels in space is straightforward: P = area × efficiency × solar constant × cos(incidence angle). During eclipse, generation drops to zero and the battery must supply all power.
- In main.py, find the end of the detect_eclipses() function (the return in_eclipse line).
- Add the following function below it, before the def main(): line:
# --- Power Subsystem Model ---
def simulate_power_system(times_min, in_eclipse, config):
"""Simulate the CubeSat power system over the orbit."""
num_steps = len(times_min)
dt_hours = (times_min[1] - times_min[0]) / 60.0
power_generated = np.zeros(num_steps)
power_consumed = np.full(num_steps, config["load_watts"])
net_power = np.zeros(num_steps)
battery_soc = np.zeros(num_steps)
battery_soc[0] = config["initial_soc"]
sun_angle_rad = config["sun_incidence_deg"] * DEG2RAD
max_power = (config["panel_area_m2"] * config["cell_efficiency"]
* SOLAR_CONSTANT * np.cos(sun_angle_rad))
for i in range(num_steps):
# Power generation: full power in sunlight, zero in eclipse
if not in_eclipse[i]:
power_generated[i] = max_power
else:
power_generated[i] = 0.0
net_power[i] = power_generated[i] - power_consumed[i]
# Update battery state of charge
if i > 0:
delta_soc = (net_power[i] * dt_hours) / config["battery_capacity_wh"]
battery_soc[i] = np.clip(battery_soc[i - 1] + delta_soc, 0.0, 1.0)
return power_generated, power_consumed, net_power, battery_socWhat does this function do?
- max_power computes the peak solar panel output using the formula P = A × η × S × cos(θ). With these parameters, that is about 22.5 W.
- During sunlight, the panels generate at full capacity. During eclipse, generation is zero.
- net_power is the difference between generation and consumption. Positive means the battery is charging. Negative means it is draining.
- np.clip clamps the battery between 0% and 100%. A real battery management system enforces the same limits to protect the cells.
- Save main.py.
Seeing a syntax error?
- Make sure the max_power calculation is split across two lines with proper indentation. The second line (* SOLAR_CONSTANT * np.cos(sun_angle_rad))) must be indented to align under the opening parenthesis.
- Check that the function is defined at the top level (no extra indentation), not nested inside another function.
Why is my simulate_power_system function throwing an error?
Run the power simulation from main()
The simulation function is ready. Now you need to call it from main() and print the results to see whether your CubeSat's battery can handle the eclipse periods.
- In main.py, scroll down to the main() function.
- Find the last line inside main():
print(f"Average eclipse duration: {avg_eclipse_min:.1f} min")- Add the following code below that line (still inside main(), at the same indentation level):
print("\nSimulating power subsystem...")
power_gen, power_con, net_power, battery_soc = simulate_power_system(
times_min, in_eclipse, CUBESAT_CONFIG
)
max_gen = CUBESAT_CONFIG["panel_area_m2"] * CUBESAT_CONFIG["cell_efficiency"] * SOLAR_CONSTANT * np.cos(CUBESAT_CONFIG["sun_incidence_deg"] * DEG2RAD)
print(f"Max solar power: {max_gen:.2f} W")
print(f"Constant load: {CUBESAT_CONFIG['load_watts']:.1f} W")
print(f"Average net power: {np.mean(net_power):.2f} W")
print(f"Minimum battery SoC: {np.min(battery_soc)*100:.1f}%")What does this code do?
- It calls simulate_power_system() with the time array, eclipse array, and hardware config you defined earlier.
- max_gen recalculates peak solar power so you can print it as a reference value.
- The minimum battery SoC is the most critical metric. If it drops below 20%, the mission is at risk of losing subsystems.
- Save main.py.
- Run the simulation by executing this command:
uv run python main.py- You should see the orbit and eclipse output from before, followed by new power subsystem results showing max solar power around 22.5 W, the 2.5 W constant load, average net power, and a minimum battery SoC percentage.
What should the numbers look like?
For the ISS orbit with these parameters, you should see the minimum battery SoC stay well above 20%. The average net power should be positive, meaning the panels generate more energy over one orbit than the subsystems consume.
If minimum SoC is suspiciously low or zero, double-check that battery_soc[0] is set to config["initial_soc"] (0.80) before the loop starts.
Getting a KeyError or NameError?
- A NameError: name 'CUBESAT_CONFIG' is not defined means the dictionary is not at the top level of the file. Make sure it is defined between DEG2RAD and the TLE lines, with no extra indentation.
- A KeyError means a key name in the dictionary does not match what the function expects. Check that all keys use underscores and match exactly (e.g. panel_area_m2, not panel_area).
Help me debug my power subsystem simulation
✔️ Awesome, I've got everything!
Great. Double check you have saved main.py and your terminal output shows the power subsystem results.
ⓧ I'd like to double check the full code
Here is the full main.py as it should look at the end of this step:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
# --- Configuration ---
EARTH_RADIUS_KM = 6371.0
SOLAR_CONSTANT = 1361.0 # W/m^2
DEG2RAD = np.pi / 180.0
# CubeSat power system parameters
CUBESAT_CONFIG = {
"panel_area_m2": 0.06, # 3U body-mounted panels
"cell_efficiency": 0.30, # GaAs triple-junction
"battery_capacity_wh": 40.0, # Watt-hours
"initial_soc": 0.80, # 80% initial charge
"load_watts": 2.5, # Average subsystem load
"sun_incidence_deg": 23.0, # Average sun incidence angle
}
# ISS TLE data (example - replace with current TLE for best accuracy)
TLE_LINE1 = "1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990"
TLE_LINE2 = "2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
# Simulation parameters
SIM_START = (2024, 2, 14, 12, 0, 0)
NUM_ORBITS = 6
TIME_STEP_SEC = 30
# --- Orbit Propagation ---
def propagate_orbit(tle_line1, tle_line2, start_time, num_orbits, dt_sec):
"""Propagate satellite orbit using SGP4 and return positions over time."""
satellite = Satrec.twoline2rv(tle_line1, tle_line2)
orbital_period_min = 2 * np.pi / satellite.no_kozai
total_time_min = orbital_period_min * num_orbits
num_steps = int(total_time_min * 60 / dt_sec)
jd_start, fr_start = jday(*start_time)
times_min = np.zeros(num_steps)
positions_km = np.zeros((num_steps, 3))
altitudes_km = np.zeros(num_steps)
jd_array = np.zeros(num_steps)
fr_array = np.zeros(num_steps)
for i in range(num_steps):
t_min = i * dt_sec / 60.0
times_min[i] = t_min
fr = fr_start + (t_min / 1440.0)
jd = jd_start
while fr >= 1.0:
fr -= 1.0
jd += 1.0
jd_array[i] = jd
fr_array[i] = fr
e, r, v = satellite.sgp4(jd, fr)
if e != 0:
print(f"SGP4 error at step {i}: error code {e}")
continue
positions_km[i] = r
altitudes_km[i] = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) - EARTH_RADIUS_KM
return times_min, positions_km, altitudes_km, jd_array, fr_array
# --- Sun Position (Simplified) ---
def compute_sun_direction(jd, fr):
"""Compute approximate sun direction vector in ECI frame."""
T = (jd + fr) - 2451545.0
mean_longitude = (280.460 + 0.9856474 * T) % 360.0
mean_anomaly = (357.528 + 0.9856003 * T) % 360.0
mean_anomaly_rad = mean_anomaly * DEG2RAD
ecliptic_lon = mean_longitude + 1.915 * np.sin(mean_anomaly_rad) + 0.020 * np.sin(2 * mean_anomaly_rad)
ecliptic_lon_rad = ecliptic_lon * DEG2RAD
obliquity_rad = 23.439 * DEG2RAD
sun_x = np.cos(ecliptic_lon_rad)
sun_y = np.cos(obliquity_rad) * np.sin(ecliptic_lon_rad)
sun_z = np.sin(obliquity_rad) * np.sin(ecliptic_lon_rad)
magnitude = np.sqrt(sun_x**2 + sun_y**2 + sun_z**2)
return np.array([sun_x, sun_y, sun_z]) / magnitude
# --- Eclipse Detection ---
def detect_eclipses(positions_km, jd_array, fr_array):
"""Detect eclipse periods using the cylindrical shadow model."""
num_steps = len(positions_km)
in_eclipse = np.zeros(num_steps, dtype=bool)
for i in range(num_steps):
sat_pos = positions_km[i]
sun_dir = compute_sun_direction(jd_array[i], fr_array[i])
projection = np.dot(sat_pos, sun_dir)
if projection >= 0:
in_eclipse[i] = False
continue
perp_vector = sat_pos - projection * sun_dir
perp_distance = np.sqrt(np.dot(perp_vector, perp_vector))
in_eclipse[i] = perp_distance < EARTH_RADIUS_KM
return in_eclipse
# --- Power Subsystem Model ---
def simulate_power_system(times_min, in_eclipse, config):
"""Simulate the CubeSat power system over the orbit."""
num_steps = len(times_min)
dt_hours = (times_min[1] - times_min[0]) / 60.0
power_generated = np.zeros(num_steps)
power_consumed = np.full(num_steps, config["load_watts"])
net_power = np.zeros(num_steps)
battery_soc = np.zeros(num_steps)
battery_soc[0] = config["initial_soc"]
sun_angle_rad = config["sun_incidence_deg"] * DEG2RAD
max_power = (config["panel_area_m2"] * config["cell_efficiency"]
* SOLAR_CONSTANT * np.cos(sun_angle_rad))
for i in range(num_steps):
if not in_eclipse[i]:
power_generated[i] = max_power
else:
power_generated[i] = 0.0
net_power[i] = power_generated[i] - power_consumed[i]
if i > 0:
delta_soc = (net_power[i] * dt_hours) / config["battery_capacity_wh"]
battery_soc[i] = np.clip(battery_soc[i - 1] + delta_soc, 0.0, 1.0)
return power_generated, power_consumed, net_power, battery_soc
def main():
print("Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f"Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f"Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
print("\nDetecting eclipse periods...")
in_eclipse = detect_eclipses(positions_km, jd_array, fr_array)
eclipse_fraction = np.sum(in_eclipse) / len(in_eclipse)
print(f"Eclipse fraction: {eclipse_fraction*100:.1f}%")
print(f"Sunlight fraction: {(1-eclipse_fraction)*100:.1f}%")
eclipse_transitions = np.diff(in_eclipse.astype(int))
num_eclipses = np.sum(eclipse_transitions == 1)
print(f"Number of eclipse periods: {num_eclipses}")
if num_eclipses > 0:
avg_eclipse_min = (np.sum(in_eclipse) * TIME_STEP_SEC / 60) / num_eclipses
print(f"Average eclipse duration: {avg_eclipse_min:.1f} min")
print("\nSimulating power subsystem...")
power_gen, power_con, net_power, battery_soc = simulate_power_system(
times_min, in_eclipse, CUBESAT_CONFIG
)
max_gen = CUBESAT_CONFIG["panel_area_m2"] * CUBESAT_CONFIG["cell_efficiency"] * SOLAR_CONSTANT * np.cos(CUBESAT_CONFIG["sun_incidence_deg"] * DEG2RAD)
print(f"Max solar power: {max_gen:.2f} W")
print(f"Constant load: {CUBESAT_CONFIG['load_watts']:.1f} W")
print(f"Average net power: {np.mean(net_power):.2f} W")
print(f"Minimum battery SoC: {np.min(battery_soc)*100:.1f}%")
if __name__ == "__main__":
main()Your power subsystem model is complete. You now know the satellite's peak generation, its average power margin, and whether the battery survives eclipse. Next up, you will visualize all of this data in a 4-panel plot that shows orbit, eclipses, power, and battery state together.
Visualize the Complete Simulation
You have a working orbit propagator, eclipse detector, and power model. The simulation produces arrays of numbers, but numbers alone don't tell you whether your CubeSat survives its mission.
In this step, you will build a 4-panel matplotlib visualization that shows how orbital altitude, eclipse windows, power generation, and battery charge relate to each other over time. You will also add a mission feasibility check that tells you whether the battery stays above the critical 20% threshold.
In this step, get ready to:
- Build a 4-panel plot showing altitude, eclipse state, power budget, and battery state of charge.
- Add annotations marking critical battery levels and mission parameters.
- Run the end-to-end simulation and generate a mission feasibility report.
Create the visualization function
A 4-panel figure with a shared x-axis lets you see exactly how orbital position drives power system behavior. Each panel tells a different part of the story, and gray eclipse shading across all panels makes the cause-and-effect relationship visually obvious.
- In main.py, find the end of the simulate_power_system() function.
- Add the following create_visualization() function below it, before the main() function:
# --- Visualization ---
def create_visualization(times_min, altitudes_km, in_eclipse,
power_generated, power_consumed, battery_soc, config):
"""Create a 4-panel visualization of the simulation results."""
times_hours = times_min / 60.0
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)
fig.suptitle(
f"CubeSat Power Simulation | Panel: {config['panel_area_m2']} m\u00b2 | "
f"Battery: {config['battery_capacity_wh']} Wh | Load: {config['load_watts']} W",
fontsize=12, fontweight="bold"
)
# Panel 1: Altitude
axes[0].plot(times_hours, altitudes_km, "b-", linewidth=0.8)
axes[0].set_ylabel("Altitude (km)")
axes[0].set_title("Orbital Altitude")
axes[0].grid(True, alpha=0.3)
# Add eclipse shading to all panels
for ax in axes:
eclipse_starts = []
eclipse_ends = []
in_ecl = False
for i in range(len(in_eclipse)):
if in_eclipse[i] and not in_ecl:
eclipse_starts.append(times_hours[i])
in_ecl = True
elif not in_eclipse[i] and in_ecl:
eclipse_ends.append(times_hours[i])
in_ecl = False
if in_ecl:
eclipse_ends.append(times_hours[-1])
for start, end in zip(eclipse_starts, eclipse_ends):
ax.axvspan(start, end, alpha=0.15, color="gray")
# Panel 2: Eclipse state
axes[1].fill_between(times_hours, 0, in_eclipse.astype(int),
color="darkgray", alpha=0.6, label="Eclipse")
axes[1].fill_between(times_hours, 0, (~in_eclipse).astype(int),
color="gold", alpha=0.4, label="Sunlight")
axes[1].set_ylabel("State")
axes[1].set_title("Sunlight / Eclipse")
axes[1].set_yticks([0, 1])
axes[1].set_yticklabels(["Eclipse", "Sunlight"])
axes[1].legend(loc="upper right")
axes[1].grid(True, alpha=0.3)
# Panel 3: Power
axes[2].plot(times_hours, power_generated, "orange", linewidth=0.8, label="Generated")
axes[2].plot(times_hours, power_consumed, "r--", linewidth=0.8, label="Consumed")
axes[2].set_ylabel("Power (W)")
axes[2].set_title("Power Budget")
axes[2].legend(loc="upper right")
axes[2].grid(True, alpha=0.3)
# Panel 4: Battery SoC
axes[3].plot(times_hours, battery_soc * 100, "g-", linewidth=1.2)
axes[3].axhline(y=20, color="red", linestyle="--", linewidth=0.8, label="Critical (20%)")
axes[3].set_ylabel("SoC (%)")
axes[3].set_xlabel("Time (hours)")
axes[3].set_title("Battery State of Charge")
axes[3].set_ylim(0, 105)
axes[3].legend(loc="upper right")
axes[3].grid(True, alpha=0.3)
# Mark minimum SoC
min_soc_idx = np.argmin(battery_soc)
axes[3].annotate(
f"Min: {battery_soc[min_soc_idx]*100:.1f}%",
xy=(times_hours[min_soc_idx], battery_soc[min_soc_idx] * 100),
xytext=(times_hours[min_soc_idx] + 0.3, battery_soc[min_soc_idx] * 100 + 10),
arrowprops=dict(arrowstyle="->", color="red"),
fontsize=9, color="red"
)
plt.tight_layout()
plt.savefig("cubesat_power_sim.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nPlot saved to cubesat_power_sim.png")What does this code do?
This function creates a stacked 4-panel figure where each panel shares the same x-axis (time in hours). Here is what each section handles:
- Panel 1 plots orbital altitude in km, showing the sinusoidal pattern as the satellite moves between perigee and apogee.
- The eclipse shading loop adds gray vertical spans to ALL four panels wherever the satellite is in shadow. This makes it immediately obvious how eclipse drives the power system.
- Panel 2 uses fill_between() to show gold for sunlight and dark gray for eclipse as a binary state.
- Panel 3 overlays generated power (orange) against consumed power (red dashed). When the orange line drops to zero, the satellite is in eclipse.
- Panel 4 shows battery state of charge with a red dashed danger line at 20%. The annotate() call marks the lowest point the battery reaches during the entire simulation.
- Save main.py.
You won't see a change yet. The next chunk updates main() to actually call this function.
Update main() with formatted output and feasibility status
The current main() function prints basic results but does not call the visualization or make a mission feasibility decision. You need to replace it with a version that formats the output as a structured report and calls create_visualization().
- In main.py, find the existing main() function. It starts with:
def main():
print("Propagating orbit with SGP4...")- Delete the entire main() function (everything from def main(): down to the last print statement before if __name__ == "__main__":).
- Replace it with the following updated version:
# --- Main ---
def main():
print("=" * 60)
print(" CubeSat Power System Simulation")
print("=" * 60)
# Step 1: Propagate orbit
print("\n[1/4] Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f" Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f" Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
# Step 2: Detect eclipses
print("\n[2/4] Detecting eclipse periods...")
in_eclipse = detect_eclipses(positions_km, jd_array, fr_array)
eclipse_fraction = np.sum(in_eclipse) / len(in_eclipse)
print(f" Eclipse fraction: {eclipse_fraction*100:.1f}%")
print(f" Sunlight fraction: {(1-eclipse_fraction)*100:.1f}%")
# Count eclipse periods
eclipse_transitions = np.diff(in_eclipse.astype(int))
num_eclipses = np.sum(eclipse_transitions == 1)
print(f" Number of eclipse periods: {num_eclipses}")
if num_eclipses > 0:
avg_eclipse_min = (np.sum(in_eclipse) * TIME_STEP_SEC / 60) / num_eclipses
print(f" Average eclipse duration: {avg_eclipse_min:.1f} min")
# Step 3: Simulate power system
print("\n[3/4] Simulating power subsystem...")
power_gen, power_con, net_power, battery_soc = simulate_power_system(
times_min, in_eclipse, CUBESAT_CONFIG
)
max_gen = CUBESAT_CONFIG["panel_area_m2"] * CUBESAT_CONFIG["cell_efficiency"] * SOLAR_CONSTANT * np.cos(CUBESAT_CONFIG["sun_incidence_deg"] * DEG2RAD)
print(f" Max solar power: {max_gen:.2f} W")
print(f" Constant load: {CUBESAT_CONFIG['load_watts']:.1f} W")
print(f" Average net power: {np.mean(net_power):.2f} W")
# Step 4: Mission feasibility summary
print("\n[4/4] Mission Feasibility Summary")
print("-" * 40)
min_soc = np.min(battery_soc)
print(f" Minimum battery SoC: {min_soc*100:.1f}%")
print(f" Average power margin: {np.mean(net_power):.2f} W")
print(f" Eclipse time: {eclipse_fraction*100:.1f}%")
if min_soc > 0.20:
print(" STATUS: MISSION FEASIBLE")
print(" Battery never drops below 20% critical threshold.")
else:
print(" STATUS: MISSION AT RISK")
print(" Battery drops below 20%! Consider larger panels or battery.")
# Generate visualization
print("\n Generating visualization...")
create_visualization(
times_min, altitudes_km, in_eclipse,
power_gen, power_con, battery_soc, CUBESAT_CONFIG
)What does this code do?
The updated main() now runs all four simulation stages in sequence with clear progress indicators.
- The [1/4] through [4/4] labels give you a progress report as the simulation runs.
- The feasibility check compares minimum battery SoC against the 20% threshold. If the battery stays above 20%, the mission is feasible. If not, you would need to increase panel area or battery capacity.
- The final call to create_visualization() generates and saves the 4-panel plot.
- Make sure the if __name__ == "__main__": block is still at the very bottom of the file, unchanged.
- Save main.py.
Run the complete simulation end-to-end
Everything is in place. Time to run the full simulation and see the results.
- Run the simulation by executing this command in your terminal:
uv run python main.pyYou should see a structured report in your terminal with all four stages completing. The output ends with a Mission Feasibility Summary showing minimum battery SoC, average power margin, eclipse time percentage, and a STATUS line reading either MISSION FEASIBLE or MISSION AT RISK.
A matplotlib window opens showing the 4-panel plot. Gray shaded regions mark eclipse periods across all panels. The bottom panel shows battery charge rising during sunlight and falling during eclipse, with a red dashed line at 20% and a red arrow marking the minimum SoC point.
- Close the matplotlib window after reviewing the plot.
- Confirm that cubesat_power_sim.png was saved in your project folder by running:
ls cubesat_power_sim.pngYou should see cubesat_power_sim.png listed in the output.
Plot window not appearing?
- If you are running on a remote server or inside a container, matplotlib may not have a display backend available. The plot still saves to the PNG file. You can open it directly from your file browser.
- If you get a NameError: name 'create_visualization' is not defined, make sure the function is defined BEFORE main() in the file.
- If the plot looks empty or flat, check that your main() passes all the correct variables to create_visualization().
help me debug my CubeSat visualization
What should the plot tell you?
For the ISS orbit at ~420 km altitude, you should see approximately 35% of the time spent in eclipse. The battery SoC dips during each eclipse period and recharges during sunlight. With 0.06 m² of panel area and a 2.5 W load, the battery should stay well above the 20% danger line.
This is the exact analysis real CubeSat teams run before committing to hardware. If the battery dropped below 20%, you would need larger panels, a bigger battery, or a lower power budget.
✔️ Awesome, I've got everything!
Great. Double check you have saved main.py and confirmed the plot saved successfully.
ⓧ I'd like to double check the full code
Here is the complete main.py file as it should look at the end of this step:
import numpy as np
import matplotlib.pyplot as plt
from sgp4.api import Satrec, jday
# --- Configuration ---
EARTH_RADIUS_KM = 6371.0
SOLAR_CONSTANT = 1361.0 # W/m^2
DEG2RAD = np.pi / 180.0
# CubeSat power system parameters
CUBESAT_CONFIG = {
"panel_area_m2": 0.06, # 3U body-mounted panels
"cell_efficiency": 0.30, # GaAs triple-junction
"battery_capacity_wh": 40.0, # Watt-hours
"initial_soc": 0.80, # 80% initial charge
"load_watts": 2.5, # Average subsystem load
"sun_incidence_deg": 23.0, # Average sun incidence angle
}
# ISS TLE data (example - replace with current TLE for best accuracy)
TLE_LINE1 = "1 25544U 98067A 24045.51780806 .00015925 00000-0 28461-3 0 9990"
TLE_LINE2 = "2 25544 51.6414 263.3740 0005650 108.6591 320.8797 15.49579897441075"
# Simulation parameters
SIM_START = (2024, 2, 14, 12, 0, 0) # Year, month, day, hour, min, sec
NUM_ORBITS = 6
TIME_STEP_SEC = 30
# --- Orbit Propagation ---
def propagate_orbit(tle_line1, tle_line2, start_time, num_orbits, dt_sec):
"""Propagate satellite orbit using SGP4 and return positions over time."""
satellite = Satrec.twoline2rv(tle_line1, tle_line2)
# Compute orbital period from mean motion (radians/minute)
orbital_period_min = 2 * np.pi / satellite.no_kozai
total_time_min = orbital_period_min * num_orbits
num_steps = int(total_time_min * 60 / dt_sec)
# Start time as Julian Date
jd_start, fr_start = jday(*start_time)
times_min = np.zeros(num_steps)
positions_km = np.zeros((num_steps, 3))
altitudes_km = np.zeros(num_steps)
jd_array = np.zeros(num_steps)
fr_array = np.zeros(num_steps)
for i in range(num_steps):
t_min = i * dt_sec / 60.0
times_min[i] = t_min
# Compute Julian Date for this timestep
fr = fr_start + (t_min / 1440.0) # Convert minutes to fraction of day
jd = jd_start
# Handle day overflow
while fr >= 1.0:
fr -= 1.0
jd += 1.0
jd_array[i] = jd
fr_array[i] = fr
e, r, v = satellite.sgp4(jd, fr)
if e != 0:
print(f"SGP4 error at step {i}: error code {e}")
continue
positions_km[i] = r
altitudes_km[i] = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) - EARTH_RADIUS_KM
return times_min, positions_km, altitudes_km, jd_array, fr_array
# --- Sun Position (Simplified) ---
def compute_sun_direction(jd, fr):
"""Compute approximate sun direction vector in ECI frame.
Uses a simplified model: the sun moves along the ecliptic
at ~360 degrees per 365.25 days. This is accurate enough
for eclipse detection purposes.
"""
# Days since J2000.0 (2000 Jan 1, 12:00 TT)
T = (jd + fr) - 2451545.0
# Mean longitude of the sun (degrees)
mean_longitude = (280.460 + 0.9856474 * T) % 360.0
# Mean anomaly (degrees)
mean_anomaly = (357.528 + 0.9856003 * T) % 360.0
mean_anomaly_rad = mean_anomaly * DEG2RAD
# Ecliptic longitude (degrees)
ecliptic_lon = mean_longitude + 1.915 * np.sin(mean_anomaly_rad) + 0.020 * np.sin(2 * mean_anomaly_rad)
ecliptic_lon_rad = ecliptic_lon * DEG2RAD
# Obliquity of ecliptic
obliquity_rad = 23.439 * DEG2RAD
# Sun direction unit vector in ECI (equatorial) frame
sun_x = np.cos(ecliptic_lon_rad)
sun_y = np.cos(obliquity_rad) * np.sin(ecliptic_lon_rad)
sun_z = np.sin(obliquity_rad) * np.sin(ecliptic_lon_rad)
magnitude = np.sqrt(sun_x**2 + sun_y**2 + sun_z**2)
return np.array([sun_x, sun_y, sun_z]) / magnitude
# --- Eclipse Detection ---
def detect_eclipses(positions_km, jd_array, fr_array):
"""Detect eclipse periods using the cylindrical shadow model.
A satellite is in eclipse when:
1. It is on the opposite side of Earth from the sun (projection < 0)
2. Its perpendicular distance from the Earth-Sun line is less than Earth's radius
"""
num_steps = len(positions_km)
in_eclipse = np.zeros(num_steps, dtype=bool)
for i in range(num_steps):
sat_pos = positions_km[i]
sun_dir = compute_sun_direction(jd_array[i], fr_array[i])
# Project satellite position onto sun direction
projection = np.dot(sat_pos, sun_dir)
# If projection is positive, satellite is on the sun-side of Earth
if projection >= 0:
in_eclipse[i] = False
continue
# Compute perpendicular distance from Earth-Sun line
perp_vector = sat_pos - projection * sun_dir
perp_distance = np.sqrt(np.dot(perp_vector, perp_vector))
# In eclipse if within Earth's shadow cylinder
in_eclipse[i] = perp_distance < EARTH_RADIUS_KM
return in_eclipse
# --- Power Subsystem Model ---
def simulate_power_system(times_min, in_eclipse, config):
"""Simulate the CubeSat power system over the orbit."""
num_steps = len(times_min)
dt_hours = (times_min[1] - times_min[0]) / 60.0 # Time step in hours
power_generated = np.zeros(num_steps)
power_consumed = np.full(num_steps, config["load_watts"])
net_power = np.zeros(num_steps)
battery_soc = np.zeros(num_steps)
battery_soc[0] = config["initial_soc"]
sun_angle_rad = config["sun_incidence_deg"] * DEG2RAD
max_power = (config["panel_area_m2"] * config["cell_efficiency"]
* SOLAR_CONSTANT * np.cos(sun_angle_rad))
for i in range(num_steps):
# Power generation: full power in sunlight, zero in eclipse
if not in_eclipse[i]:
power_generated[i] = max_power
else:
power_generated[i] = 0.0
net_power[i] = power_generated[i] - power_consumed[i]
# Update battery SoC
if i > 0:
delta_soc = (net_power[i] * dt_hours) / config["battery_capacity_wh"]
battery_soc[i] = np.clip(battery_soc[i - 1] + delta_soc, 0.0, 1.0)
return power_generated, power_consumed, net_power, battery_soc
# --- Visualization ---
def create_visualization(times_min, altitudes_km, in_eclipse,
power_generated, power_consumed, battery_soc, config):
"""Create a 4-panel visualization of the simulation results."""
times_hours = times_min / 60.0
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)
fig.suptitle(
f"CubeSat Power Simulation | Panel: {config['panel_area_m2']} m\u00b2 | "
f"Battery: {config['battery_capacity_wh']} Wh | Load: {config['load_watts']} W",
fontsize=12, fontweight="bold"
)
# Panel 1: Altitude
axes[0].plot(times_hours, altitudes_km, "b-", linewidth=0.8)
axes[0].set_ylabel("Altitude (km)")
axes[0].set_title("Orbital Altitude")
axes[0].grid(True, alpha=0.3)
# Add eclipse shading to all panels
for ax in axes:
eclipse_starts = []
eclipse_ends = []
in_ecl = False
for i in range(len(in_eclipse)):
if in_eclipse[i] and not in_ecl:
eclipse_starts.append(times_hours[i])
in_ecl = True
elif not in_eclipse[i] and in_ecl:
eclipse_ends.append(times_hours[i])
in_ecl = False
if in_ecl:
eclipse_ends.append(times_hours[-1])
for start, end in zip(eclipse_starts, eclipse_ends):
ax.axvspan(start, end, alpha=0.15, color="gray")
# Panel 2: Eclipse state
axes[1].fill_between(times_hours, 0, in_eclipse.astype(int),
color="darkgray", alpha=0.6, label="Eclipse")
axes[1].fill_between(times_hours, 0, (~in_eclipse).astype(int),
color="gold", alpha=0.4, label="Sunlight")
axes[1].set_ylabel("State")
axes[1].set_title("Sunlight / Eclipse")
axes[1].set_yticks([0, 1])
axes[1].set_yticklabels(["Eclipse", "Sunlight"])
axes[1].legend(loc="upper right")
axes[1].grid(True, alpha=0.3)
# Panel 3: Power
axes[2].plot(times_hours, power_generated, "orange", linewidth=0.8, label="Generated")
axes[2].plot(times_hours, power_consumed, "r--", linewidth=0.8, label="Consumed")
axes[2].set_ylabel("Power (W)")
axes[2].set_title("Power Budget")
axes[2].legend(loc="upper right")
axes[2].grid(True, alpha=0.3)
# Panel 4: Battery SoC
axes[3].plot(times_hours, battery_soc * 100, "g-", linewidth=1.2)
axes[3].axhline(y=20, color="red", linestyle="--", linewidth=0.8, label="Critical (20%)")
axes[3].set_ylabel("SoC (%)")
axes[3].set_xlabel("Time (hours)")
axes[3].set_title("Battery State of Charge")
axes[3].set_ylim(0, 105)
axes[3].legend(loc="upper right")
axes[3].grid(True, alpha=0.3)
# Mark minimum SoC
min_soc_idx = np.argmin(battery_soc)
axes[3].annotate(
f"Min: {battery_soc[min_soc_idx]*100:.1f}%",
xy=(times_hours[min_soc_idx], battery_soc[min_soc_idx] * 100),
xytext=(times_hours[min_soc_idx] + 0.3, battery_soc[min_soc_idx] * 100 + 10),
arrowprops=dict(arrowstyle="->", color="red"),
fontsize=9, color="red"
)
plt.tight_layout()
plt.savefig("cubesat_power_sim.png", dpi=150, bbox_inches="tight")
plt.show()
print("\nPlot saved to cubesat_power_sim.png")
# --- Main ---
def main():
print("=" * 60)
print(" CubeSat Power System Simulation")
print("=" * 60)
# Step 1: Propagate orbit
print("\n[1/4] Propagating orbit with SGP4...")
times_min, positions_km, altitudes_km, jd_array, fr_array = propagate_orbit(
TLE_LINE1, TLE_LINE2, SIM_START, NUM_ORBITS, TIME_STEP_SEC
)
print(f" Simulated {NUM_ORBITS} orbits ({times_min[-1]/60:.1f} hours)")
print(f" Altitude range: {altitudes_km.min():.1f} - {altitudes_km.max():.1f} km")
# Step 2: Detect eclipses
print("\n[2/4] Detecting eclipse periods...")
in_eclipse = detect_eclipses(positions_km, jd_array, fr_array)
eclipse_fraction = np.sum(in_eclipse) / len(in_eclipse)
print(f" Eclipse fraction: {eclipse_fraction*100:.1f}%")
print(f" Sunlight fraction: {(1-eclipse_fraction)*100:.1f}%")
# Count eclipse periods
eclipse_transitions = np.diff(in_eclipse.astype(int))
num_eclipses = np.sum(eclipse_transitions == 1)
print(f" Number of eclipse periods: {num_eclipses}")
if num_eclipses > 0:
avg_eclipse_min = (np.sum(in_eclipse) * TIME_STEP_SEC / 60) / num_eclipses
print(f" Average eclipse duration: {avg_eclipse_min:.1f} min")
# Step 3: Simulate power system
print("\n[3/4] Simulating power subsystem...")
power_gen, power_con, net_power, battery_soc = simulate_power_system(
times_min, in_eclipse, CUBESAT_CONFIG
)
max_gen = CUBESAT_CONFIG["panel_area_m2"] * CUBESAT_CONFIG["cell_efficiency"] * SOLAR_CONSTANT * np.cos(CUBESAT_CONFIG["sun_incidence_deg"] * DEG2RAD)
print(f" Max solar power: {max_gen:.2f} W")
print(f" Constant load: {CUBESAT_CONFIG['load_watts']:.1f} W")
print(f" Average net power: {np.mean(net_power):.2f} W")
# Step 4: Mission feasibility summary
print("\n[4/4] Mission Feasibility Summary")
print("-" * 40)
min_soc = np.min(battery_soc)
print(f" Minimum battery SoC: {min_soc*100:.1f}%")
print(f" Average power margin: {np.mean(net_power):.2f} W")
print(f" Eclipse time: {eclipse_fraction*100:.1f}%")
if min_soc > 0.20:
print(" STATUS: MISSION FEASIBLE")
print(" Battery never drops below 20% critical threshold.")
else:
print(" STATUS: MISSION AT RISK")
print(" Battery drops below 20%! Consider larger panels or battery.")
# Generate visualization
print("\n Generating visualization...")
create_visualization(
times_min, altitudes_km, in_eclipse,
power_gen, power_con, battery_soc, CUBESAT_CONFIG
)
if __name__ == "__main__":
main()Your simulation is complete. You have built the same kind of power budget analysis that real CubeSat teams run before committing to hardware. Next up, you can extend this simulation with intelligent load scheduling in the Secret Mission.
Secret mission
Add Mission Phase Load Scheduling
Real CubeSats don't consume power at a constant rate. They cycle through mission phases like imaging, downlink, and idle, each with wildly different power draws. Can you upgrade your simulation to model intelligent load scheduling and autonomous power protection?
Clean Up Your Resources
Clean Up Your Resources
This project runs entirely locally on your machine, so there are no ongoing costs or cloud resources to worry about.
Decide whether to keep your project files for future experimentation, or delete them to free up disk space.
Resources you used:
- cubesat-power-sim project directory (contains main.py, pyproject.toml, and configuration files).
- Python virtual environment (.venv folder managed by uv).
- Installed packages: sgp4, numpy, matplotlib.
- Generated output file: cubesat_power_sim.png.
Keep everything running
No action needed. Choose this if you want to keep experimenting with different orbital parameters, panel sizes, or battery configurations.
- Your cubesat-power-sim folder stays on your Desktop with all code and output intact.
- You can re-run the simulation anytime with uv run python main.py.
Pause - I'll come back to this later
Since this is a local project with no running processes, there is nothing to pause. Your files sit idle on disk and consume no CPU or memory when not running.
- Close your code editor and any terminal windows running the simulation.
- Your project is ready to pick up again whenever you return.
Delete - I don't want to use this again
Remove all project files and the virtual environment to free up disk space.
- Close your code editor and any open terminal windows.
- Delete the entire cubesat-power-sim folder from your Desktop.
- On macOS, drag the folder to the Trash. On Windows, right-click and select Delete.
- This removes all code, the virtual environment, installed packages, and generated output.
Nice Work!
Nice Work!
You just built the same kind of simulation that real space engineers run before every CubeSat mission. Your code propagates a real satellite orbit, detects eclipse periods, and proves whether a power system can survive in space.
You've learned how to:
- Built an orbit propagator using the same SGP4 algorithm that NASA and NORAD use to track every object in orbit.
- Implemented eclipse detection using the cylindrical shadow model and modeled a complete power subsystem with solar generation, battery storage, and load management.
- Created a 4-panel visualization that proves whether your CubeSat's battery can survive eclipse periods without hitting critical levels.
- Added intelligent load scheduling with mission phases that mirrors real flight software power management.
Ready to quiz yourself?