Aerospace MCP

""" Orbital mechanics and spacecraft trajectory analysis tools for aerospace MCP. Provides orbital mechanics calculations with manual implementations and optional advanced library integration for poliastro/astropy/spiceypy when available. """ import math from dataclasses import dataclass from datetime import UTC, datetime from typing import Any # Constants MU_EARTH = 3.986004418e14 # m³/s² - Earth's gravitational parameter R_EARTH = 6378137.0 # m - Earth's equatorial radius (WGS84) J2_EARTH = 1.08262668e-3 # Earth's J2 perturbation coefficient OMEGA_EARTH = 7.2921159e-5 # rad/s - Earth's rotation rate @dataclass class OrbitElements: """Classical orbital elements.""" semi_major_axis_m: float # Semi-major axis (m) eccentricity: float # Eccentricity (dimensionless) inclination_deg: float # Inclination (degrees) raan_deg: float # Right ascension of ascending node (degrees) arg_periapsis_deg: float # Argument of periapsis (degrees) true_anomaly_deg: float # True anomaly (degrees) epoch_utc: str # Epoch in UTC ISO format @dataclass class StateVector: """Position and velocity state vector.""" position_m: list[float] # Position vector [x, y, z] in meters velocity_ms: list[float] # Velocity vector [vx, vy, vz] in m/s epoch_utc: str # Epoch in UTC ISO format frame: str = "J2000" # Reference frame @dataclass class OrbitProperties: """Computed orbital properties.""" period_s: float # Orbital period (seconds) apoapsis_m: float # Apoapsis altitude above Earth surface (m) periapsis_m: float # Periapsis altitude above Earth surface (m) energy_j_kg: float # Specific orbital energy (J/kg) angular_momentum_m2s: float # Specific angular momentum magnitude (m²/s) @dataclass class GroundTrack: """Ground track point.""" latitude_deg: float longitude_deg: float altitude_m: float time_utc: str @dataclass class Maneuver: """Orbital maneuver definition.""" delta_v_ms: list[float] # Delta-V vector [x, y, z] in m/s time_utc: str # Maneuver execution time description: str = "" # Optional description def deg_to_rad(deg: float) -> float: """Convert degrees to radians.""" return deg * math.pi / 180.0 def rad_to_deg(rad: float) -> float: """Convert radians to degrees.""" return rad * 180.0 / math.pi def vector_magnitude(vec: list[float]) -> float: """Calculate vector magnitude.""" return math.sqrt(sum(x**2 for x in vec)) def vector_dot(a: list[float], b: list[float]) -> float: """Calculate dot product of two vectors.""" return sum(a[i] * b[i] for i in range(len(a))) def vector_cross(a: list[float], b: list[float]) -> list[float]: """Calculate cross product of two 3D vectors.""" return [ a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0], ] def vector_normalize(vec: list[float]) -> list[float]: """Normalize a vector.""" mag = vector_magnitude(vec) if mag == 0: return [0.0, 0.0, 0.0] return [x / mag for x in vec] def kepler_equation_solver( mean_anomaly_rad: float, eccentricity: float, tolerance: float = 1e-8, max_iterations: int = 50, ) -> float: """ Solve Kepler's equation for eccentric anomaly using Newton-Raphson method. Args: mean_anomaly_rad: Mean anomaly in radians eccentricity: Orbital eccentricity tolerance: Convergence tolerance max_iterations: Maximum iterations Returns: Eccentric anomaly in radians """ # Initial guess E = mean_anomaly_rad if eccentricity < 0.8 else math.pi for _ in range(max_iterations): f = E - eccentricity * math.sin(E) - mean_anomaly_rad f_prime = 1 - eccentricity * math.cos(E) if abs(f_prime) < 1e-12: break E_new = E - f / f_prime if abs(E_new - E) < tolerance: return E_new E = E_new return E def elements_to_state_vector(elements: OrbitElements) -> StateVector: """ Convert orbital elements to state vector using manual calculations. Args: elements: Orbital elements Returns: State vector in J2000 frame """ # Convert angles to radians i = deg_to_rad(elements.inclination_deg) raan = deg_to_rad(elements.raan_deg) arg_pe = deg_to_rad(elements.arg_periapsis_deg) nu = deg_to_rad(elements.true_anomaly_deg) # Semi-major axis and eccentricity a = elements.semi_major_axis_m e = elements.eccentricity # Calculate distance and flight path angle r = a * (1 - e**2) / (1 + e * math.cos(nu)) # Position in perifocal coordinates r_peri = [r * math.cos(nu), r * math.sin(nu), 0.0] # Velocity in perifocal coordinates p = a * (1 - e**2) v_peri = [ -math.sqrt(MU_EARTH / p) * math.sin(nu), math.sqrt(MU_EARTH / p) * (e + math.cos(nu)), 0.0, ] # Rotation matrices cos_raan, sin_raan = math.cos(raan), math.sin(raan) cos_i, sin_i = math.cos(i), math.sin(i) cos_arg, sin_arg = math.cos(arg_pe), math.sin(arg_pe) # Rotation matrix from perifocal to J2000 R11 = cos_raan * cos_arg - sin_raan * sin_arg * cos_i R12 = -cos_raan * sin_arg - sin_raan * cos_arg * cos_i R13 = sin_raan * sin_i R21 = sin_raan * cos_arg + cos_raan * sin_arg * cos_i R22 = -sin_raan * sin_arg + cos_raan * cos_arg * cos_i R23 = -cos_raan * sin_i R31 = sin_arg * sin_i R32 = cos_arg * sin_i R33 = cos_i # Transform to J2000 frame r_j2000 = [ R11 * r_peri[0] + R12 * r_peri[1] + R13 * r_peri[2], R21 * r_peri[0] + R22 * r_peri[1] + R23 * r_peri[2], R31 * r_peri[0] + R32 * r_peri[1] + R33 * r_peri[2], ] v_j2000 = [ R11 * v_peri[0] + R12 * v_peri[1] + R13 * v_peri[2], R21 * v_peri[0] + R22 * v_peri[1] + R23 * v_peri[2], R31 * v_peri[0] + R32 * v_peri[1] + R33 * v_peri[2], ] return StateVector( position_m=r_j2000, velocity_ms=v_j2000, epoch_utc=elements.epoch_utc, frame="J2000", ) def state_vector_to_elements(state: StateVector) -> OrbitElements: """ Convert state vector to orbital elements using manual calculations. Args: state: State vector in J2000 frame Returns: Classical orbital elements """ r_vec = state.position_m v_vec = state.velocity_ms # Position and velocity magnitudes r = vector_magnitude(r_vec) v = vector_magnitude(v_vec) # Specific angular momentum h_vec = vector_cross(r_vec, v_vec) h = vector_magnitude(h_vec) # Semi-major axis energy = v**2 / 2 - MU_EARTH / r a = -MU_EARTH / (2 * energy) # Eccentricity vector v_cross_h = vector_cross(v_vec, h_vec) e_vec = [v_cross_h[i] / MU_EARTH - r_vec[i] / r for i in range(3)] e = vector_magnitude(e_vec) # Inclination i = math.acos(h_vec[2] / h) # Node vector k_vec = [0, 0, 1] n_vec = vector_cross(k_vec, h_vec) n = vector_magnitude(n_vec) # RAAN if n > 1e-10: raan = math.acos(n_vec[0] / n) if n_vec[1] < 0: raan = 2 * math.pi - raan else: raan = 0.0 # Argument of periapsis if n > 1e-10 and e > 1e-10: cos_arg_pe = vector_dot(n_vec, e_vec) / (n * e) cos_arg_pe = max(-1, min(1, cos_arg_pe)) # Clamp to [-1, 1] arg_pe = math.acos(cos_arg_pe) if e_vec[2] < 0: arg_pe = 2 * math.pi - arg_pe else: arg_pe = 0.0 # True anomaly if e > 1e-10: cos_nu = vector_dot(e_vec, r_vec) / (e * r) cos_nu = max(-1, min(1, cos_nu)) # Clamp to [-1, 1] nu = math.acos(cos_nu) if vector_dot(r_vec, v_vec) < 0: nu = 2 * math.pi - nu else: # For circular orbits, use longitude of ascending node if n > 1e-10: cos_nu = vector_dot(n_vec, r_vec) / (n * r) cos_nu = max(-1, min(1, cos_nu)) nu = math.acos(cos_nu) if r_vec[2] < 0: nu = 2 * math.pi - nu else: nu = math.atan2(r_vec[1], r_vec[0]) if nu < 0: nu += 2 * math.pi return OrbitElements( semi_major_axis_m=a, eccentricity=e, inclination_deg=rad_to_deg(i), raan_deg=rad_to_deg(raan), arg_periapsis_deg=rad_to_deg(arg_pe), true_anomaly_deg=rad_to_deg(nu), epoch_utc=state.epoch_utc, ) def calculate_orbit_properties(elements: OrbitElements) -> OrbitProperties: """ Calculate orbital properties from elements. Args: elements: Orbital elements Returns: Orbital properties """ a = elements.semi_major_axis_m e = elements.eccentricity # Orbital period period = 2 * math.pi * math.sqrt(a**3 / MU_EARTH) # Apoapsis and periapsis altitudes r_ap = a * (1 + e) - R_EARTH r_pe = a * (1 - e) - R_EARTH # Specific orbital energy energy = -MU_EARTH / (2 * a) # Specific angular momentum h = math.sqrt(MU_EARTH * a * (1 - e**2)) return OrbitProperties( period_s=period, apoapsis_m=r_ap, periapsis_m=r_pe, energy_j_kg=energy, angular_momentum_m2s=h, ) def propagate_orbit_j2( initial_state: StateVector, time_span_s: float, time_step_s: float = 60.0 ) -> list[StateVector]: """ Propagate orbit with J2 perturbations using numerical integration. Args: initial_state: Initial state vector time_span_s: Propagation time span (seconds) time_step_s: Integration time step (seconds) Returns: List of state vectors over time """ def acceleration_j2(r_vec: list[float]) -> list[float]: """Calculate acceleration including J2 perturbations.""" r = vector_magnitude(r_vec) # Central body acceleration a_central = [-MU_EARTH * r_vec[i] / r**3 for i in range(3)] # J2 perturbation factor = 1.5 * J2_EARTH * MU_EARTH * R_EARTH**2 / r**5 z2_r2 = (r_vec[2] / r) ** 2 a_j2 = [ factor * r_vec[0] * (1 - 5 * z2_r2), factor * r_vec[1] * (1 - 5 * z2_r2), factor * r_vec[2] * (3 - 5 * z2_r2), ] return [a_central[i] + a_j2[i] for i in range(3)] # Initialize states = [initial_state] r = initial_state.position_m.copy() v = initial_state.velocity_ms.copy() # Parse initial epoch try: epoch = datetime.fromisoformat(initial_state.epoch_utc.replace("Z", "+00:00")) except (ValueError, TypeError): epoch = datetime.now(UTC) # Numerical integration (RK4) t = 0.0 while t < time_span_s: dt = min(time_step_s, time_span_s - t) # RK4 integration k1_r = v k1_v = acceleration_j2(r) r2 = [r[i] + 0.5 * dt * k1_r[i] for i in range(3)] v2 = [v[i] + 0.5 * dt * k1_v[i] for i in range(3)] k2_r = v2 k2_v = acceleration_j2(r2) r3 = [r[i] + 0.5 * dt * k2_r[i] for i in range(3)] v3 = [v[i] + 0.5 * dt * k2_v[i] for i in range(3)] k3_r = v3 k3_v = acceleration_j2(r3) r4 = [r[i] + dt * k3_r[i] for i in range(3)] v4 = [v[i] + dt * k3_v[i] for i in range(3)] k4_r = v4 k4_v = acceleration_j2(r4) # Update state r = [ r[i] + dt / 6 * (k1_r[i] + 2 * k2_r[i] + 2 * k3_r[i] + k4_r[i]) for i in range(3) ] v = [ v[i] + dt / 6 * (k1_v[i] + 2 * k2_v[i] + 2 * k3_v[i] + k4_v[i]) for i in range(3) ] t += dt # Create new state new_epoch = epoch.replace(microsecond=0) + type(epoch - epoch)(seconds=int(t)) states.append( StateVector( position_m=r.copy(), velocity_ms=v.copy(), epoch_utc=new_epoch.isoformat(), frame=initial_state.frame, ) ) return states def calculate_ground_track( orbit_states: list[StateVector], time_step_s: float = 60.0 ) -> list[GroundTrack]: """ Calculate ground track from orbit state vectors. Args: orbit_states: List of state vectors time_step_s: Time step between states (seconds) Returns: List of ground track points """ ground_track = [] for i, state in enumerate(orbit_states): r_vec = state.position_m # Convert to latitude/longitude r = vector_magnitude(r_vec) lat = math.asin(r_vec[2] / r) # Account for Earth rotation try: datetime.fromisoformat(state.epoch_utc.replace("Z", "+00:00")) t_since_epoch = i * time_step_s lon = math.atan2(r_vec[1], r_vec[0]) - OMEGA_EARTH * t_since_epoch except Exception: lon = math.atan2(r_vec[1], r_vec[0]) # Normalize longitude to [-180, 180] degrees while lon > math.pi: lon -= 2 * math.pi while lon < -math.pi: lon += 2 * math.pi altitude = r - R_EARTH ground_track.append( GroundTrack( latitude_deg=rad_to_deg(lat), longitude_deg=rad_to_deg(lon), altitude_m=altitude, time_utc=state.epoch_utc, ) ) return ground_track def hohmann_transfer(r1_m: float, r2_m: float) -> dict[str, float]: """ Calculate Hohmann transfer orbit parameters. Args: r1_m: Initial circular orbit radius (m) r2_m: Final circular orbit radius (m) Returns: Transfer parameters including delta-V requirements """ # Transfer orbit semi-major axis a_transfer = (r1_m + r2_m) / 2 # Velocities v1_circular = math.sqrt(MU_EARTH / r1_m) v2_circular = math.sqrt(MU_EARTH / r2_m) v1_transfer = math.sqrt(MU_EARTH * (2 / r1_m - 1 / a_transfer)) v2_transfer = math.sqrt(MU_EARTH * (2 / r2_m - 1 / a_transfer)) # Delta-V requirements (signed values) dv1 = v1_transfer - v1_circular # Positive for prograde, negative for retrograde dv2 = v2_circular - v2_transfer # Positive for prograde, negative for retrograde dv_total = abs(dv1) + abs(dv2) # Total magnitude # Transfer time transfer_time = math.pi * math.sqrt(a_transfer**3 / MU_EARTH) return { "delta_v1_ms": dv1, "delta_v2_ms": dv2, "total_delta_v_ms": dv_total, "transfer_time_s": transfer_time, "transfer_time_h": transfer_time / 3600, "transfer_semi_major_axis_m": a_transfer, } def lambert_solver_simple( r1_vec: list[float], r2_vec: list[float], time_flight_s: float, mu: float = MU_EARTH ) -> dict[str, Any]: """ Simple Lambert problem solver for two-body trajectory. Args: r1_vec: Initial position vector (m) r2_vec: Final position vector (m) time_flight_s: Time of flight (seconds) mu: Gravitational parameter (m³/s²) Returns: Dictionary with initial and final velocity vectors """ r1 = vector_magnitude(r1_vec) r2 = vector_magnitude(r2_vec) # Chord length c = vector_magnitude([r2_vec[i] - r1_vec[i] for i in range(3)]) # Semi-perimeter s = (r1 + r2 + c) / 2 # Minimum energy ellipse semi-major axis a_min = s / 2 # Check if transfer time is feasible t_min = math.pi * math.sqrt(a_min**3 / mu) if time_flight_s < t_min: return { "feasible": False, "reason": f"Transfer time {time_flight_s:.1f}s less than minimum {t_min:.1f}s", "v1_ms": [0, 0, 0], "v2_ms": [0, 0, 0], } # Simplified solution assuming elliptical transfer # This is an approximation - full Lambert solver requires iterative methods # Transfer angle (simplified) cos_dnu = vector_dot(r1_vec, r2_vec) / (r1 * r2) cos_dnu = max(-1, min(1, cos_dnu)) # Clamp dnu = math.acos(cos_dnu) # Approximate semi-major axis for given flight time # Using Kepler's 3rd law and approximation n_approx = 2 * math.pi / time_flight_s # Approximate mean motion a_approx = (mu / n_approx**2) ** (1 / 3) # Approximate velocities (simplified circular approximation) v1_mag = math.sqrt(mu / r1) v2_mag = math.sqrt(mu / r2) # Direction vectors (perpendicular to radius for circular approximation) r1_unit = vector_normalize(r1_vec) r2_unit = vector_normalize(r2_vec) # Simplified velocity directions (tangential) h_vec = vector_cross(r1_vec, r2_vec) h_unit = vector_normalize(h_vec) v1_dir = vector_normalize(vector_cross(h_unit, r1_unit)) v2_dir = vector_normalize(vector_cross(h_unit, r2_unit)) v1_vec = [v1_mag * v1_dir[i] for i in range(3)] v2_vec = [v2_mag * v2_dir[i] for i in range(3)] return { "feasible": True, "v1_ms": v1_vec, "v2_ms": v2_vec, "v1_vec_ms": v1_vec, # Add expected key "initial_velocity_ms": v1_vec, # Add expected key "transfer_angle_deg": rad_to_deg(dnu), "estimated_semi_major_axis_m": a_approx, } def orbital_rendezvous_planning( chaser_elements: OrbitElements, target_elements: OrbitElements ) -> dict[str, Any]: """ Plan orbital rendezvous maneuvers between two spacecraft. Args: chaser_elements: Chaser spacecraft orbital elements target_elements: Target spacecraft orbital elements Returns: Rendezvous plan with phasing and approach maneuvers """ # Convert to state vectors for analysis chaser_state = elements_to_state_vector(chaser_elements) target_state = elements_to_state_vector(target_elements) # Calculate relative position rel_pos = [ target_state.position_m[i] - chaser_state.position_m[i] for i in range(3) ] rel_distance = vector_magnitude(rel_pos) # Orbital properties chaser_props = calculate_orbit_properties(chaser_elements) target_props = calculate_orbit_properties(target_elements) # Phase angle between spacecraft chaser_r = vector_magnitude(chaser_state.position_m) target_r = vector_magnitude(target_state.position_m) cos_phase = vector_dot(chaser_state.position_m, target_state.position_m) / ( chaser_r * target_r ) cos_phase = max(-1, min(1, cos_phase)) phase_angle_rad = math.acos(cos_phase) # Estimate phasing time if abs(chaser_props.period_s - target_props.period_s) > 1: synodic_period = abs( chaser_props.period_s * target_props.period_s / (chaser_props.period_s - target_props.period_s) ) phasing_time_s = phase_angle_rad / (2 * math.pi) * synodic_period else: phasing_time_s = float("inf") # Coplanar, similar orbits # Delta-V estimates for circularization if needed chaser_circ_dv = 0.0 if chaser_elements.eccentricity > 0.01: # Circularization delta-V (rough estimate) v_ap = math.sqrt( MU_EARTH * ( 2 / ( chaser_elements.semi_major_axis_m * (1 + chaser_elements.eccentricity) ) - 1 / chaser_elements.semi_major_axis_m ) ) v_circ = math.sqrt( MU_EARTH / (chaser_elements.semi_major_axis_m * (1 + chaser_elements.eccentricity)) ) chaser_circ_dv = abs(v_circ - v_ap) # Create sample maneuvers maneuvers = [] # Add a simple phasing maneuver if chaser_circ_dv > 0: maneuvers.append( { "delta_v_ms": chaser_circ_dv, "type": "circularization", "description": "Circularize chaser orbit", } ) # Add altitude matching maneuver if needed altitude_diff = abs(chaser_props.apoapsis_m - target_props.apoapsis_m) if altitude_diff > 1000: # More than 1 km difference altitude_dv = 50.0 # Rough estimate for altitude adjustment maneuvers.append( { "delta_v_ms": altitude_dv, "type": "altitude_adjustment", "description": "Match target altitude", } ) # Calculate total delta-V total_dv = sum(m["delta_v_ms"] for m in maneuvers) return { "relative_distance_m": rel_distance, "relative_distance_km": rel_distance / 1000, "phase_angle_deg": rad_to_deg(phase_angle_rad), "phasing_time_s": phasing_time_s, "phasing_time_h": ( phasing_time_s / 3600 if phasing_time_s != float("inf") else float("inf") ), "time_to_rendezvous_s": phasing_time_s, # Add expected key "chaser_period_s": chaser_props.period_s, "target_period_s": target_props.period_s, "period_difference_s": abs(chaser_props.period_s - target_props.period_s), "altitude_difference_m": abs(chaser_props.apoapsis_m - target_props.apoapsis_m), "estimated_circularization_dv_ms": chaser_circ_dv, "total_delta_v_ms": total_dv, # Add expected key "maneuvers": maneuvers, # Add expected key "feasibility": ( "Good" if rel_distance < 100000 and abs(chaser_elements.inclination_deg - target_elements.inclination_deg) < 5.0 else "Challenging" ), } # Update availability def porkchop_plot_analysis( departure_body: str = "Earth", arrival_body: str = "Mars", departure_dates: list[str] = None, arrival_dates: list[str] = None, min_tof_days: int = 100, max_tof_days: int = 400, ) -> dict[str, Any]: """ Generate porkchop plot analysis for interplanetary transfers. Args: departure_body: Departure celestial body (default: Earth) arrival_body: Arrival celestial body (default: Mars) departure_dates: List of departure dates (ISO format) arrival_dates: List of arrival dates (ISO format) min_tof_days: Minimum time of flight (days) max_tof_days: Maximum time of flight (days) Returns: Dictionary containing transfer analysis grid """ # Default date ranges if not provided if departure_dates is None: # Earth-Mars synodic period is ~26 months, sample over 2 years departure_dates = [ "2025-01-01T00:00:00", "2025-03-01T00:00:00", "2025-05-01T00:00:00", "2025-07-01T00:00:00", "2025-09-01T00:00:00", "2025-11-01T00:00:00", "2026-01-01T00:00:00", "2026-03-01T00:00:00", "2026-05-01T00:00:00", "2026-07-01T00:00:00", "2026-09-01T00:00:00", "2026-11-01T00:00:00", ] if arrival_dates is None: # Generate arrival dates based on TOF constraints arrival_dates = [ "2025-06-01T00:00:00", "2025-08-01T00:00:00", "2025-10-01T00:00:00", "2025-12-01T00:00:00", "2026-02-01T00:00:00", "2026-04-01T00:00:00", "2026-06-01T00:00:00", "2026-08-01T00:00:00", "2026-10-01T00:00:00", "2026-12-01T00:00:00", "2027-02-01T00:00:00", "2027-04-01T00:00:00", ] # Simplified planetary positions (circular orbits for demonstration) # In production, would use SPICE kernels or JPL ephemeris earth_orbit_radius = 1.0 # AU mars_orbit_radius = 1.524 # AU earth_period_days = 365.25 mars_period_days = 686.98 def get_planet_position(body: str, date_iso: str) -> list[float]: """Get simplified planet position at given date.""" # Parse date (simplified) import datetime try: dt = datetime.datetime.fromisoformat(date_iso.replace("Z", "+00:00")) except Exception: dt = datetime.datetime.fromisoformat(date_iso) # Days since J2000 j2000 = datetime.datetime(2000, 1, 1, 12, 0, 0) days_since_j2000 = (dt - j2000).total_seconds() / 86400 if body.lower() == "earth": # Earth mean anomaly mean_anomaly = 2 * math.pi * days_since_j2000 / earth_period_days x = earth_orbit_radius * math.cos(mean_anomaly) y = earth_orbit_radius * math.sin(mean_anomaly) return [x * 1.496e11, y * 1.496e11, 0.0] # Convert AU to meters elif body.lower() == "mars": # Mars mean anomaly mean_anomaly = 2 * math.pi * days_since_j2000 / mars_period_days x = mars_orbit_radius * math.cos(mean_anomaly) y = mars_orbit_radius * math.sin(mean_anomaly) return [x * 1.496e11, y * 1.496e11, 0.0] # Convert AU to meters else: # Default to Earth position mean_anomaly = 2 * math.pi * days_since_j2000 / earth_period_days x = earth_orbit_radius * math.cos(mean_anomaly) y = earth_orbit_radius * math.sin(mean_anomaly) return [x * 1.496e11, y * 1.496e11, 0.0] # Generate porkchop data grid porkchop_data = [] for dep_date in departure_dates: for arr_date in arrival_dates: # Calculate time of flight import datetime try: dep_dt = datetime.datetime.fromisoformat( dep_date.replace("Z", "+00:00") ) except Exception: dep_dt = datetime.datetime.fromisoformat(dep_date) try: arr_dt = datetime.datetime.fromisoformat( arr_date.replace("Z", "+00:00") ) except Exception: arr_dt = datetime.datetime.fromisoformat(arr_date) tof_days = (arr_dt - dep_dt).total_seconds() / 86400 # Filter by TOF constraints if tof_days < min_tof_days or tof_days > max_tof_days: continue # Get planetary positions r1 = get_planet_position(departure_body, dep_date) r2 = get_planet_position(arrival_body, arr_date) # Calculate transfer using simplified Lambert solver try: mu_sun = 1.32712440018e20 # Sun's gravitational parameter (m³/s²) tof_seconds = tof_days * 86400 # Use simplified Lambert solver lambert_result = lambert_solver_simple(r1, r2, tof_seconds, mu_sun) # Calculate characteristic energy (C3) if "v1_vec_ms" in lambert_result: v1 = lambert_result["v1_vec_ms"] elif "initial_velocity_ms" in lambert_result: v1 = lambert_result["initial_velocity_ms"] else: v1 = [0, 0, 0] # Fallback # Earth's orbital velocity at 1 AU earth_v_orbit = math.sqrt(mu_sun / (1.496e11)) # m/s # Excess velocity magnitude v_inf_vec = [ v1[i] - earth_v_orbit if i == 1 else v1[i] for i in range(3) ] v_inf_mag = math.sqrt(sum(v**2 for v in v_inf_vec)) # Characteristic energy C3 (km²/s²) c3 = (v_inf_mag / 1000) ** 2 # Convert to km/s then square # Delta-V estimate (simplified) delta_v_ms = v_inf_mag porkchop_data.append( { "departure_date": dep_date, "arrival_date": arr_date, "time_of_flight_days": tof_days, "c3_km2_s2": c3, "delta_v_ms": delta_v_ms, "departure_position_m": r1, "arrival_position_m": r2, "transfer_feasible": c3 < 100.0, # Reasonable C3 limit } ) except Exception as e: # Add failed transfer case porkchop_data.append( { "departure_date": dep_date, "arrival_date": arr_date, "time_of_flight_days": tof_days, "c3_km2_s2": float("inf"), "delta_v_ms": float("inf"), "departure_position_m": r1, "arrival_position_m": r2, "transfer_feasible": False, "error": str(e), } ) # Find optimal transfer feasible_transfers = [t for t in porkchop_data if t["transfer_feasible"]] optimal_transfer = None if feasible_transfers: # Find minimum C3 transfer optimal_transfer = min(feasible_transfers, key=lambda x: x["c3_km2_s2"]) # Generate summary statistics if feasible_transfers: c3_values = [t["c3_km2_s2"] for t in feasible_transfers] tof_values = [t["time_of_flight_days"] for t in feasible_transfers] summary_stats = { "feasible_transfers": len(feasible_transfers), "total_transfers_computed": len(porkchop_data), "min_c3_km2_s2": min(c3_values), "max_c3_km2_s2": max(c3_values), "mean_c3_km2_s2": sum(c3_values) / len(c3_values), "min_tof_days": min(tof_values), "max_tof_days": max(tof_values), "mean_tof_days": sum(tof_values) / len(tof_values), } else: summary_stats = { "feasible_transfers": 0, "total_transfers_computed": len(porkchop_data), "min_c3_km2_s2": None, "max_c3_km2_s2": None, "mean_c3_km2_s2": None, "min_tof_days": min_tof_days, "max_tof_days": max_tof_days, "mean_tof_days": (min_tof_days + max_tof_days) / 2, } return { "departure_body": departure_body, "arrival_body": arrival_body, "transfer_grid": porkchop_data, "optimal_transfer": optimal_transfer, "summary_statistics": summary_stats, "constraints": { "min_tof_days": min_tof_days, "max_tof_days": max_tof_days, "max_feasible_c3_km2_s2": 100.0, }, "note": "Simplified analysis using circular planetary orbits. Use SPICE kernels for production applications.", } def get_ephemeris_position( body: str, epoch_utc: str, frame: str = "J2000" ) -> dict[str, Any]: """ Get ephemeris position of celestial body (stub implementation). Args: body: Celestial body name (Earth, Mars, Venus, etc.) epoch_utc: UTC epoch in ISO format frame: Reference frame (default: J2000) Returns: Position and velocity vectors, with accuracy notes Note: This is a simplified implementation using circular orbits. For production use, install SPICE kernels and use spiceypy or similar. """ # Import required for date parsing import datetime try: dt = datetime.datetime.fromisoformat(epoch_utc.replace("Z", "+00:00")) except Exception: dt = datetime.datetime.fromisoformat(epoch_utc) # Days since J2000.0 epoch j2000 = datetime.datetime(2000, 1, 1, 12, 0, 0) days_since_j2000 = (dt - j2000).total_seconds() / 86400 # Simplified planetary orbital elements (circular orbits) orbital_data = { "earth": { "semi_major_axis_au": 1.0, "period_days": 365.25, "inclination_deg": 0.0, }, "mars": { "semi_major_axis_au": 1.524, "period_days": 686.98, "inclination_deg": 1.85, }, "venus": { "semi_major_axis_au": 0.723, "period_days": 224.70, "inclination_deg": 3.39, }, "jupiter": { "semi_major_axis_au": 5.204, "period_days": 4332.82, "inclination_deg": 1.30, }, } body_lower = body.lower() if body_lower not in orbital_data: return { "error": f"Body '{body}' not supported in simplified ephemeris", "supported_bodies": list(orbital_data.keys()), "recommendation": "Use SPICE kernels with spiceypy for accurate ephemeris data", } data = orbital_data[body_lower] # Calculate mean anomaly mean_anomaly = 2 * math.pi * days_since_j2000 / data["period_days"] # Position in heliocentric frame (simplified circular orbit) r_au = data["semi_major_axis_au"] x_au = r_au * math.cos(mean_anomaly) y_au = r_au * math.sin(mean_anomaly) z_au = 0.0 # Simplified: ignore inclination # Convert to meters AU_TO_M = 1.496e11 position_m = [x_au * AU_TO_M, y_au * AU_TO_M, z_au * AU_TO_M] # Velocity (circular orbit) mu_sun = 1.32712440018e20 # m³/s² r_m = r_au * AU_TO_M v_circular = math.sqrt(mu_sun / r_m) # m/s velocity_ms = [ -v_circular * math.sin(mean_anomaly), v_circular * math.cos(mean_anomaly), 0.0, ] return { "body": body, "epoch_utc": epoch_utc, "frame": frame, "position_m": position_m, "velocity_ms": velocity_ms, "accuracy": "low", "model": "circular_orbit_approximation", "warning": "This is a simplified model. Use SPICE kernels for mission-critical applications.", "spice_recommendation": { "kernels_needed": ["de430.bsp", "pck00010.tpc", "naif0012.tls"], "library": "spiceypy", "installation": "pip install spiceypy", }, } # Optional SPICE integration (if available) SPICE_AVAILABLE = False try: import spiceypy as spice SPICE_AVAILABLE = True def get_ephemeris_position_spice( body: str, epoch_utc: str, frame: str = "J2000" ) -> dict[str, Any]: """ Get accurate ephemeris position using SPICE kernels (if loaded). Args: body: SPICE body name or ID epoch_utc: UTC epoch in ISO format frame: SPICE reference frame Returns: High-accuracy position and velocity vectors """ try: # Convert ISO time to ET (Ephemeris Time) et = spice.str2et(epoch_utc) # Get state vector (position and velocity) state, light_time = spice.spkezr(body, et, frame, "NONE", "SUN") # Convert km to m, km/s to m/s position_m = [state[i] * 1000 for i in range(3)] velocity_ms = [state[i + 3] * 1000 for i in range(3)] return { "body": body, "epoch_utc": epoch_utc, "frame": frame, "position_m": position_m, "velocity_ms": velocity_ms, "light_time_s": light_time, "accuracy": "high", "model": "SPICE_kernels", "source": "JPL_ephemeris", } except Exception as e: return { "error": f"SPICE error: {str(e)}", "fallback_available": True, "recommendation": "Check SPICE kernel loading and body names", } except ImportError: def get_ephemeris_position_spice( body: str, epoch_utc: str, frame: str = "J2000" ) -> dict[str, Any]: """Placeholder when SPICE is not available.""" return { "error": "SPICE not available", "install_command": "pip install spiceypy", "kernel_sources": [ "https://naif.jpl.nasa.gov/naif/data_generic.html", "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/", ], } try: from . import update_availability update_availability("orbits", True, {"spice_available": SPICE_AVAILABLE}) except ImportError: pass