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
This post contains a Python program developed by Google's Gemini (2.5 Flash) to help plan a flight of one of Dr. John Hunter's gas guns. In the scenario for which the program is designed, the air ahead of the projectile will be collected as payload. In other words, drag is converted to payload mass (by some means not yet known). the mass of the projectile increases as it ascends and collects molecules from the atmosphere. the force of gravity decreases but not by much. On the other hand, the mass of the air collected also decreases on the way up. this program is designed to attempt to deal with the complexity of the scenario.
import math
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m) - a good approximation for lower atmosphere
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
This model is a reasonable approximation for the altitudes of interest.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def simulate_air_harvesting_launch(
projectile_diameter_m=1.0,
initial_projectile_mass_kg=1000.0,
initial_velocity_mps=11000.0,
launch_altitude_meters_above_sea_level=0.0,
final_altitude_meters_above_sea_level=100000.0,
time_step_s=0.01
):
"""
Simulates a projectile launch from a gas gun, modeling the deceleration
from gravity and the work done to collect and accelerate atmospheric mass.
Args:
projectile_diameter_m (float): Diameter of the projectile in meters.
initial_projectile_mass_kg (float): Starting mass of the projectile in kilograms.
initial_velocity_mps (float): Initial launch velocity in meters per second.
launch_altitude_meters_above_sea_level (float): Starting altitude.
final_altitude_meters_above_sea_level (float): The target altitude for the simulation.
time_step_s (float): The time step for the numerical integration.
"""
# --- Input Validation ---
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values (velocity, diameter, mass) must be positive.")
return
# --- Initial Conditions ---
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
print(f"--- Gas Gun Air Harvesting Simulation ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print("------------------------------------------")
# --- Main Simulation Loop ---
while current_altitude_m < final_altitude_meters_above_sea_level:
# Check for projectile slowing to a stop or going backward
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
# Calculate forces and mass changes at the current state
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
# Calculate force from gravity
force_gravity = current_mass_kg * g
# Calculate the rate of mass collection and "drag" force from momentum transfer
# This is the core assumption of the "magic pipe" - all mass is collected and accelerated.
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
# Calculate work done in this time step
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
# Note: Using a simplified Euler integration for demonstration purposes.
# It's an approximation but sufficient for this conceptual model.
net_force = -force_gravity - force_air_collection
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_mass_kg += dm_dt * time_step_s
# --- Final Output ---
final_mass_of_collected_air_kg = current_mass_kg - initial_projectile_mass_kg
print("--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {final_mass_of_collected_air_kg:.2f} kg ({final_mass_of_collected_air_kg / 1000:.2f} metric tons)")
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
# For comparison, show the initial kinetic energy
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
# --- Example Usage ---
# Default values based on our discussion
simulate_air_harvesting_launch(
projectile_diameter_m=1.0,
initial_projectile_mass_kg=1000.0, # 1 metric ton payload
initial_velocity_mps=11000.0, # This is a high-end velocity, close to escape velocity
launch_altitude_meters_above_sea_level=0.0,
final_altitude_meters_above_sea_level=100000.0
)This program should run on any modern computer with one of the three major OS: Linux, Windows, Apple.
The program will run on Android if Linux is installed as a subsystem.
Note: the input values are defined in the statement at the bottom of the code. This program can be enhanced with input statements.
(th)
Offline
Like button can go here
Here is a preliminary output from the utility in Post #57
--- Gas Gun Air Harvesting Simulation ---
Projectile Diameter: 1.0 m
Initial Mass: 1000.0 kg
Initial Velocity: 11000.0 m/s (11.00 km/s)
Target Altitude: 100 km
------------------------------------------
Projectile stalled at 77.35 km. Velocity is zero or negative.
--- Simulation Complete ---
Final Altitude: 77.35 km
Final Velocity: -0.09 m/s (-0.00 km/s)
Mass of Collected Air: 8269.73 kg (8.27 metric tons)
Total Work Done on Air Mass: 104.39 GigaJoules
Total Work Done Against Gravity: 6.28 GigaJoules
Initial Kinetic Energy: 60.50 GigaJoules
The simulation is for a gas gun firing straight up from the surface of the Earth.
The scenario is capture of the air in the column above the projectile.
The capture scenario is an alternative to the normal practice of pushing the air aside.
The air above the projectile totals (about) 10.34 metric tons per square meter.
The project starting mass will increase by that amount as it ascends, in this simulation.
(th)
Offline
Like button can go here
In this post we have an update to the Python program to simulate collection of atmosphere above a projectile.
In this version, we have added input opportunities for the operator, and a parameter file that persists between runs.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
time_step_s
):
"""
Simulates a projectile launch from a gas gun, modeling the deceleration
from gravity and the work done to collect and accelerate atmospheric mass.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values (velocity, diameter, mass) must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
print(f"--- Gas Gun Air Harvesting Simulation ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
total_work_done_against_gravity_J += work_gravity_step
net_force = -force_gravity - force_air_collection
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_mass_kg += dm_dt * time_step_s
final_mass_of_collected_air_kg = current_mass_kg - initial_projectile_mass_kg
print("--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {final_mass_of_collected_air_kg:.2f} kg ({final_mass_of_collected_air_kg / 1000:.2f} metric tons)")
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"time_step_s": 0.01
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func=float):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
return type_func(user_input)
except ValueError:
print("Invalid input. Please enter a valid number.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"])
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"])
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"])
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"])
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"])
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
# Run the simulation with the final parameters
simulate_air_harvesting_launch(**params)
# Run the main function
if __name__ == "__main__":
main()(th)
Offline
Like button can go here
In this post, the code block will contain the next iteration of the Air Collection system. In this version, we add a front door to the "magic pipe" that collects air on the way up. We offer the user the opportunity to turn the switch on or off, in order to compare performance. The door will automatically close when the pressure inside the pipe matches the dynamic pressure of the air ahead of the pipe. At that point, the mass of the projectile will become constant, and the drag will become significant. However, the drag will decrease with altitude and with the falling velocity of the projectile. I am expecting the energy cost of the flight to decrease when the door closes, but have no idea how much.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K=288.0):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This is an approximation and assumes a constant temperature.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.028 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
time_step_s,
use_door_closing
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
print(f"--- Gas Gun Air Harvesting Simulation ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
if not doors_closed:
# --- Phase 1: Open Doors ---
# Check if doors should close (internal pressure >= dynamic pressure)
if use_door_closing:
collected_air_mass = current_mass_kg - initial_projectile_mass_kg
volume_of_pipe = cross_sectional_area_sqm * current_altitude_m
internal_pressure = get_gas_pressure(collected_air_mass, volume_of_pipe)
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if internal_pressure >= dynamic_pressure and collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
# In this phase, we use traditional drag from now on.
# The momentum transfer drag is now replaced by traditional drag.
if not doors_closed:
# Calculate collection force and add to net force
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
net_force -= force_air_collection
# Work done on air
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
# Update mass
current_mass_kg += dm_dt * time_step_s
else:
# --- Phase 2: Closed Doors ---
# Now we have fixed mass and traditional drag
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# Work done by drag is dissipative, so we won't add it to a 'useful' work counter
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
final_mass_of_collected_air_kg = current_mass_kg - initial_projectile_mass_kg
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {final_mass_of_collected_air_kg:.2f} kg ({final_mass_of_collected_air_kg / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"time_step_s": 0.01,
"use_door_closing": False
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
return type_func(user_input)
except ValueError:
print("Invalid input. Please enter a valid number or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
# New input for the door closing feature
while True:
use_door_input = input(f"Use door closing feature? (True/False) (default: {params['use_door_closing']}): ")
if use_door_input == "":
break
if use_door_input.lower() in ('true', 'false'):
params['use_door_closing'] = use_door_input.lower() == 'true'
break
else:
print("Invalid input. Please enter 'True' or 'False'.")
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()As always, anyone can run this program on a modern computer using Linux, Windows or Apple operating systems. The program will run under Android (such as Chromebook) if you have Linux installed as an app. Warning! Please do NOT try to add Linux if your Chromebook has less than 8 Gb.
(th)
Offline
Like button can go here
This post contains an update to the "Air Capture" simulation.
This version includes a parameter to define the volume into which we will be stuffing air on the way up.
In the earlier versions of the program, the volume of the collection tank was assumed to be infinite. That was a useful device to get the program working, but it is now time to impose a limit on the volume of the collection tank.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K=288.0):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This is an approximation and assumes a constant temperature.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
collected_air_mass = 0.0
print(f"--- Gas Gun Air Harvesting Simulation ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
if not doors_closed:
# --- Phase 1: Open Doors ---
# Check if doors should close (internal pressure >= dynamic pressure)
if use_door_closing and collected_air_mass > 0:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if internal_pressure >= dynamic_pressure:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
if not doors_closed:
# Calculate collection force and update mass
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
net_force -= force_air_collection
# Work done on air
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
# Update mass
mass_collected_in_step = dm_dt * time_step_s
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
else:
# --- Phase 2: Closed Doors ---
# Now we have fixed mass and traditional drag
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0, # A new default based on previous calculations
"time_step_s": 0.01,
"use_door_closing": False
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
while True:
use_door_input = input(f"Use door closing feature? (True/False) (default: {params['use_door_closing']}): ")
if use_door_input == "":
break
if use_door_input.lower() in ('true', 'false'):
params['use_door_closing'] = use_door_input.lower() == 'true'
break
else:
print("Invalid input. Please enter 'True' or 'False'.")
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()The default value is the volume of the 10.34 tons of air above a square meter, if it were cooled to liquid.
(th)
Offline
Like button can go here
In this version of the Air Capture scenario, we improve performance of the "Doors Open" option.
The gas collected during the early phase of flight will push against the oncoming gas and flow out of the container as dynamic pressure decreases. This may allow for deflection of oncoming gas to the side as a strategy for dealing with the gas ahead of the projectile. At the end of the flight the container would be empty, so the purpose of the exercise is to see if it has protected the projectile from the rigors of flight through the atmosphere. It is possible this version does not correctly reduce the mass of the projectile as captured gas escapes. What I'm hoping for with this version is improved management of the forces in play. The gas escaping from the capture container will cause a negative thrust against the forward momentum of the projectile.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K=288.0):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This is an approximation and assumes a constant temperature.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_full = False
collected_air_mass = 0.0
print(f"--- Gas Gun Air Harvesting Simulation ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
if use_door_closing:
# --- Scenario 1: Closed Doors ---
if not doors_closed:
# Check if doors should close (internal pressure >= dynamic pressure)
if collected_air_mass > 0:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if internal_pressure >= dynamic_pressure:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
if not doors_closed:
# Phase 1: Collect mass, increase projectile mass
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
mass_collected_in_step = dm_dt * time_step_s
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
else:
# Phase 2: Fixed mass, traditional drag
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
else:
# --- Scenario 2: Open Doors (New Logic) ---
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
if not pipe_full:
# Phase 1: Fill the pipe
mass_to_fill = (pipe_volume_m3 * rho) - collected_air_mass
if mass_to_fill > dm_dt_in * time_step_s:
mass_collected_in_step = dm_dt_in * time_step_s
force_air_collection = dm_dt_in * current_velocity_mps
else:
mass_collected_in_step = mass_to_fill
force_air_collection = dm_dt_in * current_velocity_mps
pipe_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at fill: {collected_air_mass + mass_to_fill:.2f} kg.")
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
else:
# Phase 2: Pipe is full, mass inflow = mass outflow (fixed mass)
# This is an idealized representation of the complex reality you described
# where outflow causes a braking force.
# Net force is from gravity plus the dynamic pressure of the incoming air.
# Mass inflow and outflow are equal, so total mass is constant.
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": False
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
while True:
use_door_input = input(f"Use door closing feature? (True/False) (default: {params['use_door_closing']}): ")
if use_door_input == "":
break
if use_door_input.lower() in ('true', 'false'):
params['use_door_closing'] = use_door_input.lower() == 'true'
break
else:
print("Invalid input. Please enter 'True' or 'False'.")
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()This program can be run by anyone on a modern computer, with Linux, Windows or Apple. Android can be given a Linux app that will work if you have 8 Gb on your Android device.
(th)
Offline
Like button can go here
In this version we concentrate on improving performance of the "open door" flow.
If a reader is comparing versions, we are now correcting some settings that were temporary fillers during earlier phases of development.
This latest version attempts to model the pressure of air at the entrance to the "magic pipe" more accurately, and it attempts to model the loss of mass as the projectile ascends, due to the greater pressure inside the capture chamber compared to the dynamic pressure at the entrance to the chamber.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K=288.0):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This is an approximation and assumes a constant temperature.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
print(f"--- Gas Gun Air Harvesting Simulation V6 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
if use_door_closing:
# --- Scenario 1: Closed Doors ---
if not doors_closed:
# Phase 1: Collect mass, increase projectile mass
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if internal_pressure >= dynamic_pressure and collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
if not doors_closed:
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
mass_collected_in_step = dm_dt * time_step_s
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
else:
# Phase 2: Fixed mass, traditional drag
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
else:
# --- Scenario 2: Open Doors (New Logic) ---
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
if not pipe_is_full:
# Phase 1: Fill the pipe until mass_collected equals pipe_volume * rho_at_fill
mass_needed_to_fill = (pipe_volume_m3 * rho)
mass_collected_in_step = dm_dt_in * time_step_s
if collected_air_mass + mass_collected_in_step >= mass_needed_to_fill:
mass_to_fill_remainder = mass_needed_to_fill - collected_air_mass
current_mass_kg += mass_to_fill_remainder
collected_air_mass += mass_to_fill_remainder
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
else:
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
else:
# Phase 2 & 3: Outflow and eventual empty pipe
if not pipe_is_empty:
# Mass outflow creates retro-thrust
outflow_rate = outflow_rate_factor * (get_gas_pressure(collected_air_mass, pipe_volume_m3) - (0.5 * rho * current_velocity_mps**2))
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
# Retro-thrust force
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
# Mass of projectile decreases
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
# Inflow is still happening, creating a drag force.
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": False,
"outflow_rate_factor": 1.0e-5
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()This program will run on any modern computer with Linux, Windows or Apple.
(the)
Offline
Like button can go here
In the version of the Air Capture Simulation #7 below, we are attempting to correct leftover code from early development. If anyone is following progress, they will see that initialization of the procedure for filling the capture container has been revised. The behavior of the two flows should be identical right up to the instant the door at the top of the intake closes. In Version 7, Gemini is hoping that is accomplished, but we won't be sure until the next test results are in. I am ** very ** interested in seeing how the alternative paths evolve with this change. The two modes of operation might be selected for different situations.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K=288.0):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This is an approximation and assumes a constant temperature.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
print(f"--- Gas Gun Air Harvesting Simulation V7 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
force_gravity = current_mass_kg * g
net_force = -force_gravity
if use_door_closing:
# --- Scenario 1: Closed Doors ---
if not doors_closed:
# Phase 1: Collect mass until internal pressure >= dynamic pressure
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
if internal_pressure >= dynamic_pressure and collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
if not doors_closed:
dm_dt = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
mass_collected_in_step = dm_dt * time_step_s
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
else:
# Phase 2: Fixed mass, traditional drag
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
else:
# --- Scenario 2: Open Doors (Corrected Logic) ---
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
if not pipe_is_full:
# Phase 1: Fill the pipe until internal pressure >= dynamic pressure
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
if internal_pressure >= dynamic_pressure and collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
if not pipe_is_full:
mass_collected_in_step = dm_dt_in * time_step_s
force_air_collection = dm_dt_in * current_velocity_mps
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
else:
# Phase 2 & 3: Outflow and eventual empty pipe
if not pipe_is_empty:
# Mass outflow creates retro-thrust
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
# Retro-thrust force
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
# Mass of projectile changes
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
# Inflow is still happening, creating a drag force.
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": False,
"outflow_rate_factor": 1.0e-5
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()(th)
Offline
Like button can go here
This post contains Air Capture Simulation #8
I asked Gemini to add a data output to the display. In the course of trying to do that, Gemini found an error which it corrected.
The program in the block below ** should ** report temperature as well as pressure. We'll see.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This now requires a temperature input.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
# We now track the gas temperature
gas_temperature_K = TEMP_0_K
print(f"--- Gas Gun Air Harvesting Simulation V8 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
# Calculate dynamic pressure and ambient pressure for the adiabatic compression model
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
ambient_pressure_at_altitude = PRESS_0_PA * math.exp(-current_altitude_m / H)
# --- PHASE 1: COLLECTING AIR ---
if not doors_closed and not pipe_is_full:
# We now use adiabatic compression to calculate temperature and pressure
# Use the new mass to find the new pressure
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
# Update collected mass
temp_collected_air_mass = collected_air_mass + mass_collected_in_step
# Use the new pressure and mass to find the new temperature
# This is an iterative process, so we use the previous step's values as P1 and T1
if temp_collected_air_mass > 0:
# Calculate new internal pressure based on the new mass
temp_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K)
# Use adiabatic compression formula to update temperature based on the new pressure
gas_temperature_K = TEMP_0_K * (temp_pressure / PRESS_0_PA)**((GAMMA - 1) / GAMMA)
# Now, use the corrected temperature to find the final pressure for this step
internal_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K)
else:
internal_pressure = 0.0
# This is where the two tracks diverge
if use_door_closing:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
# --- PHASE 2: CLOSED DOORS ---
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# --- PHASE 3: OPEN DOORS (OUTFLOW) ---
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# --- PHASE 4: EMPTY PIPE, OPEN DOORS ---
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg # mass reverts to original, no more collection
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": False,
"outflow_rate_factor": 1.0e-5
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()This program will run on any modern computer using Linux, Windows or Apple.
If someone has a Chromebook with 8 Gb, then the Linux app can be installed.
(th)
Offline
Like button can go here
This post will show an inspection program that ChatGPT5 developed for the ongoing Merlin Engine update that has been in progress for months.
# inspect_face_typesV08.py
# Fully robust inspector with real lookup of owner and neighbour relationships
import json
import re
def load_openfoam_file(filepath):
try:
with open(filepath, 'r') as f:
lines = f.readlines()
open_idx = None
close_idx = None
for i, line in enumerate(lines):
if line.strip() == "(":
open_idx = i
break
for i in range(len(lines)-1, -1, -1):
if lines[i].strip() == ")":
close_idx = i
break
if open_idx is None or close_idx is None:
print(f"❌ Could not find block in {filepath}")
return {}
count_line = lines[open_idx - 1]
match = re.match(r"\s*(\d+)\s*", count_line)
if match:
count = int(match.group(1))
else:
print(f"❌ Could not find count in {filepath}")
return {}
data_lines = [line.strip() for line in lines[open_idx+1:close_idx]]
if len(data_lines) != count:
print(f"⚠️ WARNING: {filepath} declared {count}, found {len(data_lines)})")
else:
print(f"✅ {filepath}: {count} entries loaded correctly.")
return {
"count": count,
"data_lines": data_lines,
}
except FileNotFoundError:
print(f"❌ File not found: {filepath}")
return {}
def inspect_face_types(cells_file="persistent_memory_cells.json", owner_file="owner", neighbour_file="neighbour"):
with open(cells_file, 'r') as f:
cell_data = json.load(f)
owner = load_openfoam_file(owner_file)
neighbour = load_openfoam_file(neighbour_file)
owner_list = [int(x) for x in owner["data_lines"]]
neighbour_list = [int(x) for x in neighbour["data_lines"]]
# Create lookup maps
face_to_owner = {i: cell_id for i, cell_id in enumerate(owner_list)}
face_to_neighbour = {i: cell_id for i, cell_id in enumerate(neighbour_list)} # internal only
total_checked = 0
total_matched = 0
total_mismatched = 0
missing_from_owner = 0
internal_count = 0
boundary_count = 0
mismatch_log = []
for cell_name, cell in cell_data.items():
cell_id = int(cell_name.split("_")[1])
for face_label, face in cell["faces"].items():
face_id = face["face_id"]
face_type = face["type"]
total_checked += 1
if face_type == "internal":
internal_count += 1
elif face_type == "boundary":
boundary_count += 1
if face_id not in face_to_owner:
print(f"⚠️ Face ID {face_id} missing from owner file")
missing_from_owner += 1
continue
# Determine actual connectivity type
is_internal = face_id in face_to_neighbour
expected_type = "internal" if is_internal else "boundary"
# Check if this cell is one of the owners
owner_id = face_to_owner[face_id]
neighbour_id = face_to_neighbour.get(face_id, None)
cell_is_connected = (cell_id == owner_id) or (cell_id == neighbour_id)
if not cell_is_connected:
print(f"⚠️ Cell ID {cell_id} not connected to face ID {face_id} (owner: {owner_id}, neighbour: {neighbour_id})")
if face_type == expected_type and cell_is_connected:
total_matched += 1
else:
total_mismatched += 1
mismatch_log.append({
"face_id": face_id,
"label": face_label,
"cell": cell_name,
"expected": expected_type,
"actual": face_type
})
# Summary
print("✅ Face Type Inspection Complete")
print(f" Total faces checked: {total_checked}")
print(f" Matches: {total_matched}")
print(f" Mismatches: {total_mismatched}")
print(f" Missing from owner: {missing_from_owner}")
print(f" Face type tally from memory:")
print(f" Boundary: {boundary_count}")
print(f" Internal: {internal_count}")
if mismatch_log:
with open("face_type_mismatches.txt", "w") as out:
for entry in mismatch_log:
out.write(
f"Face ID {entry['face_id']} ({entry['label']}) in {entry['cell']}:\n"
f" Expected: {entry['expected']} | Actual: {entry['actual']}\n\n"
)
print("⚠️ Mismatch details written to face_type_mismatches.txt")
if __name__ == "__main__":
inspect_face_types() # 139 lines in V08(th)
Offline
Like button can go here
This post contains Version 9 of the Air Capture Simulation program:
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters):
"""
Calculates air density using a simplified exponential model.
"""
return RHO_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K):
"""
Estimates pressure of collected gas using Ideal Gas Law.
This now requires a temperature input.
"""
# Assuming the gas is mostly N2, the molar mass is ~0.02896 kg/mol
molar_mass = 0.02896
R = 8.314 # Universal gas constant
n = mass_kg / molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0 # New variable to track elapsed time
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
# We now track the gas temperature
gas_temperature_K = TEMP_0_K
print(f"--- Gas Gun Air Harvesting Simulation V9 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m)
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
# --- PHASE 1: COLLECTING AIR ---
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
temp_collected_air_mass = collected_air_mass + mass_collected_in_step
if temp_collected_air_mass > 0:
# Calculate new internal pressure based on the new mass
temp_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K)
# Use adiabatic compression formula to update temperature based on the new pressure
gas_temperature_K = TEMP_0_K * (temp_pressure / PRESS_0_PA)**((GAMMA - 1) / GAMMA)
# Now, use the corrected temperature to find the final pressure for this step
internal_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K)
else:
internal_pressure = 0.0
if use_door_closing:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.") # New time reporting
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.") # New time reporting
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
# --- PHASE 2: CLOSED DOORS ---
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# --- PHASE 3: OPEN DOORS (OUTFLOW) ---
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# --- PHASE 4: EMPTY PIPE, OPEN DOORS ---
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg # mass reverts to original, no more collection
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s # Increment elapsed time
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
print(f"Total Time: {current_time_s:.2f} s") # New total time reporting
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": False,
"outflow_rate_factor": 1.0e-5
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()The program will run on any modern computer with Linux, Windows or Apple. Chromebook can run the Linux addon with 8 Gb of RAM.
(th)
Offline
Like button can go here
I am pleased to be able to offer Version 10 of the Air Capture Simulation program that Gemini and I've been working on.
This one adds a new feature to account for humidity in the air above the launcher. The program will include an input to allow the operator to set humidity at a value between 0 and 100. I assume there may be a need for a setting above 100% if the launch occurs during a rain storm, but for this first version with the new feature, just working with humidity seems like a major step forward. One detail that caught my eye as we discussed this new feature, is that density of humid air is slightly less than dry air. The connection I make is with weather patterns that feature low pressure air and high pressure air. I deduce (subject to correction) that low pressure areas are more humid than high pressure areas. My assumption is that Dr. Hunter's launcher will be located in the ocean when it reaches full size, and the air above the launcher will always be humid to some extent. Thus, the payload delivered to LEO would include hydrogen as well as oxygen and nitrogen.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
# Molar masses (kg/mol)
MOLAR_MASS_DRY_AIR = 0.028965
MOLAR_MASS_WATER = 0.018015
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_air_density(altitude_meters, relative_humidity_percent=0.0):
"""
Calculates moist air density using a simplified exponential model and humidity.
Relative humidity is a percentage from 0 to 100.
"""
# Calculate saturation vapor pressure at sea-level temperature
# Using the Arden Buck equation for T in Celsius
temp_celsius = TEMP_0_K - 273.15
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius))) # kPa
es_pa = es * 1000.0 # Pa
# Actual vapor pressure
e_pa = es_pa * (relative_humidity_percent / 100.0)
# Air density at sea level, adjusted for humidity
# See https://en.wikipedia.org/wiki/Density_of_air#Humid_air
rho_moist_0 = (PRESS_0_PA / (287.058 * TEMP_0_K)) * (1 - e_pa / PRESS_0_PA * (1 - MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR))
# Apply exponential decay for altitude
return rho_moist_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K, relative_humidity_percent=0.0):
"""
Estimates pressure of collected gas using Ideal Gas Law, adjusted for molar mass.
"""
# Calculate average molar mass of the collected gas mixture
# We'll use the ratio of partial pressures to get the molar fraction
# Calculate saturation vapor pressure at sea-level temperature
temp_celsius = TEMP_0_K - 273.15
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius))) # kPa
es_pa = es * 1000.0 # Pa
e_pa = es_pa * (relative_humidity_percent / 100.0)
# Assuming constant temperature at collection, pressure ratio is mole fraction
# This is a simplification but works for this model
mole_fraction_water = e_pa / PRESS_0_PA
mole_fraction_dry_air = 1.0 - mole_fraction_water
avg_molar_mass = (mole_fraction_dry_air * MOLAR_MASS_DRY_AIR) + (mole_fraction_water * MOLAR_MASS_WATER)
R = 8.314 # Universal gas constant
n = mass_kg / avg_molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor,
relative_humidity_percent
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
# We now track the gas temperature
gas_temperature_K = TEMP_0_K
print(f"--- Gas Gun Air Harvesting Simulation V10 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Relative Humidity: {relative_humidity_percent:.1f}%")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m, relative_humidity_percent) # Use new function
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
# --- PHASE 1: COLLECTING AIR ---
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
temp_collected_air_mass = collected_air_mass + mass_collected_in_step
if temp_collected_air_mass > 0:
# Calculate new internal pressure based on the new mass
temp_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
# Use adiabatic compression formula to update temperature based on the new pressure
gas_temperature_K = TEMP_0_K * (temp_pressure / PRESS_0_PA)**((GAMMA - 1) / GAMMA)
# Now, use the corrected temperature to find the final pressure for this step
internal_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
else:
internal_pressure = 0.0
if use_door_closing:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
# --- PHASE 2: CLOSED DOORS ---
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# --- PHASE 3: OPEN DOORS (OUTFLOW) ---
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# --- PHASE 4: EMPTY PIPE, OPEN DOORS ---
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg # mass reverts to original, no more collection
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
print(f"Total Time: {current_time_s:.2f} s")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": True,
"outflow_rate_factor": 1.0e-5,
"relative_humidity_percent": 0.0
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
params["relative_humidity_percent"] = get_user_input("Enter new relative humidity (%)", params["relative_humidity_percent"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()(th)
Offline
Like button can go here
In Version 11 of the Air Capture Simulation, Gemini has added estimates of the amount of oxygen, nitrogen and water collected by the projectile. As a reminder, the capture only works if the doors are closed after pressure inside the container equals the pressure outside the rising projectile.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
# Molar masses (kg/mol)
MOLAR_MASS_DRY_AIR = 0.028965
MOLAR_MASS_WATER = 0.018015
# Mass fractions of dry air
MASS_FRACTION_NITROGEN = 0.755
MASS_FRACTION_OXYGEN = 0.232
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_vapor_pressure(temp_celsius):
"""Calculates saturation vapor pressure in Pa using the Arden Buck equation."""
# Arden Buck equation for T in Celsius
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius))) # kPa
return es * 1000.0 # Pa
def get_air_density(altitude_meters, relative_humidity_percent=0.0):
"""
Calculates moist air density using a simplified exponential model and humidity.
Relative humidity is a percentage from 0 to 100.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
# Air density at sea level, adjusted for humidity
rho_moist_0 = (PRESS_0_PA / (287.058 * TEMP_0_K)) * (1 - e_pa / PRESS_0_PA * (1 - MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR))
# Apply exponential decay for altitude
return rho_moist_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K, relative_humidity_percent=0.0):
"""
Estimates pressure of collected gas using Ideal Gas Law, adjusted for molar mass.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
mole_fraction_water = e_pa / PRESS_0_PA
mole_fraction_dry_air = 1.0 - mole_fraction_water
avg_molar_mass = (mole_fraction_dry_air * MOLAR_MASS_DRY_AIR) + (mole_fraction_water * MOLAR_MASS_WATER)
R = 8.314 # Universal gas constant
n = mass_kg / avg_molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
initial_projectile_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor,
relative_humidity_percent
):
"""
Simulates a projectile launch from a gas gun with optional door closing.
"""
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
gas_temperature_K = TEMP_0_K
print(f"--- Gas Gun Air Harvesting Simulation V11 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Relative Humidity: {relative_humidity_percent:.1f}%")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m, relative_humidity_percent)
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
# --- PHASE 1: COLLECTING AIR ---
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
temp_collected_air_mass = collected_air_mass + mass_collected_in_step
if temp_collected_air_mass > 0:
temp_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
gas_temperature_K = TEMP_0_K * (temp_pressure / PRESS_0_PA)**((GAMMA - 1) / GAMMA)
internal_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
else:
internal_pressure = 0.0
if use_door_closing:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
# --- PHASE 2: CLOSED DOORS ---
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
# --- PHASE 3: OPEN DOORS (OUTFLOW) ---
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
# --- PHASE 4: EMPTY PIPE, OPEN DOORS ---
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg
# Work done against gravity
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
# Update state variables
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Mass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
print(f"Total Time: {current_time_s:.2f} s")
if use_door_closing:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
else:
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
initial_ke = 0.5 * initial_projectile_mass_kg * initial_velocity_mps**2
print(f"Initial Kinetic Energy: {initial_ke/1e9:.2f} GigaJoules")
# --- New Composition Analysis Section ---
if collected_air_mass > 0:
print("\n--- Payload Composition (by mass) ---")
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
# Assuming gas law applies, mole fractions are proportional to partial pressures
mass_of_water = (e_pa / PRESS_0_PA) * (MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR) * collected_air_mass
mass_of_dry_air = collected_air_mass - mass_of_water
# Recalculate based on a more accurate mass fraction of water
rho_moist_0 = get_air_density(0, relative_humidity_percent)
rho_water_vapor_0 = (e_pa / (461.5 * TEMP_0_K))
mass_of_water = collected_air_mass * (rho_water_vapor_0 / rho_moist_0)
mass_of_dry_air = collected_air_mass - mass_of_water
mass_of_nitrogen = mass_of_dry_air * MASS_FRACTION_NITROGEN
mass_of_oxygen = mass_of_dry_air * MASS_FRACTION_OXYGEN
# Calculate percentages
percent_water = (mass_of_water / collected_air_mass) * 100.0
percent_nitrogen = (mass_of_nitrogen / collected_air_mass) * 100.0
percent_oxygen = (mass_of_oxygen / collected_air_mass) * 100.0
print(f"Water: {mass_of_water:.2f} kg ({percent_water:.2f}%)")
print(f"Nitrogen: {mass_of_nitrogen:.2f} kg ({percent_nitrogen:.2f}%)")
print(f"Oxygen: {mass_of_oxygen:.2f} kg ({percent_oxygen:.2f}%)")
# The remainder is other dry air components (Argon, etc.)
other_mass = collected_air_mass - mass_of_water - mass_of_nitrogen - mass_of_oxygen
other_percent = (other_mass / collected_air_mass) * 100.0
print(f"Other: {other_mass:.2f} kg ({other_percent:.2f}%)")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"initial_projectile_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": True,
"outflow_rate_factor": 1.0e-5,
"relative_humidity_percent": 0.0
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["initial_projectile_mass_kg"] = get_user_input("Enter new initial projectile mass (kg)", params["initial_projectile_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
params["relative_humidity_percent"] = get_user_input("Enter new relative humidity (%)", params["relative_humidity_percent"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()feedback is welcome. This is an investigation of a possible way of improving commercial value of a gun launch system.
(th)
Offline
Like button can go here
This post contains Version 12 of the Air Capture Simulation program.
Advice: This program needs to create a new parameter store. If you have run a previous version, please deleted the file:
simulation_params.json
*** added 2025/08/24 after a test run
This version adds inputs for the three major subsystems of the projectile.
The first section is the Air Capture container which must include a doors mechanism. That is not yet defined.
The middle section is the payload and control section.
The third and essential section is a solid fuel rocket with a range of capability from 8 km/s for LEO to 3 km/s for GEO.
No other changes were made in this version.
As a side note: Gemini advises against running this system in rain. The solid (incompressible) form of water will destroy the system.
On the other hand, running with 100% humidity is fine, because the water in the air is a gas.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
# Molar masses (kg/mol)
MOLAR_MASS_DRY_AIR = 0.028965
MOLAR_MASS_WATER = 0.018015
# Mass fractions of dry air
MASS_FRACTION_NITROGEN = 0.755
MASS_FRACTION_OXYGEN = 0.232
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_vapor_pressure(temp_celsius):
"""Calculates saturation vapor pressure in Pa using the Arden Buck equation."""
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius)))
return es * 1000.0
def get_air_density(altitude_meters, relative_humidity_percent=0.0):
"""
Calculates moist air density using a simplified exponential model and humidity.
Relative humidity is a percentage from 0 to 100.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = (PRESS_0_PA / (287.058 * TEMP_0_K)) * (1 - e_pa / PRESS_0_PA * (1 - MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR))
return rho_moist_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K, relative_humidity_percent=0.0):
"""
Estimates pressure of collected gas using Ideal Gas Law, adjusted for molar mass.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
mole_fraction_water = e_pa / PRESS_0_PA
mole_fraction_dry_air = 1.0 - mole_fraction_water
avg_molar_mass = (mole_fraction_dry_air * MOLAR_MASS_DRY_AIR) + (mole_fraction_water * MOLAR_MASS_WATER)
R = 8.314
n = mass_kg / avg_molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
pipe_mass_kg,
payload_mass_kg,
rocket_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor,
relative_humidity_percent
):
"""
Simulates a projectile launch with subsystem mass inputs.
"""
initial_projectile_mass_kg = pipe_mass_kg + payload_mass_kg + rocket_mass_kg
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
gas_temperature_K = TEMP_0_K
print(f"--- Gas Gun Air Harvesting Simulation V12 ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Pipe Mass: {pipe_mass_kg} kg")
print(f"Payload Mass: {payload_mass_kg} kg")
print(f"Rocket Mass: {rocket_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Relative Humidity: {relative_humidity_percent:.1f}%")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m, relative_humidity_percent)
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
temp_collected_air_mass = collected_air_mass + mass_collected_in_step
if temp_collected_air_mass > 0:
temp_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
gas_temperature_K = TEMP_0_K * (temp_pressure / PRESS_0_PA)**((GAMMA - 1) / GAMMA)
internal_pressure = get_gas_pressure(temp_collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
else:
internal_pressure = 0.0
if use_door_closing:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and temp_collected_air_mass > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Total Time: {current_time_s:.2f} s")
print(f"Initial Kinetic Energy: {0.5 * initial_projectile_mass_kg * initial_velocity_mps**2 / 1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
print("\n--- Mass Breakdown at Closure ---")
print(f"Total Projectile Mass: {initial_projectile_mass_kg:.2f} kg")
print(f"Magic Pipe Mass: {pipe_mass_kg:.2f} kg")
print(f"Payload Mass: {payload_mass_kg:.2f} kg")
print(f"Rocket Mass: {rocket_mass_kg:.2f} kg")
print(f"\nMass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if collected_air_mass > 0:
print("\n--- Payload Composition (by mass) ---")
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = get_air_density(0, relative_humidity_percent)
rho_water_vapor_0 = (e_pa / (461.5 * TEMP_0_K))
mass_of_water = collected_air_mass * (rho_water_vapor_0 / rho_moist_0)
mass_of_dry_air = collected_air_mass - mass_of_water
mass_of_nitrogen = mass_of_dry_air * MASS_FRACTION_NITROGEN
mass_of_oxygen = mass_of_dry_air * MASS_FRACTION_OXYGEN
percent_water = (mass_of_water / collected_air_mass) * 100.0
percent_nitrogen = (mass_of_nitrogen / collected_air_mass) * 100.0
percent_oxygen = (mass_of_oxygen / collected_air_mass) * 100.0
print(f"Water: {mass_of_water:.2f} kg ({percent_water:.2f}%)")
print(f"Nitrogen: {mass_of_nitrogen:.2f} kg ({percent_nitrogen:.2f}%)")
print(f"Oxygen: {mass_of_oxygen:.2f} kg ({percent_oxygen:.2f}%)")
other_mass = collected_air_mass - mass_of_water - mass_of_nitrogen - mass_of_oxygen
other_percent = (other_mass / collected_air_mass) * 100.0
print(f"Other: {other_mass:.2f} kg ({other_percent:.2f}%)")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
if os.path.exists(PARAMS_FILE):
with open(PARAMS_FILE, 'r') as f:
return json.load(f)
else:
return {
"projectile_diameter_m": 1.0,
"pipe_mass_kg": 1000.0,
"payload_mass_kg": 1000.0,
"rocket_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": True,
"outflow_rate_factor": 1.0e-5,
"relative_humidity_percent": 0.0
}
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["pipe_mass_kg"] = get_user_input("Enter new magic pipe mass (kg)", params["pipe_mass_kg"], float)
params["payload_mass_kg"] = get_user_input("Enter new payload mass (kg)", params["payload_mass_kg"], float)
params["rocket_mass_kg"] = get_user_input("Enter new rocket mass (kg)", params["rocket_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
params["relative_humidity_percent"] = get_user_input("Enter new relative humidity (%)", params["relative_humidity_percent"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()(th)
Offline
Like button can go here
This post contains Version 14 of the Air Capture Simulation.
Version 13 did not run.
This version includes checking the version of the parameter file and replacing it automatically if necessary.
The key update is to add code to evaluate performance if the Air Capture subsystem is cooled before flight.
The default temperature will be set at 100 Kelvin. That may be difficult to achieve in practice, but we need a value to work with.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
SPECIFIC_HEAT_AIR_J_KG_K = 1005.0 # Specific heat of air at constant pressure
SPECIFIC_HEAT_PIPE_J_KG_K = 450.0 # Specific heat of steel (J/kg*K)
# Molar masses (kg/mol)
MOLAR_MASS_DRY_AIR = 0.028965
MOLAR_MASS_WATER = 0.018015
# Mass fractions of dry air
MASS_FRACTION_NITROGEN = 0.755
MASS_FRACTION_OXYGEN = 0.232
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
VERSION = 14
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_vapor_pressure(temp_celsius):
"""Calculates saturation vapor pressure in Pa using the Arden Buck equation."""
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius)))
return es * 1000.0
def get_air_density(altitude_meters, relative_humidity_percent=0.0):
"""
Calculates moist air density using a simplified exponential model and humidity.
Relative humidity is a percentage from 0 to 100.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = (PRESS_0_PA / (287.058 * TEMP_0_K)) * (1 - e_pa / PRESS_0_PA * (1 - MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR))
return rho_moist_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K, relative_humidity_percent=0.0):
"""
Estimates pressure of collected gas using Ideal Gas Law, adjusted for molar mass.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
mole_fraction_water = e_pa / PRESS_0_PA
mole_fraction_dry_air = 1.0 - mole_fraction_water
avg_molar_mass = (mole_fraction_dry_air * MOLAR_MASS_DRY_AIR) + (mole_fraction_water * MOLAR_MASS_WATER)
R = 8.314
n = mass_kg / avg_molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
pipe_mass_kg,
payload_mass_kg,
rocket_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor,
relative_humidity_percent,
initial_pipe_temp_K,
version # Added new parameter
):
"""
Simulates a projectile launch with subsystem mass inputs and pre-cooling.
"""
initial_projectile_mass_kg = pipe_mass_kg + payload_mass_kg + rocket_mass_kg
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
gas_temperature_K = initial_pipe_temp_K
pipe_temperature_K = initial_pipe_temp_K
print(f"--- Gas Gun Air Harvesting Simulation V{version} ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Pipe Mass: {pipe_mass_kg} kg")
print(f"Payload Mass: {payload_mass_kg} kg")
print(f"Rocket Mass: {rocket_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Relative Humidity: {relative_humidity_percent:.1f}%")
print(f"Initial Pipe Temp: {initial_pipe_temp_K:.1f} K")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m, relative_humidity_percent)
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
# Simplified heat exchange model
total_mass_in_pipe = collected_air_mass + mass_collected_in_step
if total_mass_in_pipe > 0:
# Energy of existing air + energy of new air + kinetic energy converted to heat
# Energy of new air is at ambient temp (TEMP_0_K)
energy_in = (collected_air_mass * SPECIFIC_HEAT_AIR_J_KG_K * gas_temperature_K) + \
(mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K * TEMP_0_K) + \
(0.5 * mass_collected_in_step * current_velocity_mps**2)
# Energy of pipe
energy_pipe = pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K * pipe_temperature_K
# Distribute total energy across combined system (air + pipe)
total_energy_combined = energy_in + energy_pipe
total_thermal_mass = (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K) + (pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K)
# New equilibrium temperature
new_temp = total_energy_combined / total_thermal_mass
gas_temperature_K = new_temp
pipe_temperature_K = new_temp
internal_pressure = get_gas_pressure(total_mass_in_pipe, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
else:
internal_pressure = 0.0
if use_door_closing:
if internal_pressure >= dynamic_pressure and total_mass_in_pipe > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and total_mass_in_pipe > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Total Time: {current_time_s:.2f} s")
print(f"Initial Kinetic Energy: {0.5 * initial_projectile_mass_kg * initial_velocity_mps**2 / 1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
print("\n--- Mass Breakdown at Closure ---")
print(f"Total Projectile Mass: {initial_projectile_mass_kg:.2f} kg")
print(f"Magic Pipe Mass: {pipe_mass_kg:.2f} kg")
print(f"Payload Mass: {payload_mass_kg:.2f} kg")
print(f"Rocket Mass: {rocket_mass_kg:.2f} kg")
print(f"\nMass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if collected_air_mass > 0:
print("\n--- Payload Composition (by mass) ---")
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = get_air_density(0, relative_humidity_percent)
rho_water_vapor_0 = (e_pa / (461.5 * TEMP_0_K))
mass_of_water = collected_air_mass * (rho_water_vapor_0 / rho_moist_0)
mass_of_dry_air = collected_air_mass - mass_of_water
mass_of_nitrogen = mass_of_dry_air * MASS_FRACTION_NITROGEN
mass_of_oxygen = mass_of_dry_air * MASS_FRACTION_OXYGEN
percent_water = (mass_of_water / collected_air_mass) * 100.0
percent_nitrogen = (mass_of_nitrogen / collected_air_mass) * 100.0
percent_oxygen = (mass_of_oxygen / collected_air_mass) * 100.0
print(f"Water: {mass_of_water:.2f} kg ({percent_water:.2f}%)")
print(f"Nitrogen: {mass_of_nitrogen:.2f} kg ({percent_nitrogen:.2f}%)")
print(f"Oxygen: {mass_of_oxygen:.2f} kg ({percent_oxygen:.2f}%)")
other_mass = collected_air_mass - mass_of_water - mass_of_nitrogen - mass_of_oxygen
other_percent = (other_mass / collected_air_mass) * 100.0
print(f"Other: {other_mass:.2f} kg ({other_percent:.2f}%)")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
default_params = {
"projectile_diameter_m": 1.0,
"pipe_mass_kg": 1000.0,
"payload_mass_kg": 1000.0,
"rocket_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": True,
"outflow_rate_factor": 1.0e-5,
"relative_humidity_percent": 0.0,
"initial_pipe_temp_K": 100.0,
"version": VERSION
}
if os.path.exists(PARAMS_FILE):
try:
with open(PARAMS_FILE, 'r') as f:
params = json.load(f)
if params.get("version") != VERSION:
print(f"Parameter file from old version found. Deleting and creating new file for V{VERSION}.")
os.remove(PARAMS_FILE)
return default_params
return params
except (json.JSONDecodeError, KeyError):
print("Corrupted or invalid parameter file found. Deleting and creating a new one.")
if os.path.exists(PARAMS_FILE):
os.remove(PARAMS_FILE)
return default_params
else:
return default_params
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["pipe_mass_kg"] = get_user_input("Enter new magic pipe mass (kg)", params["pipe_mass_kg"], float)
params["payload_mass_kg"] = get_user_input("Enter new payload mass (kg)", params["payload_mass_kg"], float)
params["rocket_mass_kg"] = get_user_input("Enter new rocket mass (kg)", params["rocket_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
params["relative_humidity_percent"] = get_user_input("Enter new relative humidity (%)", params["relative_humidity_percent"], float)
params["initial_pipe_temp_K"] = get_user_input("Enter new initial pipe temperature (K)", params["initial_pipe_temp_K"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()(th)
Offline
Like button can go here
This post is for Version 15 of the Air Capture Simulation ...
Gemini explained that the physics model was incomplete in Version 14. The model has been extended to spread incoming energy over time as is appropriate for the situation we are studying.
import math
import json
import os
# --- Constants for Earth and Atmosphere ---
G = 6.67430e-11 # Gravitational constant (N * m^2 / kg^2)
M_EARTH = 5.972e24 # Mass of Earth (kg)
R_EARTH = 6.371e6 # Radius of Earth (m)
RHO_0 = 1.225 # Sea-level air density (kg/m^3)
H = 8500.0 # Atmospheric scale height (m)
CD = 0.2 # Drag coefficient for a streamlined body (assumed)
V_EXHAUST_ASSUMED = 1000.0 # Assumed exhaust velocity for outflow (m/s)
# --- Constants for Thermodynamic Calculations ---
TEMP_0_K = 288.0 # Sea-level standard temperature (Kelvin)
PRESS_0_PA = 101325.0 # Sea-level standard pressure (Pascals)
GAMMA = 1.4 # Heat capacity ratio for air
SPECIFIC_HEAT_AIR_J_KG_K = 1005.0 # Specific heat of air at constant pressure
SPECIFIC_HEAT_PIPE_J_KG_K = 450.0 # Specific heat of steel (J/kg*K)
# Molar masses (kg/mol)
MOLAR_MASS_DRY_AIR = 0.028965
MOLAR_MASS_WATER = 0.018015
# Mass fractions of dry air
MASS_FRACTION_NITROGEN = 0.755
MASS_FRACTION_OXYGEN = 0.232
# --- File for saving parameters ---
PARAMS_FILE = "simulation_params.json"
VERSION = 15
def get_gravity(altitude_meters):
"""Calculates the acceleration due to gravity at a given altitude."""
r = R_EARTH + altitude_meters
return G * M_EARTH / r**2
def get_vapor_pressure(temp_celsius):
"""Calculates saturation vapor pressure in Pa using the Arden Buck equation."""
es = 0.61121 * math.exp((18.678 - (temp_celsius / 234.5)) * (temp_celsius / (257.14 + temp_celsius)))
return es * 1000.0
def get_air_density(altitude_meters, relative_humidity_percent=0.0):
"""
Calculates moist air density using a simplified exponential model and humidity.
Relative humidity is a percentage from 0 to 100.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = (PRESS_0_PA / (287.058 * TEMP_0_K)) * (1 - e_pa / PRESS_0_PA * (1 - MOLAR_MASS_WATER / MOLAR_MASS_DRY_AIR))
return rho_moist_0 * math.exp(-altitude_meters / H)
def get_gas_pressure(mass_kg, volume_m3, temp_K, relative_humidity_percent=0.0):
"""
Estimates pressure of collected gas using Ideal Gas Law, adjusted for molar mass.
"""
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
mole_fraction_water = e_pa / PRESS_0_PA
mole_fraction_dry_air = 1.0 - mole_fraction_water
avg_molar_mass = (mole_fraction_dry_air * MOLAR_MASS_DRY_AIR) + (mole_fraction_water * MOLAR_MASS_WATER)
R = 8.314
n = mass_kg / avg_molar_mass
if volume_m3 <= 0:
return float('inf')
return (n * R * temp_K) / volume_m3
def simulate_air_harvesting_launch(
projectile_diameter_m,
pipe_mass_kg,
payload_mass_kg,
rocket_mass_kg,
initial_velocity_mps,
launch_altitude_meters_above_sea_level,
final_altitude_meters_above_sea_level,
pipe_volume_m3,
time_step_s,
use_door_closing,
outflow_rate_factor,
relative_humidity_percent,
initial_pipe_temp_K,
version
):
"""
Simulates a projectile launch with subsystem mass inputs and pre-cooling.
"""
initial_projectile_mass_kg = pipe_mass_kg + payload_mass_kg + rocket_mass_kg
if initial_velocity_mps <= 0 or projectile_diameter_m <= 0 or initial_projectile_mass_kg <= 0:
print("Error: All initial values must be positive.")
return
current_mass_kg = initial_projectile_mass_kg
current_velocity_mps = initial_velocity_mps
current_altitude_m = launch_altitude_meters_above_sea_level
current_time_s = 0.0
cross_sectional_area_sqm = math.pi * (projectile_diameter_m / 2)**2
total_work_done_on_air_J = 0.0
total_work_done_against_gravity_J = 0.0
doors_closed = False
pipe_is_full = False
pipe_is_empty = False
collected_air_mass = 0.0
gas_temperature_K = TEMP_0_K
pipe_temperature_K = initial_pipe_temp_K
print(f"--- Gas Gun Air Harvesting Simulation V{version} ---")
print(f"Projectile Diameter: {projectile_diameter_m} m")
print(f"Initial Mass: {initial_projectile_mass_kg} kg")
print(f"Pipe Mass: {pipe_mass_kg} kg")
print(f"Payload Mass: {payload_mass_kg} kg")
print(f"Rocket Mass: {rocket_mass_kg} kg")
print(f"Initial Velocity: {initial_velocity_mps} m/s ({initial_velocity_mps/1000:.2f} km/s)")
print(f"Target Altitude: {final_altitude_meters_above_sea_level / 1000:.0f} km")
print(f"Pipe Volume: {pipe_volume_m3} m^3")
print(f"Relative Humidity: {relative_humidity_percent:.1f}%")
print(f"Initial Pipe Temp: {initial_pipe_temp_K:.1f} K")
print(f"Use Door Closing: {use_door_closing}")
print("------------------------------------------")
while current_altitude_m < final_altitude_meters_above_sea_level:
if current_velocity_mps <= 0:
print(f"Projectile stalled at {current_altitude_m / 1000:.2f} km. Velocity is zero or negative.")
break
g = get_gravity(current_altitude_m)
rho = get_air_density(current_altitude_m, relative_humidity_percent)
force_gravity = current_mass_kg * g
net_force = -force_gravity
dynamic_pressure = 0.5 * rho * current_velocity_mps**2
if not doors_closed and not pipe_is_full:
mass_collected_in_step = rho * cross_sectional_area_sqm * current_velocity_mps * time_step_s
# --- CORRECTED THERMAL MODEL ---
# Step 1: Calculate the work done on the incoming air
work_air_step = 0.5 * mass_collected_in_step * current_velocity_mps**2
# Step 2: Determine the potential temperature of the new air before heat exchange
temp_new_air = TEMP_0_K + (work_air_step / (mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K))
# Step 3: Calculate the heat transfer between the new air and the pipe
# Simplified heat transfer model: assumes perfect mixing and heat exchange
total_energy_air = (collected_air_mass * SPECIFIC_HEAT_AIR_J_KG_K * gas_temperature_K) + \
(mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K * temp_new_air)
total_mass_in_pipe = collected_air_mass + mass_collected_in_step
if total_mass_in_pipe > 0:
# Equilibrium temperature between the new and old air in the pipe
temp_air_mixed = total_energy_air / (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K)
# Energy transfer from the mixed air to the cold pipe
heat_transferred_to_pipe = SPECIFIC_HEAT_PIPE_J_KG_K * pipe_mass_kg * (temp_air_mixed - pipe_temperature_K)
# New temperatures of the air and the pipe after heat exchange
final_gas_energy = total_energy_air - heat_transferred_to_pipe
final_pipe_energy = SPECIFIC_HEAT_PIPE_J_KG_K * pipe_mass_kg * pipe_temperature_K + heat_transferred_to_pipe
gas_temperature_K = final_gas_energy / (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K)
pipe_temperature_K = final_pipe_energy / (pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K)
internal_pressure = get_gas_pressure(total_mass_in_pipe, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
# --- END CORRECTED THERMAL MODEL ---
if use_door_closing:
if internal_pressure >= dynamic_pressure and total_mass_in_pipe > 0:
doors_closed = True
print(f"Doors closed at {current_altitude_m / 1000:.2f} km.")
print(f"Time at closure: {current_time_s:.2f} s.")
print(f"Mass collected at closure: {collected_air_mass:.2f} kg.")
print(f"Dynamic pressure at closure: {dynamic_pressure:.2f} Pa.")
print(f"Temperature at closure: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
else:
if internal_pressure >= dynamic_pressure and total_mass_in_pipe > 0:
pipe_is_full = True
print(f"Pipe filled at {current_altitude_m / 1000:.2f} km.")
print(f"Time at fill: {current_time_s:.2f} s.")
print(f"Mass collected at fill: {collected_air_mass:.2f} kg.")
print(f"Temperature at fill: {gas_temperature_K:.2f} K ({gas_temperature_K - 273.15:.2f} °C)")
if not doors_closed and not pipe_is_full:
force_air_collection = mass_collected_in_step / time_step_s * current_velocity_mps
net_force -= force_air_collection
work_air_step = force_air_collection * (current_velocity_mps * time_step_s)
total_work_done_on_air_J += work_air_step
current_mass_kg += mass_collected_in_step
collected_air_mass += mass_collected_in_step
elif use_door_closing and doors_closed:
force_drag = 0.5 * rho * current_velocity_mps**2 * CD * cross_sectional_area_sqm
net_force -= force_drag
elif not use_door_closing and pipe_is_full:
if not pipe_is_empty:
internal_pressure = get_gas_pressure(collected_air_mass, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
outflow_rate = outflow_rate_factor * (internal_pressure - dynamic_pressure)
outflow_rate = max(outflow_rate, 0)
if collected_air_mass - outflow_rate * time_step_s <= 0:
outflow_rate = collected_air_mass / time_step_s
pipe_is_empty = True
print(f"Pipe emptied at {current_altitude_m / 1000:.2f} km.")
force_retro = outflow_rate * V_EXHAUST_ASSUMED
net_force -= force_retro
current_mass_kg -= outflow_rate * time_step_s
collected_air_mass -= outflow_rate * time_step_s
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
elif not use_door_closing and pipe_is_empty:
dm_dt_in = rho * cross_sectional_area_sqm * current_velocity_mps
force_air_collection = dm_dt_in * current_velocity_mps
net_force -= force_air_collection
current_mass_kg = initial_projectile_mass_kg
work_gravity_step = force_gravity * (current_velocity_mps * time_step_s)
total_work_done_against_gravity_J += work_gravity_step
acceleration = net_force / current_mass_kg
current_velocity_mps += acceleration * time_step_s
current_altitude_m += current_velocity_mps * time_step_s
current_time_s += time_step_s
print("\n--- Simulation Complete ---")
print(f"Final Altitude: {current_altitude_m / 1000:.2f} km")
print(f"Final Velocity: {current_velocity_mps:.2f} m/s ({current_velocity_mps/1000:.2f} km/s)")
print(f"Total Time: {current_time_s:.2f} s")
print(f"Initial Kinetic Energy: {0.5 * initial_projectile_mass_kg * initial_velocity_mps**2 / 1e9:.2f} GigaJoules")
print(f"Total Work Done Against Gravity: {total_work_done_against_gravity_J/1e9:.2f} GigaJoules")
print(f"Total Work Done on Air Mass: {total_work_done_on_air_J/1e9:.2f} GigaJoules (before doors closed)")
print("\n--- Mass Breakdown at Closure ---")
print(f"Total Projectile Mass: {initial_projectile_mass_kg:.2f} kg")
print(f"Magic Pipe Mass: {pipe_mass_kg:.2f} kg")
print(f"Payload Mass: {payload_mass_kg:.2f} kg")
print(f"Rocket Mass: {rocket_mass_kg:.2f} kg")
print(f"\nMass of Collected Air: {collected_air_mass:.2f} kg ({collected_air_mass / 1000:.2f} metric tons)")
if collected_air_mass > 0:
print("\n--- Payload Composition (by mass) ---")
temp_celsius = TEMP_0_K - 273.15
es_pa = get_vapor_pressure(temp_celsius)
e_pa = es_pa * (relative_humidity_percent / 100.0)
rho_moist_0 = get_air_density(0, relative_humidity_percent)
rho_water_vapor_0 = (e_pa / (461.5 * TEMP_0_K))
mass_of_water = collected_air_mass * (rho_water_vapor_0 / rho_moist_0)
mass_of_dry_air = collected_air_mass - mass_of_water
mass_of_nitrogen = mass_of_dry_air * MASS_FRACTION_NITROGEN
mass_of_oxygen = mass_of_dry_air * MASS_FRACTION_OXYGEN
percent_water = (mass_of_water / collected_air_mass) * 100.0
percent_nitrogen = (mass_of_nitrogen / collected_air_mass) * 100.0
percent_oxygen = (mass_of_oxygen / collected_air_mass) * 100.0
print(f"Water: {mass_of_water:.2f} kg ({percent_water:.2f}%)")
print(f"Nitrogen: {mass_of_nitrogen:.2f} kg ({percent_nitrogen:.2f}%)")
print(f"Oxygen: {mass_of_oxygen:.2f} kg ({percent_oxygen:.2f}%)")
other_mass = collected_air_mass - mass_of_water - mass_of_nitrogen - mass_of_oxygen
other_percent = (other_mass / collected_air_mass) * 100.0
print(f"Other: {other_mass:.2f} kg ({other_percent:.2f}%)")
def load_parameters():
"""Loads parameters from a JSON file, or returns defaults if the file doesn't exist."""
default_params = {
"projectile_diameter_m": 1.0,
"pipe_mass_kg": 1000.0,
"payload_mass_kg": 1000.0,
"rocket_mass_kg": 1000.0,
"initial_velocity_mps": 11000.0,
"launch_altitude_meters_above_sea_level": 0.0,
"final_altitude_meters_above_sea_level": 100000.0,
"pipe_volume_m3": 15.0,
"time_step_s": 0.01,
"use_door_closing": True,
"outflow_rate_factor": 1.0e-5,
"relative_humidity_percent": 0.0,
"initial_pipe_temp_K": 100.0,
"version": VERSION
}
if os.path.exists(PARAMS_FILE):
try:
with open(PARAMS_FILE, 'r') as f:
params = json.load(f)
if params.get("version") != VERSION:
print(f"Parameter file from old version found. Deleting and creating new file for V{VERSION}.")
os.remove(PARAMS_FILE)
return default_params
return params
except (json.JSONDecodeError, KeyError):
print("Corrupted or invalid parameter file found. Deleting and creating a new one.")
if os.path.exists(PARAMS_FILE):
os.remove(PARAMS_FILE)
return default_params
else:
return default_params
def save_parameters(params):
"""Saves the given parameters to a JSON file."""
with open(PARAMS_FILE, 'w') as f:
json.dump(params, f, indent=4)
print(f"\nParameters saved to {PARAMS_FILE}")
def get_user_input(prompt, default_value, type_func):
"""Gets user input with a default option."""
while True:
user_input = input(f"{prompt} (default: {default_value}): ")
if user_input == "":
return default_value
try:
if type_func is bool:
return user_input.lower() == 'true'
return type_func(user_input)
except ValueError:
print(f"Invalid input. Please enter a valid {type_func.__name__} or 'True'/'False'.")
def main():
"""Main function to handle user interaction and run the simulation."""
params = load_parameters()
print("\n--- Current Simulation Parameters ---")
for key, value in params.items():
print(f" {key}: {value}")
change_params = input("\nDo you want to change any parameters? (y/n): ").lower()
if change_params == 'y':
params["projectile_diameter_m"] = get_user_input("Enter new projectile diameter (m)", params["projectile_diameter_m"], float)
params["pipe_mass_kg"] = get_user_input("Enter new magic pipe mass (kg)", params["pipe_mass_kg"], float)
params["payload_mass_kg"] = get_user_input("Enter new payload mass (kg)", params["payload_mass_kg"], float)
params["rocket_mass_kg"] = get_user_input("Enter new rocket mass (kg)", params["rocket_mass_kg"], float)
params["initial_velocity_mps"] = get_user_input("Enter new initial velocity (m/s)", params["initial_velocity_mps"], float)
params["final_altitude_meters_above_sea_level"] = get_user_input("Enter new target altitude (m)", params["final_altitude_meters_above_sea_level"], float)
params["pipe_volume_m3"] = get_user_input("Enter new magic pipe volume (m^3)", params["pipe_volume_m3"], float)
params["time_step_s"] = get_user_input("Enter new time step (s)", params["time_step_s"], float)
params["use_door_closing"] = get_user_input("Use door closing feature? (True/False)", params["use_door_closing"], bool)
params["outflow_rate_factor"] = get_user_input("Enter new outflow rate factor", params["outflow_rate_factor"], float)
params["relative_humidity_percent"] = get_user_input("Enter new relative humidity (%)", params["relative_humidity_percent"], float)
params["initial_pipe_temp_K"] = get_user_input("Enter new initial pipe temperature (K)", params["initial_pipe_temp_K"], float)
save_option = input("Save these new parameters as defaults for future runs? (y/n): ").lower()
if save_option == 'y':
save_parameters(params)
simulate_air_harvesting_launch(**params)
if __name__ == "__main__":
main()This version includes correction of the parameters file management, so there should be no further glitches on that front.
(th)
Offline
Like button can go here
Looks like some functions should be smaller.
And also, some lines could be put into a function.
Have u guys looked into object oriented programming.
For those are very object-y things we are dealing with.
We make class projectile,
then an instance of that class and then thats a projectile.
Something like that.
When we are having as many arguments as simulate_air_harvesting_launch has.
Makes me wonder.
We already think in objects.
And everything about a projectile for instance, should just be locked to that specific projectile.
Then we can just send that instance as an argument to the function.
One arg which can contain any number of member variables.
I think for instance, we should have class planet in that code.
And then make an instance of it, earth.
Then
M_EARTH = 5.972e24 # Mass of Earth (kg)
could become
earth.mass = 5.972e24
We could make it with something like this.
# Here we decide, what we do when we create new planets. Its called a constructor.
# In the constructor, self is a reference to the object being created each time.
class planet():
# Here we make a member function.
# Thats a function all the instances of a planet, have access to.
# When this function is being called, self always means, whatever object is calling each time.
def add_mass(self,new_mass):
self.mass=new_mass;
#Lets make two planets.
earth = planet();
mars = planet();
#Lets give them both some mass.
earth.set_mass(300);
mars.set_mass(500);
And now..
M_EARTH = 5.972e24 # Mass of Earth (kg)
could become something like:
earth.set_mass(5.972e24)
And after that we can just talk about earth.mass.
The bigger the code the bigger the benefit of thinking like this.
We could have hundreds of variables and functions belonging to class planet.
We already think about planets this way in the real world.
And most objects for that matter.
We are not just reusing code we are also reusing our own thinking process.
And if we make another system in the future we can just reuse class planet.
Class planet is not really this computer program as such.
Its just a concept.
Last edited by offtherock (Yesterday 17:57:59)
Offline
Like button can go here
For offtherock re #73
Thanks for taking a look at the current version of the Air Capture Simulation, and suggesting ways we could improve the logic.
This post is intended to hold the results of a test of Version 15. This version is not yet handling cooling of incoming air correctly, when we cool the collection container. I will post the transcript of the run, and then a diff file showing the updates to version 14.
python3 20250826AirCaptureSimulationV15.py
--- Current Simulation Parameters ---
projectile_diameter_m: 1.0
pipe_mass_kg: 1000.0
payload_mass_kg: 1000.0
rocket_mass_kg: 1000.0
initial_velocity_mps: 11000.0
launch_altitude_meters_above_sea_level: 0.0
final_altitude_meters_above_sea_level: 100000.0
pipe_volume_m3: 15.0
time_step_s: 0.01
use_door_closing: True
outflow_rate_factor: 1e-05
relative_humidity_percent: 0.0
initial_pipe_temp_K: 100.0
version: 15Do you want to change any parameters? (y/n): n
--- Gas Gun Air Harvesting Simulation V15 ---
Projectile Diameter: 1.0 m
Initial Mass: 3000.0 kg
Pipe Mass: 1000.0 kg
Payload Mass: 1000.0 kg
Rocket Mass: 1000.0 kg
Initial Velocity: 11000.0 m/s (11.00 km/s)
Target Altitude: 100 km
Pipe Volume: 15.0 m^3
Relative Humidity: 0.0%
Initial Pipe Temp: 100.0 K
Use Door Closing: True
------------------------------------------
Doors closed at 0.11 km.
Time at closure: 0.01 s.
Mass collected at closure: 105.89 kg.
Dynamic pressure at closure: 68319558.87 Pa.
Temperature at closure: 214939.09 K (214665.94 °C)--- Simulation Complete ---
Final Altitude: 100.00 km
Final Velocity: 8087.26 m/s (8.09 km/s)
Total Time: 12.02 s
Initial Kinetic Energy: 181.50 GigaJoules
Total Work Done Against Gravity: 3.00 GigaJoules
Total Work Done on Air Mass: 12.81 GigaJoules (before doors closed)--- Mass Breakdown at Closure ---
Total Projectile Mass: 3000.00 kg
Magic Pipe Mass: 1000.00 kg
Payload Mass: 1000.00 kg
Rocket Mass: 1000.00 kgMass of Collected Air: 105.89 kg (0.11 metric tons)
--- Payload Composition (by mass) ---
Water: 0.00 kg (0.00%)
Nitrogen: 79.94 kg (75.50%)
Oxygen: 24.57 kg (23.20%)
Other: 1.38 kg (1.30%)
And here is a list of differences between version 15 on the left and 14 on the right:
diff 20250826AirCaptureSimulationV15.py archive/20250826AirCaptureSimulationV14.py
1c1
< # 20250826AirCaptureSimulationV15.py
---
> # 20250826AirCaptureSimulationV14.py
32c32
< VERSION = 15
---
> VERSION = 14
85c85
< version
---
> version # Added new parameter
109c109
< gas_temperature_K = TEMP_0_K
---
> gas_temperature_K = initial_pipe_temp_K
142,153c142
< # --- CORRECTED THERMAL MODEL ---
< # Step 1: Calculate the work done on the incoming air
< work_air_step = 0.5 * mass_collected_in_step * current_velocity_mps**2
<
< # Step 2: Determine the potential temperature of the new air before heat exchange
< temp_new_air = TEMP_0_K + (work_air_step / (mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K))
<
< # Step 3: Calculate the heat transfer between the new air and the pipe
< # Simplified heat transfer model: assumes perfect mixing and heat exchange
< total_energy_air = (collected_air_mass * SPECIFIC_HEAT_AIR_J_KG_K * gas_temperature_K) + \
< (mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K * temp_new_air)
<
---
> # Simplified heat exchange model
157,158c146,150
< # Equilibrium temperature between the new and old air in the pipe
< temp_air_mixed = total_energy_air / (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K)
---
> # Energy of existing air + energy of new air + kinetic energy converted to heat
> # Energy of new air is at ambient temp (TEMP_0_K)
> energy_in = (collected_air_mass * SPECIFIC_HEAT_AIR_J_KG_K * gas_temperature_K) + \
> (mass_collected_in_step * SPECIFIC_HEAT_AIR_J_KG_K * TEMP_0_K) + \
> (0.5 * mass_collected_in_step * current_velocity_mps**2)
160,161c152,153
< # Energy transfer from the mixed air to the cold pipe
< heat_transferred_to_pipe = SPECIFIC_HEAT_PIPE_J_KG_K * pipe_mass_kg * (temp_air_mixed - pipe_temperature_K)
---
> # Energy of pipe
> energy_pipe = pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K * pipe_temperature_K
163,165c155,157
< # New temperatures of the air and the pipe after heat exchange
< final_gas_energy = total_energy_air - heat_transferred_to_pipe
< final_pipe_energy = SPECIFIC_HEAT_PIPE_J_KG_K * pipe_mass_kg * pipe_temperature_K + heat_transferred_to_pipe
---
> # Distribute total energy across combined system (air + pipe)
> total_energy_combined = energy_in + energy_pipe
> total_thermal_mass = (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K) + (pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K)
167,168c159,160
< gas_temperature_K = final_gas_energy / (total_mass_in_pipe * SPECIFIC_HEAT_AIR_J_KG_K)
< pipe_temperature_K = final_pipe_energy / (pipe_mass_kg * SPECIFIC_HEAT_PIPE_J_KG_K)
---
> # New equilibrium temperature
> new_temp = total_energy_combined / total_thermal_mass
170,172c162,167
< internal_pressure = get_gas_pressure(total_mass_in_pipe, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
<
< # --- END CORRECTED THERMAL MODEL ---
---
> gas_temperature_K = new_temp
> pipe_temperature_K = new_temp
>
> internal_pressure = get_gas_pressure(total_mass_in_pipe, pipe_volume_m3, gas_temperature_K, relative_humidity_percent)
> else:
> internal_pressure = 0.0
357c352
< # lines in V15 357
---
> # lines in v14 352
(th)
Offline
Like button can go here
ur project reminds me of this nasa program they were
working on in the 90s.
it was super cool.
he heated the air above this tiny metal. and that created this upwards
momentum, which propelled the thing upwards.
and he could shoot the laster from elsewhere even.
off course, the thing that killed it, is that it doesnt work outside
the atmosphere.
still cool though, and interesting.
https://www.youtube.com/watch?v=LAdj6vpYppA
when an object is in the air,
the molecules are hitting it from all sides.
if the molecules above it, increase in temperature for instance.
two things will start to happen.
the molecules will be traveling faster, thus, hitting the object faster.
which is bad for they will tend to be on the way downwards.. pushing it down.
but also, the density of the air/molecules. should decrease.
meaning fewer molecules start hitting it downwards.
which is good.
i have no idea nor ability to evaluate the magnitude of each effect.
if we could magically just, decrease the density above an object.
shouldnt that lead to an upwards force on it.
i would think so.
because i would think the decrease in density effect would overshadow
the increase in speed effect.
i would think if a tunnel of vacuum, would, magically form over an object.
it would just, shoot right up.
cant we somehow use those lasers,
to change that density somehow.
to lead to motion.
in any direction.
speaking about lasers,
they seem to have been advancing in leaps and bounds lately.
https://www.youtube.com/watch?v=UBVlL0FNbSE
and apparently we can also create motion with sound.
https://www.youtube.com/watch?v=uENITui5_jU
https://www.bbc.com/news/science-environment-34647921
its like a riddle.
a riddle i think i will never solve.
but maybe somebody somewhere, will.
Last edited by offtherock (Today 04:18:08)
Offline
Like button can go here