ENH: enhance modeling and usage of orbital features.

Author: phmbressan

docs/notebooks/orbital_flight_playground_hohman_geo.ipynb

@@ -6,3 +6,4 @@ requests
 pytz
 simplekml
 dill
+spiceypy>=6.0

docs/notebooks/orbital_flight_playground_ssto.ipynb

@@ -2,39 +2,48 @@
 from .control import _Controller
 from .environment import (
     Atmosphere,
-    AnalyticalEphemeris,
+    AtmosphereLayer,
     AtmosphericState,
-    CelestialBody,
     DefaultGravity,
+    Earth,
+    EarthRadiationPressure,
     Environment,
     EnvironmentAnalysis,
-    ExponentialAtmosphere,
+    ExponentialAtmosphereLayer,
     Gravity,
-    HarrisPriesterAtmosphere,
-    LayeredAtmosphere,
-    NRLMSISE00,
+    HarrisPriesterAtmosphereLayer,
+    NRLMSISE00AtmosphereLayer,
+    SomiglianaGravity,
+    Space,
     SphericalGravity,
     SphericalHarmonicGravity,
     SpiceEphemeris,
-    VacuumAtmosphere,
+    ThirdBody,
     VerticalGravity,
+    ZeroAtmosphereLayer,
     ZeroGravity,
     ZonalGravity,
 )
 from .mathutils import (
+    FLAT_WGS84,
+    GRS80,
+    NUMBA_AVAILABLE,
+    SIMPLE_WGS84,
+    WGS72,
+    WGS84,
+    Datum,
     EarthDatum,
     Epoch,
+    FlatEarthDatum,
     FlightState,
+    FrameVector,
     Function,
-    NUMBA_AVAILABLE,
     OrbitalElements,
     PiecewiseFunction,
     ReferenceFrame,
-    WGS84,
+    SimpleDatum,
     VectorFunction,
     funcify_method,
-    gcrf_to_rtn_matrix,
-    itrf_to_topocentric,
     numbify,
     reset_funcified_methods,
 )
@@ -75,23 +84,21 @@
     PointMassRocket,
     RailButtons,
     Rocket,
+    Spacecraft,
     Tail,
     TrapezoidalFin,
     TrapezoidalFins,
+    Vehicle,
 )
 from .sensitivity import SensitivityModel
 from .sensors import Accelerometer, Barometer, GnssReceiver, Gyroscope
 from .simulation import (
-    EarthRadiationPressure,
     Event,
     Flight,
     FlightOrbit,
+    Mission,
     MonteCarlo,
     MultivariateRejectionSampler,
-    PlanetaryRadiationPressure,
-    RelativisticCorrection,
-    SolarRadiationPressure,
-    ThirdBodyGravity,
 )
 from .stochastic import (
     CustomSampler,

requirements.txt

@@ -140,6 +140,12 @@ def object_hook(self, obj):
 
 
 def set_minimal_flight_attributes(flight, obj):
+    from rocketpy.environment.environment import Environment
+    from rocketpy.environment.models import Earth
+    from rocketpy.mathutils.epoch import Epoch
+    from rocketpy.mathutils.flight_state import FlightState
+    from rocketpy.mathutils.reference_frame import FlatEarthDatum
+
     attributes = (
         "rocket",
         "env",
@@ -200,6 +206,42 @@ def set_minimal_flight_attributes(flight, obj):
     flight.ode_solver = obj.get("ode_solver", "LSODA")
     flight.simulation_mode = obj.get("simulation_mode", "6DOF")
 
+    # Restore the shared Flight contract for decoded objects that deliberately
+    # bypass ``Flight.__init__``. Frame and domain ownership can be derived
+    # from the serialized environment, avoiding duplicated datum state.
+    flight.initial_environment = flight.env
+    flight.earth = flight.env if isinstance(flight.env, Earth) else None
+    flight.space = obj.get("space")
+    flight.vehicle = flight.rocket
+    if isinstance(flight.env, Earth):
+        flight.datum = flight.env.datum
+    elif isinstance(flight.env, Environment):
+        weather_datum = flight.env.earth_datum
+        flight.datum = FlatEarthDatum(
+            name=f"Flat {weather_datum.name}",
+            semi_major_axis=weather_datum.semi_major_axis,
+            flattening=weather_datum.flattening,
+            angular_velocity=weather_datum.angular_velocity,
+            gravitational_parameter=weather_datum.gravitational_parameter,
+        )
+    else:  # pragma: no cover - defensive support for old custom environments
+        flight.datum = FlatEarthDatum()
+    flight.reference_frame = flight.datum.integration_frame
+    flight.initial_state = obj.get("initial_state")
+    flight.start_epoch = (
+        flight.initial_state.epoch
+        if isinstance(flight.initial_state, FlightState)
+        else getattr(flight.env, "epoch", Epoch.relative_origin())
+    )
+    flight.mission_time_offset = (
+        flight.initial_state.elapsed_time
+        if isinstance(flight.initial_state, FlightState)
+        else 0.0
+    )
+    flight.forces = list(obj.get("forces") or ())
+    flight._state_type = FlightState
+    flight._launch_with_earth = False
+
     # Controllers, sensors and custom events are derived by Flight.__init__
     # from the (already-decoded) rocket and constructor args; they are not
     # serialized separately. Mirror that derivation here so prints/plots that
@@ -208,6 +250,10 @@ def set_minimal_flight_attributes(flight, obj):
     flight._controllers = getattr(flight.rocket, "_controllers", [])[:]
     flight.sensors = flight.rocket.sensors.get_components()
     flight.sensors_by_name = flight.rocket.sensors_by_name
+    flight.sensor_data = {
+        sensor: list(getattr(sensor, "measured_data", ()))
+        for sensor in flight.sensors
+    }
     # TODO: custom_events are lost when loading from .rpy because they hold
     # user-defined callables that are not currently serialized. Add proper
     # serialization/deserialization for custom_events and restore them here

rocketpy/__init__.py

@@ -3,47 +3,56 @@
 considered private and should be used with caution.
 """
 
+from .albedo import EarthRadiationPressure
 from .atmosphere import (
     Atmosphere,
+    AtmosphereLayer,
     AtmosphericState,
-    ExponentialAtmosphere,
-    FunctionAtmosphere,
-    HarrisPriesterAtmosphere,
-    LayeredAtmosphere,
-    NRLMSISE00,
-    VacuumAtmosphere,
+    ExponentialAtmosphereLayer,
+    FunctionAtmosphereLayer,
+    HarrisPriesterAtmosphereLayer,
+    NRLMSISE00AtmosphereLayer,
+    ZeroAtmosphereLayer,
 )
-from .celestial_body import AnalyticalEphemeris, CelestialBody, SpiceEphemeris
 from .environment import Environment
 from .environment_analysis import EnvironmentAnalysis
 from .gravity import (
     DefaultGravity,
     Gravity,
+    SomiglianaGravity,
     SphericalGravity,
     SphericalHarmonicGravity,
     VerticalGravity,
     ZeroGravity,
     ZonalGravity,
 )
+from .models import Earth, Space
+from .third_body import (
+    SpiceEphemeris,
+    ThirdBody,
+)
 
 __all__ = [
-    "AnalyticalEphemeris",
     "Atmosphere",
+    "AtmosphereLayer",
     "AtmosphericState",
-    "CelestialBody",
     "DefaultGravity",
     "Environment",
+    "Earth",
+    "EarthRadiationPressure",
     "EnvironmentAnalysis",
-    "ExponentialAtmosphere",
-    "FunctionAtmosphere",
+    "ExponentialAtmosphereLayer",
+    "FunctionAtmosphereLayer",
     "Gravity",
-    "HarrisPriesterAtmosphere",
-    "LayeredAtmosphere",
-    "NRLMSISE00",
+    "HarrisPriesterAtmosphereLayer",
+    "NRLMSISE00AtmosphereLayer",
     "SphericalGravity",
+    "SomiglianaGravity",
     "SphericalHarmonicGravity",
     "SpiceEphemeris",
-    "VacuumAtmosphere",
+    "Space",
+    "ThirdBody",
+    "ZeroAtmosphereLayer",
     "VerticalGravity",
     "ZeroGravity",
     "ZonalGravity",

rocketpy/_encoders.py

@@ -0,0 +1,157 @@
+"""Earth radiation-pressure models."""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+
+from rocketpy.environment.third_body import ASTRONOMICAL_UNIT, ThirdBody
+from rocketpy.mathutils.compilation import numbify
+from rocketpy.mathutils.frame_vector import FrameVector
+from rocketpy.mathutils.reference_frame import ReferenceFrame
+
+
+class EarthRadiationPressure:
+    """Knocke finite-element Earth albedo and infrared pressure model.
+
+    The visible terrestrial cap is integrated over a fixed 10-by-20 surface
+    grid. Reflected shortwave radiation and emitted longwave radiation are
+    both included.
+    """
+
+    name = "earth_radiation_pressure"
+
+    def __init__(self, sun):
+        if not isinstance(sun, ThirdBody):
+            raise TypeError("sun must be an explicit ThirdBody.")
+        self.sun = sun
+
+    def acceleration(self, epoch, state, vehicle, earth):
+        """Return Earth radiation-pressure acceleration in GCRF."""
+        sun_position = self.sun.position(epoch)
+        incident_direction = FrameVector(-state.position, ReferenceFrame.GCRF)
+        coefficient = vehicle.evaluate_radiation_coefficient(
+            epoch, state, incident_direction
+        )
+        area = vehicle.evaluate_radiation_area(epoch, state, incident_direction)
+        mass = vehicle.mass_at(state.elapsed_time)
+        return self._evaluate_knocke_erp(
+            epoch.jd_tdb,
+            np.array(state.position, dtype=float, copy=True),
+            np.array(sun_position, dtype=float, copy=True),
+            float(area),
+            float(coefficient),
+            mass,
+        )
+
+    @staticmethod
+    @numbify(
+        nopython=True,
+        signature=(
+            "float64[:](float64, float64[:], float64[:], float64, float64, float64)"
+        ),
+    )
+    def _evaluate_knocke_erp(jd, r_sat, r_sun, area_sat, coefficient, mass):
+        """Evaluate the finite-element Earth-radiation integration."""
+        earth_radius = 6378136.0
+        satellite_norm = math.sqrt(r_sat[0] ** 2 + r_sat[1] ** 2 + r_sat[2] ** 2)
+        sun_norm = math.sqrt(r_sun[0] ** 2 + r_sun[1] ** 2 + r_sun[2] ** 2)
+        if satellite_norm <= earth_radius:
+            return np.zeros(3, dtype=np.float64)
+
+        sun_unit = r_sun / sun_norm
+        satellite_unit = r_sat / satellite_norm
+        solar_pressure = 4.56e-6 / (sun_norm / ASTRONOMICAL_UNIT) ** 2
+        cap_angle = math.asin(earth_radius / satellite_norm)
+        local_z = satellite_unit
+        auxiliary = np.zeros(3, dtype=np.float64)
+        if abs(local_z[2]) < 0.99:
+            auxiliary[2] = 1.0
+        else:
+            auxiliary[0] = 1.0
+        local_x = np.array(
+            [
+                local_z[1] * auxiliary[2] - local_z[2] * auxiliary[1],
+                local_z[2] * auxiliary[0] - local_z[0] * auxiliary[2],
+                local_z[0] * auxiliary[1] - local_z[1] * auxiliary[0],
+            ],
+            dtype=np.float64,
+        )
+        local_x /= math.sqrt(local_x[0] ** 2 + local_x[1] ** 2 + local_x[2] ** 2)
+        local_y = np.array(
+            [
+                local_z[1] * local_x[2] - local_z[2] * local_x[1],
+                local_z[2] * local_x[0] - local_z[0] * local_x[2],
+                local_z[0] * local_x[1] - local_z[1] * local_x[0],
+            ],
+            dtype=np.float64,
+        )
+
+        cosine_phase = math.cos(math.tau / 365.25 * (jd - 2451535.0))
+        rings = 10
+        sectors = 20
+        delta_alpha = cap_angle / rings
+        delta_beta = math.tau / sectors
+        pressure_vector = np.zeros(3, dtype=np.float64)
+        for ring in range(rings):
+            alpha = (ring + 0.5) * delta_alpha
+            sin_alpha = math.sin(alpha)
+            cos_alpha = math.cos(alpha)
+            area_element = earth_radius**2 * sin_alpha * delta_alpha * delta_beta
+            for sector in range(sectors):
+                beta = (sector + 0.5) * delta_beta
+                normal = (
+                    sin_alpha * math.cos(beta) * local_x
+                    + sin_alpha * math.sin(beta) * local_y
+                    + cos_alpha * local_z
+                )
+                element_to_satellite = r_sat - earth_radius * normal
+                distance = math.sqrt(
+                    element_to_satellite[0] ** 2
+                    + element_to_satellite[1] ** 2
+                    + element_to_satellite[2] ** 2
+                )
+                direction = element_to_satellite / distance
+                satellite_cosine = (
+                    normal[0] * direction[0]
+                    + normal[1] * direction[1]
+                    + normal[2] * direction[2]
+                )
+                if satellite_cosine <= 0.0:
+                    continue
+                latitude_sine = normal[2]
+                legendre_p2 = 0.5 * (3.0 * latitude_sine**2 - 1.0)
+                emissivity = max(
+                    0.0,
+                    min(
+                        1.0,
+                        0.68 - 0.07 * cosine_phase * latitude_sine - 0.18 * legendre_p2,
+                    ),
+                )
+                surface_pressure = emissivity * solar_pressure / 4.0
+                sun_cosine = (
+                    normal[0] * sun_unit[0]
+                    + normal[1] * sun_unit[1]
+                    + normal[2] * sun_unit[2]
+                )
+                if sun_cosine > 0.0:
+                    albedo = max(
+                        0.0,
+                        min(
+                            1.0,
+                            0.34
+                            + 0.10 * cosine_phase * latitude_sine
+                            + 0.29 * legendre_p2,
+                        ),
+                    )
+                    surface_pressure += albedo * solar_pressure * sun_cosine
+                pressure_vector += (
+                    surface_pressure
+                    / math.pi
+                    * area_element
+                    * satellite_cosine
+                    / distance**2
+                    * direction
+                )
+        return pressure_vector * (coefficient * area_sat / mass)