Skip to content

dynamics

cosmica.dynamics

__all__ module-attribute 

__all__ = [
    "CircularSatelliteOrbit",
    "CircularSatelliteOrbitPropagator",
    "EllipticalSatelliteOrbit",
    "EllipticalSatelliteOrbitPropagator",
    "MOPCSatelliteKey",
    "MultiOrbitalPlaneConstellation",
    "SatelliteConstellation",
    "SatelliteOrbit",
    "SatelliteOrbitState",
    "get_sun_direction_eci",
    "make_satellite_orbit",
]

CircularSatelliteOrbit dataclass 

CircularSatelliteOrbit(
    *,
    semi_major_axis: float,
    inclination: float,
    raan: float,
    phase_at_epoch: float,
    epoch: datetime64
)

Bases: SatelliteOrbit

dcm_orbit_to_eci cached property 

dcm_orbit_to_eci: NDArray[floating]

epoch instance-attribute 

epoch: datetime64

inclination instance-attribute 

inclination: float

mean_motion cached property 

mean_motion: float

phase_at_epoch instance-attribute 

phase_at_epoch: float

raan instance-attribute 

raan: float

semi_major_axis instance-attribute 

semi_major_axis: float

__post_init__ 

__post_init__() -> None
Source code in src/cosmica/dynamics/orbit.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
def __post_init__(self) -> None:
    object.__setattr__(
        self,
        "_model",
        CircularSatelliteOrbitModel(
            semi_major_axis=self.semi_major_axis,
            inclination=self.inclination,
            raan=self.raan,
            phase_at_epoch=self.phase_at_epoch,
            epoch=self.epoch,
        ),
    )
    object.__setattr__(self, "_propagator", CircularSatelliteOrbitPropagator(self._model))

propagate 

propagate(time: NDArray[datetime64]) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
210
211
def propagate(self, time: npt.NDArray[np.datetime64]) -> SatelliteOrbitState:
    return self._propagator.propagate(time)

propagate_from_epoch 

propagate_from_epoch(
    time_from_epoch: NDArray[timedelta64],
) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
213
214
def propagate_from_epoch(self, time_from_epoch: npt.NDArray[np.timedelta64]) -> SatelliteOrbitState:
    return self._propagator.propagate_from_epoch(time_from_epoch)

CircularSatelliteOrbitPropagator 

CircularSatelliteOrbitPropagator(
    model: CircularSatelliteOrbitModel,
)

Bases: SatelliteOrbitPropagator

Source code in src/cosmica/dynamics/orbit.py 
67
68
def __init__(self, model: CircularSatelliteOrbitModel) -> None:
    self._model = model

dcm_orbit_to_eci cached property 

dcm_orbit_to_eci: NDArray[floating]

mean_motion cached property 

mean_motion: float

propagate 

propagate(time: NDArray[datetime64]) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
92
93
94
def propagate(self, time: npt.NDArray[np.datetime64]) -> SatelliteOrbitState:
    time_from_epoch = time - self._model.epoch
    return self.propagate_from_epoch(time_from_epoch)

propagate_from_epoch 

propagate_from_epoch(
    time_from_epoch: NDArray[timedelta64],
) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def propagate_from_epoch(self, time_from_epoch: npt.NDArray[np.timedelta64]) -> SatelliteOrbitState:
    time_from_epoch_in_seconds = time_from_epoch / np.timedelta64(1, "s")
    phase = self._model.phase_at_epoch + self.mean_motion * time_from_epoch_in_seconds
    in_plane_position = self._model.semi_major_axis * np.array([np.cos(phase), np.sin(phase), np.zeros_like(phase)])
    in_plane_velocity = (
        self.mean_motion
        * self._model.semi_major_axis
        * np.array([-np.sin(phase), np.cos(phase), np.zeros_like(phase)])
    )

    position_eci: npt.NDArray[np.floating] = rowwise_matmul(self.dcm_orbit_to_eci, in_plane_position.T)
    velocity_eci: npt.NDArray[np.floating] = rowwise_matmul(self.dcm_orbit_to_eci, in_plane_velocity.T)

    return SatelliteOrbitState(position_eci=position_eci, velocity_eci=velocity_eci)

EllipticalSatelliteOrbit dataclass 

EllipticalSatelliteOrbit(
    *,
    semi_major_axis: float,
    inclination: float,
    raan: float,
    phase_at_epoch: float,
    epoch: datetime64,
    satnum: int,
    gravity_model: GravityModel = WGS84,
    drag_coeff: float,
    eccentricity: float,
    argpo: float
)

Bases: SatelliteOrbit

argpo instance-attribute 

argpo: float

drag_coeff instance-attribute 

drag_coeff: float

eccentricity instance-attribute 

eccentricity: float

epoch instance-attribute 

epoch: datetime64

gravity_model class-attribute instance-attribute 

gravity_model: GravityModel = WGS84

inclination instance-attribute 

inclination: float

mean_motion cached property 

mean_motion: float

phase_at_epoch instance-attribute 

phase_at_epoch: float

raan instance-attribute 

raan: float

satellite cached property 

satellite: EarthSatellite

satnum instance-attribute 

satnum: int

semi_major_axis instance-attribute 

semi_major_axis: float

ts cached property 

ts: Timescale

__post_init__ 

__post_init__() -> None
Source code in src/cosmica/dynamics/orbit.py 
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
def __post_init__(self) -> None:
    object.__setattr__(
        self,
        "_model",
        EllipticalSatelliteOrbitModel(
            semi_major_axis=self.semi_major_axis,
            inclination=self.inclination,
            raan=self.raan,
            phase_at_epoch=self.phase_at_epoch,
            epoch=self.epoch,
            satnum=self.satnum,
            gravity_model=self.gravity_model,
            drag_coeff=self.drag_coeff,
            eccentricity=self.eccentricity,
            argpo=self.argpo,
        ),
    )
    object.__setattr__(self, "_propagator", EllipticalSatelliteOrbitPropagator(self._model))

datetime64_utc_to_skytime 

datetime64_utc_to_skytime(
    t_datetime64: NDArray[datetime64],
) -> Time
Source code in src/cosmica/dynamics/orbit.py 
263
264
def datetime64_utc_to_skytime(self, t_datetime64: npt.NDArray[np.datetime64]) -> Time:
    return self._propagator.datetime64_utc_to_skytime(t_datetime64)

propagate 

propagate(time: NDArray[datetime64]) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
270
271
def propagate(self, time: npt.NDArray[np.datetime64]) -> SatelliteOrbitState:
    return self._propagator.propagate(time)

EllipticalSatelliteOrbitPropagator 

EllipticalSatelliteOrbitPropagator(
    model: EllipticalSatelliteOrbitModel,
    reference_frame: ReferenceFrame = "gcrs",
)

Bases: SatelliteOrbitPropagator

Source code in src/cosmica/dynamics/orbit.py 
113
114
115
116
117
118
119
120
def __init__(self, model: EllipticalSatelliteOrbitModel, reference_frame: ReferenceFrame = "gcrs") -> None:
    self._model = model
    rf = reference_frame.lower()
    if rf not in {"teme", "j2000", "gcrs"}:
        logger.error("Invalid reference frame: %s", reference_frame)
        msg = f"Invalid reference_frame: {reference_frame!r}"
        raise ValueError(msg)
    self.reference_frame = rf

mean_motion cached property 

mean_motion: float

reference_frame instance-attribute 

reference_frame = rf

satellite cached property 

satellite: EarthSatellite

ts cached property 

ts: Timescale

datetime64_utc_to_skytime 

datetime64_utc_to_skytime(
    t_datetime64: NDArray[datetime64],
) -> Time
Source code in src/cosmica/dynamics/orbit.py 
130
131
132
def datetime64_utc_to_skytime(self, t_datetime64: npt.NDArray[np.datetime64]) -> Time:
    t_datetimes = [time.astype(datetime).replace(tzinfo=utc) for time in t_datetime64]
    return self.ts.from_datetimes(t_datetimes)

propagate 

propagate(time: NDArray[datetime64]) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
157
158
159
160
161
162
163
164
165
166
167
@override
def propagate(self, time: npt.NDArray[np.datetime64]) -> SatelliteOrbitState:
    time_from_epoch = self.datetime64_utc_to_skytime(time)
    if self.reference_frame in {"j2000", "gcrs"}:
        geoc = self.satellite.at(time_from_epoch)  # Geocentric position in GCRS/ICRF axes
        position_eci_km = geoc.position.km
        velocity_eci_km_s = geoc.velocity.km_per_s
        return SatelliteOrbitState(position_eci=position_eci_km.T * 1e3, velocity_eci=velocity_eci_km_s.T * 1e3)
    else:  # TEME
        [position_eci, velocity_eci, _] = self.satellite._position_and_velocity_TEME_km(time_from_epoch)  # noqa: SLF001
        return SatelliteOrbitState(position_eci=position_eci.T * 1e3, velocity_eci=velocity_eci.T * 1e3)

MOPCSatelliteKey

Bases: NamedTuple

plane_id instance-attribute 

plane_id: int

satellite_id instance-attribute 

satellite_id: int

__str__ 

__str__() -> str
Source code in src/cosmica/dynamics/constellation.py 
50
51
def __str__(self) -> str:
    return f"({self.plane_id},{self.satellite_id})"

MultiOrbitalPlaneConstellation 

MultiOrbitalPlaneConstellation(
    satellite_orbits: Mapping[
        ConstellationSatellite[MOPCSatelliteKey], TOrbit
    ],
)

Bases: SatelliteConstellation[MOPCSatelliteKey, TOrbit]

A constellation of satellites in multiple orbital planes.

The satellite key is a tuple of the plane ID and the satellite ID.

Source code in src/cosmica/dynamics/constellation.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def __init__(
    self,
    satellite_orbits: Mapping[ConstellationSatellite[MOPCSatelliteKey], TOrbit],
) -> None:
    self.satellite_orbits = satellite_orbits

    self.satellites = tuple(self.satellite_orbits.keys())
    self.plane_ids = sorted({sat.id.plane_id for sat in self.satellites})
    self.satellite_ids = sorted({sat.id.satellite_id for sat in self.satellites})

    self.plane_id_to_satellites: dict[int, list[ConstellationSatellite[MOPCSatelliteKey]]] = {
        plane_id: [sat for sat in self.satellites if sat.id.plane_id == plane_id] for plane_id in self.plane_ids
    }

    # Check if all planes have the same number of satellites
    self._all_planes_have_same_n_satellites = (
        len({len(self.plane_id_to_satellites[plane_id]) for plane_id in self.plane_ids}) == 1
    )
    if not self._all_planes_have_same_n_satellites:
        logger.warning("Not all planes have the same number of satellites.")

n_satellites_per_plane property 

n_satellites_per_plane: int

Number of satellites per orbital plane.

It is assumed that all planes have the same number of satellites. Otherwise, it will raise an exception.

plane_id_to_satellites instance-attribute 

plane_id_to_satellites: dict[
    int, list[ConstellationSatellite[MOPCSatelliteKey]]
] = {
    plane_id: [
        sat for sat in (satellites) if plane_id == plane_id
    ]
    for plane_id in (plane_ids)
}

plane_ids instance-attribute 

plane_ids = sorted({(plane_id) for sat in (satellites)})

satellite_ids instance-attribute 

satellite_ids = sorted(
    {(satellite_id) for sat in (satellites)}
)

satellite_orbits instance-attribute 

satellite_orbits = satellite_orbits

satellites instance-attribute 

satellites = tuple(keys())

from_toml_file classmethod 

from_toml_file(
    toml_file_path: Path | str,
) -> MultiOrbitalPlaneConstellation
Source code in src/cosmica/dynamics/constellation.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
@classmethod
@deprecated("Construction of objects from TOML files is deprecated and will be removed in future versions.")
def from_toml_file(cls, toml_file_path: Path | str) -> MultiOrbitalPlaneConstellation:
    toml_file_path = Path(toml_file_path)
    with toml_file_path.open("rb") as f:
        toml_data = tomllib.load(f)

    epoch = np.datetime64(toml_data["epoch"].astimezone(tz=UTC).replace(tzinfo=None))

    @deprecated("Construction of objects from TOML files is deprecated and will be removed in future versions.")
    def parse_satellite_item(
        item: dict[str, Any],
    ) -> tuple[ConstellationSatellite[MOPCSatelliteKey], TOrbit]:
        plane_id = item.pop("plane_id")
        satellite_id = item.pop("id")

        # Convert degrees to radians
        item["semi_major_axis"] = item.pop("sma_m")
        item["inclination"] = np.radians(item.pop("inc_deg"))
        item["raan"] = np.radians(item.pop("raan_deg"))
        item["phase_at_epoch"] = np.radians(item.pop("phase_at_epoch_deg"))
        item["epoch"] = epoch

        orbit_type = item.pop("orbit_type")
        return ConstellationSatellite(  # type: ignore[return-value]
            id=MOPCSatelliteKey(plane_id, satellite_id),
        ), make_satellite_orbit(
            orbit_type=orbit_type,
            **item,
        )

    satellite_orbits = dict(map(parse_satellite_item, toml_data["satellites"]))

    return cls(satellite_orbits)

propagate 

propagate(
    time: NDArray[datetime64],
) -> dict[ConstellationSatellite[T], SatelliteOrbitState]
Source code in src/cosmica/dynamics/constellation.py 
42
43
def propagate(self, time: npt.NDArray[np.datetime64]) -> dict[ConstellationSatellite[T], SatelliteOrbitState]:
    return {sat: orbit.propagate(time) for sat, orbit in self.satellite_orbits.items()}

SatelliteConstellation

Bases: ABC

A constellation of satellites.

satellite_orbits is a dictionary mapping satellite keys to satellite orbits.

satellite_orbits instance-attribute 

satellite_orbits: Mapping[
    ConstellationSatellite[T], TOrbit
]

satellites instance-attribute 

satellites: Sequence[ConstellationSatellite[T]]

propagate 

propagate(
    time: NDArray[datetime64],
) -> dict[ConstellationSatellite[T], SatelliteOrbitState]
Source code in src/cosmica/dynamics/constellation.py 
42
43
def propagate(self, time: npt.NDArray[np.datetime64]) -> dict[ConstellationSatellite[T], SatelliteOrbitState]:
    return {sat: orbit.propagate(time) for sat, orbit in self.satellite_orbits.items()}

SatelliteOrbit

Bases: ABC

propagate abstractmethod 

propagate(time: NDArray[datetime64]) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
172
173
@abstractmethod
def propagate(self, time: npt.NDArray[np.datetime64]) -> SatelliteOrbitState: ...

SatelliteOrbitState dataclass 

SatelliteOrbitState(
    *,
    position_eci: NDArray[floating],
    velocity_eci: NDArray[floating]
)

position_eci instance-attribute 

position_eci: NDArray[floating]

velocity_eci instance-attribute 

velocity_eci: NDArray[floating]

__getitem__ 

__getitem__(item: int | slice) -> SatelliteOrbitState
Source code in src/cosmica/dynamics/orbit.py 
51
52
53
54
55
def __getitem__(self, item: int | slice) -> SatelliteOrbitState:
    return SatelliteOrbitState(
        position_eci=self.position_eci[item],
        velocity_eci=self.velocity_eci[item],
    )

__post_init__ 

__post_init__() -> None
Source code in src/cosmica/dynamics/orbit.py 
43
44
45
46
47
48
49
def __post_init__(self) -> None:
    assert self.position_eci.shape[-1] == 3, self.position_eci.shape
    assert self.velocity_eci.shape[-1] == 3, self.velocity_eci.shape
    assert self.position_eci.shape[:-1] == self.velocity_eci.shape[:-1], (
        self.position_eci.shape,
        self.velocity_eci.shape,
    )

calc_position_ecef 

calc_position_ecef(
    dcm_eci_to_ecef: NDArray[floating],
) -> NDArray[floating]
Source code in src/cosmica/dynamics/orbit.py 
57
58
def calc_position_ecef(self, dcm_eci_to_ecef: npt.NDArray[np.floating]) -> npt.NDArray[np.floating]:
    return np.einsum("ijk,ik->ij", dcm_eci_to_ecef, self.position_eci)

get_sun_direction_eci 

get_sun_direction_eci(
    time: NDArray[datetime64],
) -> NDArray[float64]

get_sun_direction_eci(time: datetime64) -> NDArray[float64]

get_sun_direction_eci(
    time: datetime64 | NDArray[datetime64],
) -> NDArray[float64]

Sun direction approximation function.

The approximation considers an unperturbed motion of the Earth around the Sun with mean orbital elements that approximate the Sun's elliptic orbit wrt Earth around the year 2000. The algorithm is based on Oliver Montenbruck, Eberhard Grill; "Satellite Orbits: Models, Methods and Applications"; Springer-Verlag Berlin. The alterations to the original model follow JAXA's Attitude Control Handbook JERG-2-510-HB001.

PARAMETER DESCRIPTION
time

single element or array of UTC time values in np.datetime64 format.

TYPE: datetime64 | NDArray[datetime64]

RETURNS DESCRIPTION
NDArray[float64]

Corresponding Sun direction vector(s) (x,y,z) in eci. If the input is a single element, the output is a 1D array. If the input is an array, the output is a 2D array with the shape (len(time), 3).
Source code in src/cosmica/dynamics/sun_dynamics.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def get_sun_direction_eci(
    time: Annotated[
        np.datetime64 | npt.NDArray[np.datetime64],
        Doc("single element or array of UTC time values in np.datetime64 format."),
    ],
) -> Annotated[
    npt.NDArray[np.float64],
    Doc(
        "Corresponding Sun direction vector(s) (x,y,z) in eci."
        " If the input is a single element, the output is a 1D array."
        " If the input is an array, the output is a 2D array with the shape (len(time), 3).",
    ),
]:
    """Sun direction approximation function.

    The approximation considers an unperturbed motion of the Earth around the Sun with mean orbital elements
    that approximate the Sun's elliptic orbit wrt Earth around the year 2000.
    The algorithm is based on Oliver Montenbruck, Eberhard Grill; "Satellite Orbits: Models, Methods and Applications";
    Springer-Verlag Berlin. The alterations to the original model follow JAXA's Attitude Control Handbook
    JERG-2-510-HB001.
    """
    # tdt: Number of Julian Centuries. J2000 epoch -> 2451545.0 TT
    tdt = (juliandate(time) - 2451545.0) / 36525.0

    # m: mean anomaly
    m = np.deg2rad(357.5256 + 35999.045 * tdt)

    # right-ascension + argument of perigee = 282.94 deg
    # s: Sun's ecliptic longitude
    s = m + np.deg2rad(282.94 + (6892.0 * np.sin(m) + 72.0 * np.sin(2.0 * m) + 1250.09115 * tdt) / 3600.0 - 0.002652)

    # Obliquity of the ecliptic in radians
    epsilon_deg = 23.43929111
    epsilon = np.deg2rad(epsilon_deg)

    # Compute Sun Vector @ J2000. Shape: (3,) or (3, len(time))
    sun_vec_eci = np.array([np.cos(s), np.sin(s) * np.cos(epsilon), np.sin(s) * np.sin(epsilon)], dtype=np.float64)

    # normalize
    sun_vec_eci = sun_vec_eci.T  # shape: (3,) or (len(time), 3)

    # shape: (3,) or (len(time), 3)
    return normalize(sun_vec_eci, axis=-1)  # type: ignore[return-value]

make_satellite_orbit 

make_satellite_orbit(
    orbit_type: Literal["circular"],
    *args: Any,
    **kwargs: Any
) -> SatelliteOrbit
Source code in src/cosmica/dynamics/orbit.py 
274
275
276
277
278
279
280
281
282
283
def make_satellite_orbit(
    orbit_type: Literal["circular"],
    *args: Any,
    **kwargs: Any,
) -> SatelliteOrbit:
    if orbit_type == "circular":
        return CircularSatelliteOrbit(*args, **kwargs)
    else:
        msg = f"Unknown orbit type: {orbit_type}"
        raise ValueError(msg)