Source code for ADCS.orbits.helpers.orbit_factory

# Without __all__, autodoc treats the imported Orbit / Orbital_State below as
# members of this module and documents them a second time, which Sphinx reports
# as a duplicate object description -- fatal under the docs build's -W.
__all__ = ["create_random_circular_os", "create_random_circular_orbit"]

import numpy as np
from typing import Optional

from ADCS.orbits.ephemeris import Ephemeris
from ADCS.orbits.orbital_state import Orbital_State
from ADCS.orbits.orbit import Orbit
from ADCS.orbits.universal_constants import TimeConstants, EarthConstants
from ADCS.helpers.math_helpers import normalize



def create_random_circular_os(
    radius_km: float,
    J2000: float = 0.22,
    ephem: Optional[Ephemeris] = None,
    rng: Optional[np.random.Generator] = None
) -> Orbital_State:
    r"""
    Generates a random orbital state for a circular orbit with a specified radius.

    The position vector :math:`\mathbf{R}` is generated by sampling a random vector :math:`\mathbf{x}`
    from a standard normal distribution :math:`\mathcal{N}(0, I_3)` and normalizing it to lie on the unit sphere,
    then scaling by the radius :math:`r`:

    .. math::

        \hat{\mathbf{r}} = \frac{\mathbf{x}}{\|\mathbf{x}\|}

        \mathbf{R} = r \cdot \hat{\mathbf{r}}

    The velocity vector :math:`\mathbf{V}` is generated by creating a random tangent direction.
    A second random vector :math:`\mathbf{v}_{rand}` is sampled and projected onto the plane
    tangent to :math:`\hat{\mathbf{r}}`:

    .. math::

        \mathbf{v}_{tan} = \mathbf{v}_{rand} - (\mathbf{v}_{rand} \cdot \hat{\mathbf{r}}) \hat{\mathbf{r}}

        \hat{\mathbf{v}} = \frac{\mathbf{v}_{tan}}{\|\mathbf{v}_{tan}\|}

    The magnitude is set to the circular velocity :math:`v_{circ}`:

    .. math::

        v_{circ} = \sqrt{\frac{\mu}{r}}

        \mathbf{V} = v_{circ} \cdot \hat{\mathbf{v}}

    :param radius_km: The orbital radius in kilometers.
    :type radius_km: float
    :param J2000: The current time in J2000 centuries. Defaults to 0.22.
    :type J2000: float
    :param ephem: Optional ephemeris object. If None, a new one is created.
    :type ephem: :class:`~ADCS.orbits.ephemeris.Ephemeris`, optional
    :return: The generated orbital state.
    :rtype: :class:`~ADCS.orbits.orbital_state.Orbital_State`

    """
    if ephem is None:
        ephem = Ephemeris()
    
    # Fallback to default generator if none provided
    if rng is None:
        rng = np.random.default_rng()

    mu = EarthConstants.mu_e

    # Use rng.standard_normal instead of np.random.standard_normal
    r_hat = normalize(rng.standard_normal(3))
    R = radius_km * r_hat

    v_rand = rng.standard_normal(3)
    v_tan = v_rand - np.dot(v_rand, r_hat) * r_hat
    v_hat = normalize(v_tan)

    v_circ = float(np.sqrt(mu / radius_km))
    V = v_circ * v_hat

    return Orbital_State(
        ephem=ephem,
        J2000=J2000,
        R=R,
        V=V,
    )





def create_random_circular_orbit(
    radius_km: float,
    dt: float,
    tf: float,
    J2000: float = 0.22,
    fast: bool = False,
    rng: Optional[np.random.Generator] = None,
    zonal_J: int = 2,
) -> Orbit:
    r"""
    Creates an initialized Orbit object based on a random circular orbital state.

    This function utilizes :func:`~create_random_circular_os` to generate the initial
    conditions :math:`\mathbf{R}` and :math:`\mathbf{V}`. It then initializes the orbit propagator.

    The end time is calculated as:

    .. math::

        t_{end} = t_{start} + t_{f} \cdot C_{sec2cent}

    where :math:`C_{sec2cent}` is the conversion factor from seconds to centuries.

    :param radius_km: The orbital radius in kilometers.
    :type radius_km: float
    :param dt: The simulation time step in seconds.
    :type dt: float
    :param tf: The total simulation duration in seconds.
    :type tf: float
    :param J2000: The start time in J2000 centuries. Defaults to 0.22.
    :type J2000: float
    :param fast: Flag to enable fast propagation mode (reduced precision).
    :type fast: bool
    :param zonal_J: Highest zonal harmonic degree to include. ``0`` disables
        zonals, ``2`` includes only J2, and larger values include every zonal
        term up to that degree.
    :type zonal_J: int
    :return: The fully initialized orbit object.
    :rtype: :class:`~ADCS.orbits.orbit.Orbit`

    """
    os0 = create_random_circular_os(radius_km=radius_km, J2000=J2000, rng=rng)

    return Orbit(
        os0=os0,
        end_time=J2000 + tf * TimeConstants.sec2cent,
        dt=dt,
        fast=fast,
        verbose=False,
        zonal_J=zonal_J,
    )