You are not logged in.
This post is about calculation of horizontal velocity needed to orbitize a payload launched by a vertical launcher.
Dr. John Hunter is the leading developer of such a system, and there may well be others. Dr. Hunter's work has been reported at some length in the forum.
My interest in requesting this program, was to learn how the amount of horizontal velocity changes with altitude. The effect we see in the output of the program may be due to the inverse square law. The decrease in velocity changes dramatically from the program starting point of 100 km above the surface of Earth, to a nearly constant value at GEO.
In practice, it appears that solid fuel rockets may be the most practical solution for orbitizing the payloads, but some (kbd512 in particular) advocate investigation of liquid systems for the purpose.
Anyone should be able to run this program if you are using Linux, Windows or Apple. Python can run under Chrome OS using the Linux addon.
import math G = 6.6743e-11 # gravitational constant in m^3/kg/s² M = 5.972e24 # mass of Earth in kg R_earth = 6371e3 # Earth's radius in meters (6,371 km) geosynchronous_distance = R_earth + 35786e3 # geosynchronous orbit distance from Earth's center altitudes = [] velocities = [] current_altitude_km = 100 # starting altitude is 100 km above Earth's surface max_altitude_km = 35786 # maximum altitude (geosynchronous orbit) while current_altitude_km <= max_altitude_km: distance_from_center = R_earth + current_altitude_km * 1e3 # convert to meters velocity = math.sqrt(G * M / distance_from_center) altitudes.append(current_altitude_km) velocities.append(velocity) # in m/s current_altitude_km += 100 # Convert velocities from m/s to km/s and prepare for output velocities_km_s = [v / 1e3 for v in velocities] # Print the results formatted as a table print("Altitude (km)\tVelocity (km/s)") for i in range(len(altitudes)): print(f"{altitudes[i]:.2f}\t\t{velocities_km_s[i]:.4f}")
Copy the output into Excel or Calc to plot the curve.
Output looks like this (top):
python3 20250615Version02.py
Altitude (km) Velocity (km/s)
100.00 7.8483
200.00 7.7884
300.00 7.7298
400.00 7.6725
500.00 7.6165
And at GEO:
35400.00 3.0891
35500.00 3.0854
35600.00 3.0817
35700.00 3.0780
So if your business case is to create a depot in GEO, then your horizontal velocity increment would be about 3 km/s.
A depot in GEO might be most attractive for out-system shipments.
I'd appreciate someone running the numbers.
I'm assuming a shipper would appreciate the location so far above the gravity well. In addition, the location would be easy to find when navigating from out system.
(th)
Offline
Like button can go here
This program was written by Google's Gemini. It produced a set of points that can be fed into Excel or Calc to show the Earth and an ellipse created by a point launcher such as SpinLaunch or a Gas Gun.
import math import csv def generate_orbital_data(): """ Calculates orbital parameters for a launch to a specified apogee and generates points for Earth and the elliptical trajectory. The points are saved to a CSV file for plotting. """ # --- Physical Constants --- # Gravitational constant (N m^2 / kg^2) G = 6.67430e-11 # Mass of Earth (kg) M_e = 5.972e24 # Mean radius of Earth (m) R_e = 6371e3 # Target apogee altitude above surface (m) h_apogee = 300e3 # Standard gravitational parameter (mu = G * M_e) mu = G * M_e # --- Launch and Orbital Parameters --- # Radius at launch (from Earth's center) r_launch = R_e # Radius at apogee (from Earth's center) r_apogee = R_e + h_apogee # Launch elevation angle (45 degrees converted to radians) launch_angle_rad = math.radians(45) # --- Calculate Absolute Launch Velocity (v_launch_absolute) --- # This formula is derived from the conservation of energy and specific angular momentum # between the launch point and the apogee of the elliptical trajectory. # It solves for the initial velocity (v_launch_absolute) required at a given # radius (r_launch) and flight path angle (launch_angle_rad) to achieve a # specific apogee radius (r_apogee). numerator = 2 * h_apogee * mu * r_apogee denominator_term_1 = r_launch * (r_apogee**2 - r_launch**2 * math.cos(launch_angle_rad)**2) # Check if denominator_term_1 is too close to zero to avoid division by zero or very large numbers if abs(denominator_term_1) < 1e-10: # A small threshold raise ValueError("Denominator for v_launch_squared is too close to zero. Check input parameters.") v_launch_squared = numerator / denominator_term_1 # Ensure v_launch_squared is not negative, which would indicate an impossible trajectory if v_launch_squared < 0: raise ValueError("Calculated v_launch_squared is negative. The specified apogee may be unreachable with this launch angle.") v_launch_absolute = math.sqrt(v_launch_squared) # --- Calculate Velocity at Apogee of Elliptical Path (v_apogee_ellipse) --- # Using conservation of specific angular momentum (h = r * v * cos(phi)). # At apogee, the flight path angle (phi) is 0 degrees, so cos(0) = 1. v_apogee_ellipse = (r_launch * v_launch_absolute * math.cos(launch_angle_rad)) / r_apogee # --- Calculate Circular Orbit Velocity at Apogee Altitude (v_circular_at_apogee) --- # For a circular orbit, velocity is simply sqrt(mu / r). v_circular_at_apogee = math.sqrt(mu / r_apogee) # --- Calculate Delta-V (dV) for Circularization at Apogee --- # The change in velocity needed is the difference between the desired circular # velocity and the current elliptical velocity at apogee. This assumes a prograde burn. delta_v_circularization = v_circular_at_apogee - v_apogee_ellipse # --- Calculate Earth's Rotational Velocity at Equator --- # Using the sidereal day for Earth's rotation period. T_sidereal_day = 23 * 3600 + 56 * 60 + 4 # 23 hours, 56 minutes, 4 seconds v_earth_rotation = (2 * math.pi * R_e) / T_sidereal_day # --- Calculate Velocity Launcher Needs to Provide (v_launcher_relative_to_ground) --- # Since the launch is eastward from the Equator, the Earth's rotation provides a boost. # The velocity imparted by the launcher is the absolute required velocity minus # the velocity already contributed by Earth's rotation. v_launcher_relative_to_ground = v_launch_absolute - v_earth_rotation # --- Calculate Orbital Elements for the Ellipse --- # Semi-major axis (a) from the vis-viva equation: v^2 = mu * (2/r - 1/a) # Rearranging for 'a': a = r / (2 - r*v^2/mu) denominator_semi_major_axis = (2 - (r_launch * v_launch_absolute**2) / mu) if abs(denominator_semi_major_axis) < 1e-10: raise ValueError("Denominator for semi-major axis is too close to zero.") semi_major_axis = r_launch / denominator_semi_major_axis # Specific angular momentum (h) - calculated at launch specific_angular_momentum = r_launch * v_launch_absolute * math.cos(launch_angle_rad) # Eccentricity (e) # e = sqrt(1 - h^2 / (mu * a)) eccentricity_squared_term = 1 - (specific_angular_momentum**2 / (mu * semi_major_axis)) if eccentricity_squared_term < 0: # This could happen due to floating point inaccuracies, or if it's a hyperbolic trajectory. # Cap it at zero to avoid math.sqrt(negative number) eccentricity = 0.0 print("Warning: Eccentricity squared term was negative, likely due to rounding or hyperbolic trajectory. Set eccentricity to 0.") else: eccentricity = math.sqrt(eccentricity_squared_term) # True anomaly at launch (nu_0) # The true anomaly (nu) is the angle measured from perigee to the current position. # Using the polar equation of an ellipse: r = a(1-e^2) / (1 + e*cos(nu)) # Solving for cos(nu_0): cos(nu_0) = (a(1-e^2)/r_launch - 1) / e numerator_cos_nu_0 = ((semi_major_axis * (1 - eccentricity**2)) / r_launch - 1) # Check for eccentricity being zero to avoid division by zero if eccentricity == 0: nu_0 = 0.0 # For circular orbit, true anomaly doesn't matter much or is undefined print("Warning: Eccentricity is zero. True anomaly at launch set to 0.") else: cos_nu_0 = numerator_cos_nu_0 / eccentricity # Ensure cos_nu_0 is within the valid range [-1, 1] for arccos if cos_nu_0 > 1.0: cos_nu_0 = 1.0 elif cos_nu_0 < -1.0: cos_nu_0 = -1.0 # The launch point is moving towards apogee (positive flight path angle), # so nu_0 should be between 0 (perigee) and pi (apogee). nu_0 = math.acos(cos_nu_0) # --- Generate Points for Earth --- earth_points = [] num_earth_points = 360 # Points per degree for a smooth circle for i in range(num_earth_points): angle = math.radians(i) # Convert to kilometers for plotting in spreadsheet software x = R_e * math.cos(angle) / 1000 y = R_e * math.sin(angle) / 1000 earth_points.append(["Earth", x, y]) # --- Generate Points for Ellipse from Launch to Apogee --- ellipse_points = [] num_ellipse_points = 100 # Number of points for the elliptical trajectory # Iterate true anomaly from the launch point's true anomaly (nu_0) to apogee (pi). # Apogee is at true anomaly pi relative to perigee. for i in range(num_ellipse_points + 1): # Calculate current true anomaly (nu_current) along the path nu_current = nu_0 + (math.pi - nu_0) * (i / num_ellipse_points) # Calculate the radial distance at the current true anomaly # Check if denominator is zero for highly eccentric orbits or specific angles denominator_r_current = (1 + eccentricity * math.cos(nu_current)) if abs(denominator_r_current) < 1e-10: # Skip this point or handle as appropriate for a hyperbolic escape # For this problem, it should generally not be an issue for elliptic paths continue r_current = (semi_major_axis * (1 - eccentricity**2)) / denominator_r_current # Calculate X, Y coordinates for plotting. # We align the plot such that the launch point is at (R_e, 0) on the x-axis. # This means we rotate the entire ellipse by -nu_0 so the launch point is at angle 0. plot_angle = nu_current - nu_0 x = r_current * math.cos(plot_angle) / 1000 # Convert to km y = r_current * math.sin(plot_angle) / 1000 # Convert to km ellipse_points.append(["Ellipse", x, y]) # --- Write Points to CSV File --- output_filename = "orbital_points.csv" try: with open(output_filename, 'w', newline='') as csvfile: csv_writer = csv.writer(csvfile) csv_writer.writerow(["Type", "X_km", "Y_km"]) # Write header row csv_writer.writerows(earth_points) # Write Earth points csv_writer.writerows(ellipse_points) # Write Ellipse points print(f"Orbital data saved successfully to '{output_filename}'") except IOError as e: print(f"Error saving CSV file: {e}") # --- Return Calculated Results --- results = { "v_launch_absolute_mps": v_launch_absolute, "v_apogee_ellipse_mps": v_apogee_ellipse, "v_circular_at_apogee_mps": v_circular_at_apogee, "delta_v_circularization_mps": delta_v_circularization, "v_earth_rotation_mps": v_earth_rotation, "v_launcher_relative_to_ground_mps": v_launcher_relative_to_ground, "semi_major_axis_m": semi_major_axis, "eccentricity": eccentricity, "nu_0_degrees": math.degrees(nu_0), "output_filename": output_filename } return results # To run the function and get the results: results = generate_orbital_data() # You can then print the results: print("\n--- Orbital Calculation Results ---") for key, value in results.items(): if isinstance(value, (int, float)): print(f"{key.replace('_', ' ').title()}: {value:.2f}") else: print(f"{key.replace('_', ' ').title()}: {value}")
Here is the output from a test run:
python3 20250619PointsEarthLaunchV03.py
Orbital data saved successfully to 'orbital_points.csv'--- Orbital Calculation Results ---
V Launch Absolute Mps: 3216.29
V Apogee Ellipse Mps: 2171.99
V Circular At Apogee Mps: 7729.78
Delta V Circularization Mps: 5557.79 << Circularization dV is more than I had expected
V Earth Rotation Mps: 464.58 << Flight begins at the Equator. Elevation is 45 degrees
V Launcher Relative To Ground Mps: 2751.71
Semi Major Axis M: 3472589.37
Eccentricity: 0.92
Nu 0 Degrees: 174.85
Output Filename: orbital_points.csv
and here is a bit of the csv:
Type,X_km,Y_km
Earth,6371.0,0.0
Earth,6370.029665841369,111.18928141193325
Earth,6367.1189589386595,222.34469349163365
Earth,6362.26876592139,333.4323772237952
Earth,6355.48056420534,444.41849422382234
Earth,6346.756421542511,555.2692370453302
Earth,6336.098995391269,665.9508394782202
Earth,6323.511532106862,776.4295868341947
Earth,6308.997865952545,886.6718262165768
and ...
Ellipse,6646.161695217225,544.8162684888754
Ellipse,6646.201504446706,550.8337896685388
Ellipse,6646.173332272528,556.8465832168943
Ellipse,6646.077165101551,562.8544733124439
Ellipse,6645.912993236166,568.8572843678355
Ellipse,6645.6808108764335,574.8548410467276
Ellipse,6645.380616121692,580.8469682806107
Ellipse,6645.012410971555,586.833491285592
Ellipse,6644.57620132636,592.8142355791575
Ellipse,6644.071996987098,598.7890269968686
(th)
Offline
Like button can go here
This version of the Python program goes a step further, by showing the ellipse of the SpinLaunch or Gas Gun launcher, if circularization does not occur for some reason. Anyone should be able to run the program on Linux, Windows or Apple. Chrome OS would require the Linux add-in.
# Version 4 with extended ellipse import math import csv def generate_orbital_data(): """ Calculates orbital parameters for a launch to a specified apogee and generates points for Earth, the initial elliptical trajectory, and the full unburned elliptical orbit. The points are saved to a CSV file for plotting. """ # --- Physical Constants --- # Gravitational constant (N m^2 / kg^2) G = 6.67430e-11 # Mass of Earth (kg) M_e = 5.972e24 # Mean radius of Earth (m) R_e = 6371e3 # Target apogee altitude above surface (m) h_apogee = 300e3 # Standard gravitational parameter (mu = G * M_e) mu = G * M_e # --- Launch and Orbital Parameters --- # Radius at launch (from Earth's center) r_launch = R_e # Radius at apogee (from Earth's center) r_apogee = R_e + h_apogee # Launch elevation angle (45 degrees converted to radians) launch_angle_rad = math.radians(45) # --- Calculate Absolute Launch Velocity (v_launch_absolute) --- # This formula is derived from the conservation of energy and specific angular momentum # between the launch point and the apogee of the elliptical trajectory. # It solves for the initial velocity (v_launch_absolute) required at a given # radius (r_launch) and flight path angle (launch_angle_rad) to achieve a # specific apogee radius (r_apogee). numerator = 2 * h_apogee * mu * r_apogee denominator_term_1 = r_launch * (r_apogee**2 - r_launch**2 * math.cos(launch_angle_rad)**2) # Check if denominator_term_1 is too close to zero to avoid division by zero or very large numbers if abs(denominator_term_1) < 1e-10: # A small threshold raise ValueError("Denominator for v_launch_squared is too close to zero. Check input parameters.") v_launch_squared = numerator / denominator_term_1 # Ensure v_launch_squared is not negative, which would indicate an impossible trajectory if v_launch_squared < 0: raise ValueError("Calculated v_launch_squared is negative. The specified apogee may be unreachable with this launch angle.") v_launch_absolute = math.sqrt(v_launch_squared) # --- Calculate Velocity at Apogee of Elliptical Path (v_apogee_ellipse) --- # Using conservation of specific angular momentum (h = r * v * cos(phi)). # At apogee, the flight path angle (phi) is 0 degrees, so cos(0) = 1. v_apogee_ellipse = (r_launch * v_launch_absolute * math.cos(launch_angle_rad)) / r_apogee # --- Calculate Circular Orbit Velocity at Apogee Altitude (v_circular_at_apogee) --- # For a circular orbit, velocity is simply sqrt(mu / r). v_circular_at_apogee = math.sqrt(mu / r_apogee) # --- Calculate Delta-V (dV) for Circularization at Apogee --- # The change in velocity needed is the difference between the desired circular # velocity and the current elliptical velocity at apogee. This assumes a prograde burn. delta_v_circularization = v_circular_at_apogee - v_apogee_ellipse # --- Calculate Earth's Rotational Velocity at Equator --- # Using the sidereal day for Earth's rotation period. T_sidereal_day = 23 * 3600 + 56 * 60 + 4 # 23 hours, 56 minutes, 4 seconds v_earth_rotation = (2 * math.pi * R_e) / T_sidereal_day # --- Calculate Velocity Launcher Needs to Provide (v_launcher_relative_to_ground) --- # Since the launch is eastward from the Equator, the Earth's rotation provides a boost. # The velocity imparted by the launcher is the absolute required velocity minus # the velocity already contributed by Earth's rotation. v_launcher_relative_to_ground = v_launch_absolute - v_earth_rotation # --- Calculate Orbital Elements for the Ellipse --- # Semi-major axis (a) from the vis-viva equation: v^2 = mu * (2/r - 1/a) # Rearranging for 'a': a = r / (2 - r*v^2/mu) denominator_semi_major_axis = (2 - (r_launch * v_launch_absolute**2) / mu) if abs(denominator_semi_major_axis) < 1e-10: raise ValueError("Denominator for semi-major axis is too close to zero.") semi_major_axis = r_launch / denominator_semi_major_axis # Specific angular momentum (h) - calculated at launch specific_angular_momentum = r_launch * v_launch_absolute * math.cos(launch_angle_rad) # Eccentricity (e) # e = sqrt(1 - h^2 / (mu * a)) eccentricity_squared_term = 1 - (specific_angular_momentum**2 / (mu * semi_major_axis)) if eccentricity_squared_term < 0: # This could happen due to floating point inaccuracies, or if it's a hyperbolic trajectory. # Cap it at zero to avoid math.sqrt(negative number) eccentricity = 0.0 print("Warning: Eccentricity squared term was negative, likely due to rounding or hyperbolic trajectory. Set eccentricity to 0.") else: eccentricity = math.sqrt(eccentricity_squared_term) # True anomaly at launch (nu_0) # The true anomaly (nu) is the angle measured from perigee to the current position. # Using the polar equation of an ellipse: r = a(1-e^2) / (1 + e*cos(nu)) # Solving for cos(nu_0): cos(nu_0) = (a(1-e^2)/r_launch - 1) / e numerator_cos_nu_0 = ((semi_major_axis * (1 - eccentricity**2)) / r_launch - 1) # Check for eccentricity being zero to avoid division by zero if eccentricity == 0: nu_0 = 0.0 # For circular orbit, true anomaly doesn't matter much or is undefined print("Warning: Eccentricity is zero. True anomaly at launch set to 0.") else: cos_nu_0 = numerator_cos_nu_0 / eccentricity # Ensure cos_nu_0 is within the valid range [-1, 1] for arccos if cos_nu_0 > 1.0: cos_nu_0 = 1.0 elif cos_nu_0 < -1.0: cos_nu_0 = -1.0 # The launch point is moving towards apogee (positive flight path angle), # so nu_0 should be between 0 (perigee) and pi (apogee). nu_0 = math.acos(cos_nu_0) # --- Generate Points for Earth --- earth_points = [] num_earth_points = 360 # Points per degree for a smooth circle for i in range(num_earth_points): angle = math.radians(i) # Convert to kilometers for plotting in spreadsheet software x = R_e * math.cos(angle) / 1000 y = R_e * math.sin(angle) / 1000 earth_points.append(["Earth", x, y]) # --- Generate Points for Ellipse from Launch to Apogee --- # This segment only goes from the launch point to the apogee. ellipse_launch_to_apogee_points = [] num_segment_points = 100 # Number of points for the elliptical trajectory segment # Iterate true anomaly from the launch point's true anomaly (nu_0) to apogee (pi). # Apogee is at true anomaly pi relative to perigee. for i in range(num_segment_points + 1): # Calculate current true anomaly (nu_current) along the path nu_current = nu_0 + (math.pi - nu_0) * (i / num_segment_points) # Calculate the radial distance at the current true anomaly # Check if denominator is zero for highly eccentric orbits or specific angles denominator_r_current = (1 + eccentricity * math.cos(nu_current)) if abs(denominator_r_current) < 1e-10: # Skip this point, as it would represent an infinitely far point in a parabolic/hyperbolic case continue r_current = (semi_major_axis * (1 - eccentricity**2)) / denominator_r_current # Calculate X, Y coordinates for plotting. # We align the plot such that the launch point is at (R_e, 0) on the x-axis. # This means we rotate the entire ellipse by -nu_0 so the launch point is at angle 0. plot_angle = nu_current - nu_0 x = r_current * math.cos(plot_angle) / 1000 # Convert to km y = r_current * math.sin(plot_angle) / 1000 # Convert to km ellipse_launch_to_apogee_points.append(["Ellipse_Launch_To_Apogee", x, y]) # --- Generate Points for Full Unburned Ellipse --- # This generates points for the entire 360-degree elliptical orbit. full_unburned_ellipse_points = [] num_full_ellipse_points = 200 # More points for a smoother full ellipse for i in range(num_full_ellipse_points): # True anomaly from 0 to 2*pi nu_full_ellipse = (2 * math.pi * i) / num_full_ellipse_points denominator_r_full = (1 + eccentricity * math.cos(nu_full_ellipse)) if abs(denominator_r_full) < 1e-10: # Skip points if denominator is zero (e.g., at certain points in parabolic/hyperbolic orbits) continue r_full_ellipse = (semi_major_axis * (1 - eccentricity**2)) / denominator_r_full # Apply the same rotation as before to align the launch point with the positive X-axis plot_angle_full_ellipse = nu_full_ellipse - nu_0 x_full = r_full_ellipse * math.cos(plot_angle_full_ellipse) / 1000 # Convert to km y_full = r_full_ellipse * math.sin(plot_angle_full_ellipse) / 1000 # Convert to km full_unburned_ellipse_points.append(["Full_Unburned_Ellipse", x_full, y_full]) # --- Write Points to CSV File --- output_filename = "orbital_points.csv" try: with open(output_filename, 'w', newline='') as csvfile: csv_writer = csv.writer(csvfile) csv_writer.writerow(["Type", "X_km", "Y_km"]) # Write header row csv_writer.writerows(earth_points) # Write Earth points csv_writer.writerows(ellipse_launch_to_apogee_points) # Write initial segment csv_writer.writerows(full_unburned_ellipse_points) # Write full ellipse points print(f"Orbital data saved successfully to '{output_filename}'") except IOError as e: print(f"Error saving CSV file: {e}") # --- Return Calculated Results --- results = { "v_launch_absolute_mps": v_launch_absolute, "v_apogee_ellipse_mps": v_apogee_ellipse, "v_circular_at_apogee_mps": v_circular_at_apogee, "delta_v_circularization_mps": delta_v_circularization, "v_earth_rotation_mps": v_earth_rotation, "v_launcher_relative_to_ground_mps": v_launcher_relative_to_ground, "semi_major_axis_m": semi_major_axis, "eccentricity": eccentricity, "nu_0_degrees": math.degrees(nu_0), "output_filename": output_filename } return results # To run the function and get the results: results = generate_orbital_data() # You can then print the results: print("\n--- Orbital Calculation Results ---") for key, value in results.items(): if isinstance(value, (int, float)): print(f"{key.replace('_', ' ').title()}: {value:.2f}") else: print(f"{key.replace('_', ' ').title()}: {value}")
Top of CSV file;
Type,X_km,Y_km
Earth,6371.0,0.0
Earth,6370.029665841369,111.18928141193325
Earth,6367.1189589386595,222.34469349163365
Earth,6362.26876592139,333.4323772237952
Earth,6355.48056420534,444.41849422382234
Earth,6346.756421542511,555.2692370453302
Earth,6336.098995391269,665.9508394782202
Earth,6323.511532106862,776.4295868341947
Earth,6308.997865952545,886.6718262165768
Note for novices (like I am regarding Excel)... the data is separated into three ranges:
Range 1 is for Earth so it has to be specified in both X and Y
Range 2 is for the stub of the ellipse if circularization occurs
Range 3 is for the full ellipse if circularization does not occur
Naturally in the Real Universe, the payload would impact the Earth, but this plot is for educational purposes, showing what the ellipse would look like if the Earth were a point mass like a black hole.
(th)
Offline
Like button can go here
This post is about an extension of the initial Python program designed to reach LEO and to circularize there.
The extension appears to confirm my hypothesis (which I've been working on for a couple of years at least) that a launch from a point source that reaches GEO can achieve orbit at LEO at a significant savings of dV for the payload. This program provides data and the source code will be uploaded so all who might be interested can study the procedure. Output is a set of data that can be loaded into Excel or Calc to produce graphs showing the flight.
Here is output at run time:
python3 GEOtoLEOV02.py
Orbital data saved successfully to 'orbital_points_multi_burn.csv'--- Multi-Burn Orbital Calculation Results ---
H Apogee Initial Target Km: 35786.000
H Leo Target Km: 300.000
V Launch Absolute Mps: 10365.486
V Launcher Relative To Ground Mps: 9900.905
V At Geo Apogee Initial Ellipse Mps: 1107.676
Delta V Burn1 At Geo Mps: 499.650
V At Leo Perigee Transfer Ellipse Mps: 10157.403
Delta V Burn2 At Leo Mps: -2427.620
V Circular Leo Mps: 7729.783
Total Delta V Burns Mps: 2927.270
A Initial Ellipse Km: 22541.060
E Initial Ellipse: 0.870
Nu 0 Initial Ellipse Degrees: 99.346
A Transfer Ellipse Km: 24414.000
E Transfer Ellipse: 0.727
Output Filename: orbital_points_multi_burn.csv
The output is (happily) in (approximate) agreement with figures produced with an Excel tool.
more will be added ...
(th)
Offline
Like button can go here
This post will contain the entire code for a Python program to demonstrate the advantage of launching all the way to GEO in order to arrive at LEO with a modest cost in secondary propulsion compared to a direct flight to LEO.
This version writes plot points to a CSV file for import to Excel or Calc.
# Version 3 of GEOtoLEO program import math import csv def generate_orbital_data_multi_burn(): """ Calculates orbital parameters for a multi-burn scenario: 1. Initial launch to GEO apogee. 2. Burn at GEO apogee to raise perigee to LEO. 3. Burn at LEO perigee to circularize. Generates points for Earth, initial ellipse, transfer ellipse, and final circular orbit. The points are saved to a CSV file for plotting. """ # --- Physical Constants --- # Gravitational constant (N m^2 / kg^2) G = 6.67430e-11 # Mass of Earth (kg) M_e = 5.972e24 # Mean radius of Earth (m) R_e = 6371e3 # Standard gravitational parameter (mu = G * M_e) mu = G * M_e # --- Target Orbital Altitudes (above surface) --- # Target apogee altitude for initial launch (GEO altitude, m) h_apogee_initial_target = 35786e3 # Target perigee altitude for transfer ellipse / final LEO altitude (m) h_leo_target = 300e3 # --- Radii (from Earth's center) --- r_launch = R_e r_apogee_initial = R_e + h_apogee_initial_target # Radius at GEO altitude r_leo_circular = R_e + h_leo_target # Radius at 300 km LEO altitude # Launch elevation angle (45 degrees converted to radians) launch_angle_rad = math.radians(45) # --- PART 1: Initial Launch from Earth to GEO Apogee --- # Calculate Absolute Launch Velocity (v_launch_absolute) # This formula is derived from the conservation of energy and specific angular momentum # between the launch point and the apogee of the initial elliptical trajectory. # It solves for the initial velocity (v_launch_absolute) required at a given # radius (r_launch) and flight path angle (launch_angle_rad) to achieve a # specific apogee radius (r_apogee_initial). numerator = 2 * h_apogee_initial_target * mu * r_apogee_initial denominator_term_1 = r_launch * (r_apogee_initial**2 - r_launch**2 * math.cos(launch_angle_rad)**2) if abs(denominator_term_1) < 1e-10: raise ValueError("Denominator for v_launch_squared is too close to zero. Check input parameters.") v_launch_squared = numerator / denominator_term_1 if v_launch_squared < 0: raise ValueError("Calculated v_launch_squared is negative. The specified apogee may be unreachable with this launch angle.") v_launch_absolute = math.sqrt(v_launch_squared) # Calculate Velocity at Apogee of Initial Elliptical Path (v_at_geo_apogee_initial_ellipse) # Using conservation of specific angular momentum (h = r * v * cos(phi)). # At apogee, the flight path angle (phi) is 0 degrees, so cos(0) = 1. v_at_geo_apogee_initial_ellipse = (r_launch * v_launch_absolute * math.cos(launch_angle_rad)) / r_apogee_initial # Calculate Earth's Rotational Velocity at Equator # Using the sidereal day for Earth's rotation period. T_sidereal_day = 23 * 3600 + 56 * 60 + 4 # 23 hours, 56 minutes, 4 seconds v_earth_rotation = (2 * math.pi * R_e) / T_sidereal_day # Calculate Velocity Launcher Needs to Provide (v_launcher_relative_to_ground) # Since the launch is eastward from the Equator, Earth's rotation provides a boost. # The velocity imparted by the launcher is the absolute required velocity minus # the velocity already contributed by Earth's rotation. v_launcher_relative_to_ground = v_launch_absolute - v_earth_rotation # Calculate Orbital Elements for the Initial Ellipse (from Earth to GEO Apogee) # Semi-major axis (a_initial_ellipse) from the vis-viva equation: v^2 = mu * (2/r - 1/a) # Rearranging for 'a': a = r / (2 - r*v^2/mu) denominator_a_initial = (2 - (r_launch * v_launch_absolute**2) / mu) if abs(denominator_a_initial) < 1e-10: raise ValueError("Denominator for initial semi-major axis is too close to zero.") a_initial_ellipse = r_launch / denominator_a_initial # Specific angular momentum (h) - calculated at launch h_initial = r_launch * v_launch_absolute * math.cos(launch_angle_rad) # Eccentricity (e_initial_ellipse) eccentricity_squared_term_initial = 1 - (h_initial**2 / (mu * a_initial_ellipse)) if eccentricity_squared_term_initial < 0: eccentricity_initial_ellipse = 0.0 print("Warning: Initial ellipse eccentricity squared term was negative, likely due to rounding or hyperbolic trajectory. Set eccentricity to 0.") else: eccentricity_initial_ellipse = math.sqrt(eccentricity_squared_term_initial) # True anomaly at launch (nu_0_initial_ellipse) # This aligns the launch point on the positive X-axis for plotting convenience. numerator_cos_nu_0_initial = ((a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / r_launch - 1) if eccentricity_initial_ellipse == 0: nu_0_initial_ellipse = 0.0 print("Warning: Initial ellipse eccentricity is zero. True anomaly at launch set to 0.") else: cos_nu_0_initial = numerator_cos_nu_0_initial / eccentricity_initial_ellipse cos_nu_0_initial = max(-1.0, min(1.0, cos_nu_0_initial)) # Clamp to [-1, 1] nu_0_initial_ellipse = math.acos(cos_nu_0_initial) # --- PART 2: Burn 1 at GEO Apogee to Raise Perigee to LEO Altitude --- # This burn will create a new transfer ellipse whose apogee is at GEO altitude # and whose perigee is at LEO altitude (300 km). # Semi-major axis of the transfer ellipse a_transfer_ellipse = (r_apogee_initial + r_leo_circular) / 2 # Velocity required at GEO (apogee of transfer ellipse) # From vis-viva equation: v = sqrt(mu * (2/r - 1/a)) v_transfer_apogee = math.sqrt(mu * (2/r_apogee_initial - 1/a_transfer_ellipse)) # Delta-V for the first burn (at GEO apogee) # We are increasing velocity slightly (prograde burn) to raise the perigee from R_e to r_leo_circular. delta_v_burn1_at_geo = v_transfer_apogee - v_at_geo_apogee_initial_ellipse # Calculate Orbital Elements for the Transfer Ellipse (GEO Apogee to LEO Perigee) # Eccentricity of the transfer ellipse eccentricity_transfer_ellipse = (r_apogee_initial - r_leo_circular) / (r_apogee_initial + r_leo_circular) # Velocity at LEO perigee of the transfer ellipse v_transfer_perigee = math.sqrt(mu * (2/r_leo_circular - 1/a_transfer_ellipse)) # --- PART 3: Burn 2 at LEO Perigee to Circularize Orbit --- # This burn will occur when the payload reaches the 300 km LEO altitude, # changing its velocity to match the circular orbit velocity at that altitude. # Velocity required for a circular orbit at 300 km LEO v_circular_leo = math.sqrt(mu / r_leo_circular) # Delta-V for the second burn (at LEO perigee) # This will be a retrograde burn, decreasing velocity to circularize. delta_v_burn2_at_leo = v_circular_leo - v_transfer_perigee # --- GENERATE POINTS FOR PLOTTING --- # List to store all generated points all_points = [] # --- Generate Points for Earth --- num_earth_points = 361 # 360 + 1 to close the circle for i in range(num_earth_points): angle = (2 * math.pi * i) / (num_earth_points - 1) # Ensure last point is 2pi x = R_e * math.cos(angle) / 1000 # Convert to kilometers y = R_e * math.sin(angle) / 1000 # Convert to kilometers all_points.append(["Earth", x, y]) # --- Generate Points for Initial Launch Ellipse (Earth to GEO Apogee) --- num_segment_points = 100 for i in range(num_segment_points + 1): # True anomaly from launch point (nu_0_initial_ellipse) to apogee (pi) nu_current = nu_0_initial_ellipse + (math.pi - nu_0_initial_ellipse) * (i / num_segment_points) # Calculate radial distance denominator_r_current = (1 + eccentricity_initial_ellipse * math.cos(nu_current)) if abs(denominator_r_current) < 1e-10: continue r_current = (a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / denominator_r_current # Calculate X, Y coordinates for plotting, aligning launch point with +X axis plot_angle = nu_current - nu_0_initial_ellipse x = r_current * math.cos(plot_angle) / 1000 y = r_current * math.sin(plot_angle) / 1000 all_points.append(["Initial_Launch_Ellipse", x, y]) # --- Generate Points for Transfer Ellipse (GEO Apogee to LEO Perigee) --- num_transfer_points = 100 # Calculate the plot angle of the initial ellipse's apogee (where transfer starts) apogee_angle_in_plot = math.pi - nu_0_initial_ellipse # The target plot angle for the LEO perigee (opposite launch, on X-axis) target_perigee_angle_in_plot = math.pi # 180 degrees for i in range(num_transfer_points + 1): # Linearly interpolate the plotting angle between the two end points # This is a visual adjustment to ensure the path appears as desired, # flowing smoothly from the initial apogee to the LEO perigee on the far side. interp_factor = i / num_transfer_points plot_angle_transfer = apogee_angle_in_plot * (1 - interp_factor) + target_perigee_angle_in_plot * interp_factor # The true anomaly for the transfer ellipse, sweeping from apogee (pi) to perigee (0) nu_transfer = math.pi - (math.pi * i / num_transfer_points) # Calculate radial distance based on the transfer ellipse's elements and its true anomaly denominator_r_transfer = (1 + eccentricity_transfer_ellipse * math.cos(nu_transfer)) if abs(denominator_r_transfer) < 1e-10: continue r_transfer = (a_transfer_ellipse * (1 - eccentricity_transfer_ellipse**2)) / denominator_r_transfer # Calculate X, Y coordinates using the interpolated plot angle and the radial distance x = r_transfer * math.cos(plot_angle_transfer) / 1000 y = r_transfer * math.sin(plot_angle_transfer) / 1000 all_points.append(["Transfer_Ellipse", x, y]) # --- Generate Points for Circular LEO Orbit (300 km) --- num_circular_points = 201 # +1 to close the circle # The LEO circular orbit will be centered at (0,0) with radius r_leo_circular. for i in range(num_circular_points): angle = (2 * math.pi * i) / (num_circular_points - 1) # Ensure last point is 2pi x = r_leo_circular * math.cos(angle) / 1000 y = r_leo_circular * math.sin(angle) / 1000 all_points.append(["Circular_LEO_Orbit", x, y]) # --- Generate Points for Full Unburned Initial Ellipse (for comparison) --- num_full_ellipse_points = 201 # +1 to close the ellipse for i in range(num_full_ellipse_points): # True anomaly from 0 to 2*pi nu_full_ellipse = (2 * math.pi * i) / (num_full_ellipse_points - 1) # Ensure last point is 2pi denominator_r_full = (1 + eccentricity_initial_ellipse * math.cos(nu_full_ellipse)) if abs(denominator_r_full) < 1e-10: continue r_full_ellipse = (a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / denominator_r_full # Apply the same rotation as before to align the launch point with the positive X-axis plot_angle_full_ellipse = nu_full_ellipse - nu_0_initial_ellipse x_full = r_full_ellipse * math.cos(plot_angle_full_ellipse) / 1000 y_full = r_full_ellipse * math.sin(plot_angle_full_ellipse) / 1000 all_points.append(["Full_Unburned_Ellipse", x_full, y_full]) # --- Write Points to CSV File --- output_filename = "orbital_points_multi_burn.csv" # New filename to distinguish try: with open(output_filename, 'w', newline='') as csvfile: csv_writer = csv.writer(csvfile) csv_writer.writerow(["Type", "X_km", "Y_km"]) # Write header row csv_writer.writerows(all_points) print(f"Orbital data saved successfully to '{output_filename}'") except IOError as e: print(f"Error saving CSV file: {e}") # --- Return Calculated Results --- results = { "h_apogee_initial_target_km": h_apogee_initial_target / 1000, "h_leo_target_km": h_leo_target / 1000, "v_launch_absolute_mps": v_launch_absolute, "v_launcher_relative_to_ground_mps": v_launcher_relative_to_ground, "v_at_geo_apogee_initial_ellipse_mps": v_at_geo_apogee_initial_ellipse, "delta_v_burn1_at_geo_mps": delta_v_burn1_at_geo, # This will be positive for prograde burn "v_at_leo_perigee_transfer_ellipse_mps": v_transfer_perigee, "delta_v_burn2_at_leo_mps": delta_v_burn2_at_leo, # This will be negative for retrograde burn "v_circular_leo_mps": v_circular_leo, "total_delta_v_burns_mps": abs(delta_v_burn1_at_geo) + abs(delta_v_burn2_at_leo), "a_initial_ellipse_km": a_initial_ellipse / 1000, "e_initial_ellipse": eccentricity_initial_ellipse, "nu_0_initial_ellipse_degrees": math.degrees(nu_0_initial_ellipse), "a_transfer_ellipse_km": a_transfer_ellipse / 1000, "e_transfer_ellipse": eccentricity_transfer_ellipse, "output_filename": output_filename } return results # To run the function and get the results: results = generate_orbital_data_multi_burn() # You can then print the results: print("\n--- Multi-Burn Orbital Calculation Results ---") for key, value in results.items(): if isinstance(value, (int, float)): print(f"{key.replace('_', ' ').title()}: {value:.3f}") # More precision for velocities else: print(f"{key.replace('_', ' ').title()}: {value}") # 266 lines in program V03
(th)
Offline
Like button can go here
This post contains a version of the GEO to LEO program that outputs to the monitor. Output shows:
1) The Earth
2) Launch trajectory to GEO
3) Trajectory on Transfer orbit to LEO
4) Orbit at LEO
5) Default ellipse if dV is not applied at GEO
cat TkinterVersion04.py # Version 4 has Close in Header import math import csv import matplotlib.pyplot as plt from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg import tkinter as tk import sys # Import the sys module def generate_orbital_data_and_plot(): """ Calculates orbital parameters for a multi-burn scenario: 1. Initial launch to GEO apogee. 2. Burn at GEO apogee to raise perigee to LEO. 3. Burn at LEO perigee to circularize. Generates points for Earth, initial ellipse, transfer ellipse, and final circular orbit. Displays the plot in a Tkinter window using Matplotlib. """ # --- Physical Constants --- # Gravitational constant (N m^2 / kg^2) G = 6.67430e-11 # Mass of Earth (kg) M_e = 5.972e24 # Mean radius of Earth (m) R_e = 6371e3 # Standard gravitational parameter (mu = G * M_e) mu = G * M_e # --- Target Orbital Altitudes (above surface) --- # Target apogee altitude for initial launch (GEO altitude, m) h_apogee_initial_target = 35786e3 # Target perigee altitude for transfer ellipse / final LEO altitude (m) h_leo_target = 300e3 # --- Radii (from Earth's center) --- r_launch = R_e r_apogee_initial = R_e + h_apogee_initial_target # Radius at GEO altitude r_leo_circular = R_e + h_leo_target # Radius at 300 km LEO altitude # Launch elevation angle (45 degrees converted to radians) launch_angle_rad = math.radians(45) # --- PART 1: Initial Launch from Earth to GEO Apogee --- # Calculate Absolute Launch Velocity (v_launch_absolute) # This formula is derived from the conservation of energy and specific angular momentum # between the launch point and the apogee of the initial elliptical trajectory. # It solves for the initial velocity (v_launch_absolute) required at a given # radius (r_launch) and flight path angle (launch_angle_rad) to achieve a # specific apogee radius (r_apogee_initial). numerator = 2 * h_apogee_initial_target * mu * r_apogee_initial denominator_term_1 = r_launch * (r_apogee_initial**2 - r_launch**2 * math.cos(launch_angle_rad)**2) if abs(denominator_term_1) < 1e-10: raise ValueError("Denominator for v_launch_squared is too close to zero. Check input parameters.") v_launch_squared = numerator / denominator_term_1 if v_launch_squared < 0: raise ValueError("Calculated v_launch_squared is negative. The specified apogee may be unreachable with this launch angle.") v_launch_absolute = math.sqrt(v_launch_squared) # Calculate Velocity at Apogee of Initial Elliptical Path (v_at_geo_apogee_initial_ellipse) # Using conservation of specific angular momentum (h = r * v * cos(phi)). # At apogee, the flight path angle (phi) is 0 degrees, so cos(0) = 1. v_at_geo_apogee_initial_ellipse = (r_launch * v_launch_absolute * math.cos(launch_angle_rad)) / r_apogee_initial # Calculate Earth's Rotational Velocity at Equator # Using the sidereal day for Earth's rotation period. T_sidereal_day = 23 * 3600 + 56 * 60 + 4 # 23 hours, 56 minutes, 4 seconds v_earth_rotation = (2 * math.pi * R_e) / T_sidereal_day # Calculate Velocity Launcher Needs to Provide (v_launcher_relative_to_ground) # Since the launch is eastward from the Equator, Earth's rotation provides a boost. # The velocity imparted by the launcher is the absolute required velocity minus # the velocity already contributed by Earth's rotation. v_launcher_relative_to_ground = v_launch_absolute - v_earth_rotation # Calculate Orbital Elements for the Initial Ellipse (from Earth to GEO Apogee) # Semi-major axis (a_initial_ellipse) from the vis-viva equation: v^2 = mu * (2/r - 1/a) # Rearranging for 'a': a = r / (2 - r*v^2/mu) denominator_a_initial = (2 - (r_launch * v_launch_absolute**2) / mu) if abs(denominator_a_initial) < 1e-10: raise ValueError("Denominator for initial semi-major axis is too close to zero.") a_initial_ellipse = r_launch / denominator_a_initial # Specific angular momentum (h) - calculated at launch h_initial = r_launch * v_launch_absolute * math.cos(launch_angle_rad) # Eccentricity (e_initial_ellipse) eccentricity_squared_term_initial = 1 - (h_initial**2 / (mu * a_initial_ellipse)) if eccentricity_squared_term_initial < 0: eccentricity_initial_ellipse = 0.0 print("Warning: Initial ellipse eccentricity squared term was negative, likely due to rounding or hyperbolic trajectory. Set eccentricity to 0.") else: eccentricity_initial_ellipse = math.sqrt(eccentricity_squared_term_initial) # True anomaly at launch (nu_0_initial_ellipse) # This aligns the launch point on the positive X-axis for plotting convenience. numerator_cos_nu_0_initial = ((a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / r_launch - 1) if eccentricity_initial_ellipse == 0: nu_0_initial_ellipse = 0.0 print("Warning: Initial ellipse eccentricity is zero. True anomaly at launch set to 0.") else: cos_nu_0_initial = numerator_cos_nu_0_initial / eccentricity_initial_ellipse cos_nu_0_initial = max(-1.0, min(1.0, cos_nu_0_initial)) # Clamp to [-1, 1] nu_0_initial_ellipse = math.acos(cos_nu_0_initial) # --- PART 2: Burn 1 at GEO Apogee to Raise Perigee to LEO Altitude --- # This burn will create a new transfer ellipse whose apogee is at GEO altitude # and whose perigee is at LEO altitude (300 km). # Semi-major axis of the transfer ellipse a_transfer_ellipse = (r_apogee_initial + r_leo_circular) / 2 # Velocity required at GEO (apogee of transfer ellipse) # From vis-viva equation: v = sqrt(mu * (2/r - 1/a)) v_transfer_apogee = math.sqrt(mu * (2/r_apogee_initial - 1/a_transfer_ellipse)) # Delta-V for the first burn (at GEO apogee) # We are increasing velocity slightly (prograde burn) to raise the perigee from R_e to r_leo_circular. delta_v_burn1_at_geo = v_transfer_apogee - v_at_geo_apogee_initial_ellipse # Calculate Orbital Elements for the Transfer Ellipse (GEO Apogee to LEO Perigee) # Eccentricity of the transfer ellipse eccentricity_transfer_ellipse = (r_apogee_initial - r_leo_circular) / (r_apogee_initial + r_leo_circular) # Velocity at LEO perigee of the transfer ellipse v_transfer_perigee = math.sqrt(mu * (2/r_leo_circular - 1/a_transfer_ellipse)) # --- PART 3: Burn 2 at LEO Perigee to Circularize Orbit --- # This burn will occur when the payload reaches the 300 km LEO altitude, # changing its velocity to match the circular orbit velocity at that altitude. # Velocity required for a circular orbit at 300 km LEO v_circular_leo = math.sqrt(mu / r_leo_circular) # Delta-V for the second burn (at LEO perigee) # This will be a retrograde burn, decreasing velocity to circularize. delta_v_burn2_at_leo = v_circular_leo - v_transfer_perigee # --- GENERATE POINTS FOR PLOTTING --- # Initialize lists to store X and Y coordinates for each series earth_x, earth_y = [], [] initial_launch_ellipse_x, initial_launch_ellipse_y = [], [] transfer_ellipse_x, transfer_ellipse_y = [], [] circular_leo_orbit_x, circular_leo_orbit_y = [], [] full_unburned_ellipse_x, full_unburned_ellipse_y = [], [] # --- Generate Points for Earth --- num_earth_points = 361 # 360 + 1 to close the circle for i in range(num_earth_points): angle = (2 * math.pi * i) / (num_earth_points - 1) # Ensure last point is 2pi x = R_e * math.cos(angle) / 1000 # Convert to kilometers y = R_e * math.sin(angle) / 1000 # Convert to kilometers earth_x.append(x) earth_y.append(y) # --- Generate Points for Initial Launch Ellipse (Earth to GEO Apogee) --- num_segment_points = 100 for i in range(num_segment_points + 1): # True anomaly from launch point (nu_0_initial_ellipse) to apogee (pi) nu_current = nu_0_initial_ellipse + (math.pi - nu_0_initial_ellipse) * (i / num_segment_points) # Calculate radial distance denominator_r_current = (1 + eccentricity_initial_ellipse * math.cos(nu_current)) if abs(denominator_r_current) < 1e-10: continue r_current = (a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / denominator_r_current # Calculate X, Y coordinates for plotting, aligning launch point with +X axis plot_angle = nu_current - nu_0_initial_ellipse x = r_current * math.cos(plot_angle) / 1000 y = r_current * math.sin(plot_angle) / 1000 initial_launch_ellipse_x.append(x) initial_launch_ellipse_y.append(y) # --- Generate Points for Transfer Ellipse (GEO Apogee to LEO Perigee) --- num_transfer_points = 100 # Get the X, Y coordinates of the GEO apogee from the initial ellipse for plotting start geo_apogee_x_plot_start = r_apogee_initial * math.cos(math.pi - nu_0_initial_ellipse) / 1000 geo_apogee_y_plot_start = r_apogee_initial * math.sin(math.pi - nu_0_initial_ellipse) / 1000 # The target LEO perigee visual point for the end of the transfer ellipse plotting segment # We want this to be at (-r_leo_circular, 0.0) for visual alignment leo_perigee_x_plot_end = -r_leo_circular / 1000 leo_perigee_y_plot_end = 0.0 # Calculate the initial and final angles in the plotting frame for linear interpolation angle_start_plot_transfer = math.atan2(geo_apogee_y_plot_start, geo_apogee_x_plot_start) angle_end_plot_transfer = math.atan2(leo_perigee_y_plot_end, leo_perigee_x_plot_end) # Ensure angles go in the correct direction (e.g., if we cross the -X axis from +Y to -Y) # The initial angle will be in Q1 or Q2. The end angle is math.pi. # If the start angle is greater than pi, it means it's in Q2 and we want to sweep counter-clockwise towards pi. # If the start angle is less than pi but still positive, it's also a counter-clockwise sweep. # The original interpolation `angle_start_plot_transfer * (1 - interp_factor) + angle_end_plot_transfer * interp_factor` # works for a direct interpolation. for i in range(num_transfer_points + 1): # Linearly interpolate the true anomaly (nu_transfer) from pi (apogee) down to 0 (perigee) nu_transfer = math.pi - (math.pi * i / num_transfer_points) # Calculate radial distance for the transfer ellipse based on its true orbital elements denominator_r_transfer = (1 + eccentricity_transfer_ellipse * math.cos(nu_transfer)) if abs(denominator_r_transfer) < 1e-10: continue r_transfer = (a_transfer_ellipse * (1 - eccentricity_transfer_ellipse**2)) / denominator_r_transfer # Linearly interpolate the *plotting angle* between the start and end visual points. interp_factor = i / num_transfer_points current_plot_angle = angle_start_plot_transfer * (1 - interp_factor) + angle_end_plot_transfer * interp_factor x = r_transfer * math.cos(current_plot_angle) / 1000 y = r_transfer * math.sin(current_plot_angle) / 1000 transfer_ellipse_x.append(x) transfer_ellipse_y.append(y) # --- Generate Points for Circular LEO Orbit (300 km) --- num_circular_points = 201 # +1 to close the circle # The LEO circular orbit will be centered at (0,0) with radius r_leo_circular. for i in range(num_circular_points): angle = (2 * math.pi * i) / (num_circular_points - 1) # Ensure last point is 2pi x = r_leo_circular * math.cos(angle) / 1000 y = r_leo_circular * math.sin(angle) / 1000 circular_leo_orbit_x.append(x) circular_leo_orbit_y.append(y) # --- Generate Points for Full Unburned Initial Ellipse (for comparison) --- num_full_ellipse_points = 201 # +1 to close the ellipse for i in range(num_full_ellipse_points): # True anomaly from 0 to 2*pi nu_full_ellipse = (2 * math.pi * i) / (num_full_ellipse_points - 1) # Ensure last point is 2pi denominator_r_full = (1 + eccentricity_initial_ellipse * math.cos(nu_full_ellipse)) if abs(denominator_r_full) < 1e-10: continue r_full_ellipse = (a_initial_ellipse * (1 - eccentricity_initial_ellipse**2)) / denominator_r_full # Apply the same rotation as before to align the launch point with the positive X-axis plot_angle_full_ellipse = nu_full_ellipse - nu_0_initial_ellipse x_full = r_full_ellipse * math.cos(plot_angle_full_ellipse) / 1000 y_full = r_full_ellipse * math.sin(plot_angle_full_ellipse) / 1000 full_unburned_ellipse_x.append(x_full) full_unburned_ellipse_y.append(y_full) # --- Plotting with Matplotlib and Tkinter --- root = tk.Tk() root.title("Orbital Trajectories") # Function to cleanly close the program def close_program(): root.destroy() # Destroy the Tkinter window sys.exit() # Exit the Python script # Set the window close protocol to call our custom close_program function root.protocol("WM_DELETE_WINDOW", close_program) # Create a frame for controls (e.g., Close button) control_frame = tk.Frame(root) control_frame.pack(side=tk.TOP, fill=tk.X, padx=5, pady=5) # pack at top, fill horizontally # Add a Close Button within the control frame close_button = tk.Button(control_frame, text="Close Plot", command=close_program) close_button.pack(side=tk.LEFT) # Pack to the left within the frame # Add a title label to the control frame (optional, for aesthetics) title_label = tk.Label(control_frame, text="Orbital Trajectories Simulation", font=("Arial", 12, "bold")) title_label.pack(side=tk.LEFT, padx=10) # Pack to the left, next to the button fig, ax = plt.subplots(figsize=(10, 10)) # Create a figure and a set of subplots # Plotting each series ax.plot(earth_x, earth_y, color='blue', label='Earth Surface', linewidth=2) ax.plot(initial_launch_ellipse_x, initial_launch_ellipse_y, color='green', linestyle='-', label='Initial Launch Ellipse', linewidth=2) ax.plot(transfer_ellipse_x, transfer_ellipse_y, color='purple', linestyle='--', label='Transfer Ellipse', linewidth=2) ax.plot(circular_leo_orbit_x, circular_leo_orbit_y, color='red', linestyle='-', label='Circular LEO Orbit (300 km)', linewidth=2) ax.plot(full_unburned_ellipse_x, full_unburned_ellipse_y, color='orange', linestyle=':', label='Full Unburned Initial Ellipse', linewidth=1) # Add markers for key points for clarity # Launch point (on Earth) ax.plot([R_e/1000], [0], 'go', markersize=8, label='Launch Point') # GEO Apogee (start of transfer) ax.plot([geo_apogee_x_plot_start], [geo_apogee_y_plot_start], 'g*', markersize=10, label='GEO Apogee (Burn 1)') # LEO Perigee (end of transfer, start of circular) - This marker should be at (-R_LEO, 0) ax.plot([leo_perigee_x_plot_end], [leo_perigee_y_plot_end], 'rX', markersize=10, label='LEO Perigee (Burn 2)') ax.set_aspect('equal', adjustable='box') # Ensures circles are circular ax.set_xlabel("X (km)") ax.set_ylabel("Y (km)") ax.set_title("Orbital Trajectories Overview") # Title within the plot area ax.legend(loc='upper right') ax.grid(True) # Set plot limits dynamically to comfortably fit all orbits max_radius_km = r_apogee_initial / 1000 padding = max_radius_km * 0.1 # Add 10% padding ax.set_xlim(-(max_radius_km + padding), (max_radius_km + padding)) ax.set_ylim(-(max_radius_km + padding), (max_radius_km + padding)) canvas = FigureCanvasTkAgg(fig, master=root) canvas_widget = canvas.get_tk_widget() canvas_widget.pack(side=tk.TOP, fill=tk.BOTH, expand=1) # Plot takes up remaining space # Display calculated results in the console/terminal results = { "h_apogee_initial_target_km": h_apogee_initial_target / 1000, "h_leo_target_km": h_leo_target / 1000, "v_launch_absolute_mps": v_launch_absolute, "v_launcher_relative_to_ground_mps": v_launcher_relative_to_ground, "v_at_geo_apogee_initial_ellipse_mps": v_at_geo_apogee_initial_ellipse, "delta_v_burn1_at_geo_mps": delta_v_burn1_at_geo, # This will be positive for prograde burn "v_at_leo_perigee_transfer_ellipse_mps": v_transfer_perigee, "delta_v_burn2_at_leo_mps": delta_v_burn2_at_leo, # This will be negative for retrograde burn "v_circular_leo_mps": v_circular_leo, "total_delta_v_burns_mps": abs(delta_v_burn1_at_geo) + abs(delta_v_burn2_at_leo), "a_initial_ellipse_km": a_initial_ellipse / 1000, "e_initial_ellipse": eccentricity_initial_ellipse, "nu_0_initial_ellipse_degrees": math.degrees(nu_0_initial_ellipse), "a_transfer_ellipse_km": a_transfer_ellipse / 1000, "e_transfer_ellipse": eccentricity_transfer_ellipse, } print("\n--- Multi-Burn Orbital Calculation Results ---") for key, value in results.items(): if isinstance(value, (int, float)): print(f"{key.replace('_', ' ').title()}: {value:.3f}") # More precision for velocities else: print(f"{key.replace('_', ' ').title()}: {value}") root.mainloop() # Start the Tkinter event loop # Run the function to generate data and display the plot if __name__ == "__main__": generate_orbital_data_and_plot() # 341 lines in Version 4
The program will show calculation results in the terminal window.
To close the program click the "Close" button in upper left.
(th)
Offline
Like button can go here