Satellite transfer

A satellite moved between Kepler orbits of equal angular momentum by radial thrust.

• Name: satellite_transfer
• Field: engineering/control
• Formalism: ordinary differential equations
• References: Dilão et al. 2022; Dilão 2023

Factory

@register_model(
    field='engineering/control',
    tags=('control', 'orbits', 'celestial mechanics', 'ode'),
    references=('Dilão et al. 2022', 'Dilão 2023'),
    defaults={'until': 20.0, 'dt': 0.05, 'solver': RK4Solver(dt=0.005)},
    observables=lambda system: [ObservableSpec('energy', _energy_of_satellite)],
)
def satellite_transfer(
    s: float = 0.314,
    velocity: float = 0.0,
    theta: float = 0.0,
    *,
    Lz: float = 0.7,
    epsilon: float = 0.1,
    H_final: float = -0.2,
) -> System:
    """Return a satellite moved between Kepler orbits of equal angular momentum by radial thrust.

    In normalised polar coordinates the satellite follows `s'' = L_z^2 / s^3 - 1 / s^2 - ε sigma` and
    `θ' = L_z / s^2`. Radial thrust keeps the angular momentum and changes the effective energy
    `H = s'^2 / 2 + L_z^2 / (2 s^2) - 1 / s` at the rate `-ε sigma s'`, so the controller pushes inwards, `sigma = 1`,
    while the satellite falls and outwards, `sigma = -1`, while it rises, until `H` reaches the final energy. From the
    periapsis of the orbit with `H = -0.7`, the transfer to `H = -0.2` takes about 9.81.

    Args:
        s:
            The initial distance to the planet.

        velocity:
            The initial radial velocity.

        theta:
            The initial angle.

        Lz:
            The angular momentum.

        epsilon:
            The thrust.

        H_final:
            The effective energy of the final orbit.

    Returns:
        system:
            The model, with the distance `satellite.s`, the radial velocity `satellite.s_d1`, the angle
            `satellite.theta`, and the observable `energy`.

    Examples:
        >>> from skcomplex.models.engineering.control import satellite_transfer
        >>> satellite_transfer().elements[0].Lz
        0.7
    """
    satellite = SatelliteElement(
        'satellite',
        s=s,
        s_d1=velocity,
        theta=theta,
        Lz=Lz,
        epsilon=epsilon,
        H_final=H_final,
    )
    return System(elements=[satellite])

Components

class SatelliteElement(Element):
    """Satellite on a planar Kepler orbit whose radial thrust moves it to an orbit of a final energy."""

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        s: Annotated[object, StateSpec(domain=RealSet.open(0, oo), order=2)] = 0.314,
        s_d1: float = 0.0,
        theta: Annotated[object, StateSpec()] = 0.0,
        Lz: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 0.7,
        epsilon: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        H_final: Annotated[object, ParameterSpec()] = -0.2,
    ) -> None:
        super().__init__(label, visual=visual)
        self.s = s
        self.s_d1 = s_d1
        self.theta = theta
        self.Lz = Lz
        self.epsilon = epsilon
        self.H_final = H_final

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the radial rate under thrust and the rate of the angle."""
        return [
            RateLaw(s.s, s.Lz**2 / s.s**3 - 1 / s.s**2 - s.epsilon * SatelliteElement._thrust_direction(s)),
            RateLaw(s.theta, s.Lz / s.s**2),
        ]

    @staticmethod
    def _satellite_energy(s: SimpleNamespace) -> Expression:
        """Return the total effective energy of the satellite, `s'^2 / 2 + L_z^2 / (2 s^2) - 1 / s`."""
        return derivative(s.s, 1) ** 2 / 2 + s.Lz**2 / (2 * s.s**2) - 1 / s.s

    @staticmethod
    def _thrust_direction(s: SimpleNamespace) -> Expression:
        """Return the thrust `sigma`: inwards while the satellite falls, outwards while it rises, to the end energy."""
        return where(SatelliteElement._satellite_energy(s) < s.H_final, where(derivative(s.s, 1) < 0, 1, -1), 0)

Equations

satellite.s: diff(s_satellite(t), t, t) == -epsilon_satellite*where(1/2*diff(s_satellite(t), t)^2 + 1/2*Lz_satellite^2/s_satellite(t)^2 - 1/s_satellite(t) < H_final_satellite, where(diff(s_satellite(t), t) < 0, 1, -1), 0) + Lz_satellite^2/s_satellite(t)^3 - 1/s_satellite(t)^2
satellite.theta: diff(theta_satellite(t), t) == Lz_satellite/s_satellite(t)^2