Coverage for src/gncpy/orbital_mechanics.py: 0%
32 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-07-19 05:48 +0000
1"""Utility functions for performing orbital mechanics calculations."""
2import numpy as np
3from warnings import warn
5import gncpy.wgs84 as wgs84
8def ecc_anom_from_mean(mean_anom, ecc, tol=10**-7, max_iter=10**3):
9 """Calculate the eccentric anomaly from the mean anomaly.
11 Parameters
12 ----------
13 mean_anom : float
14 Mean anomaly in radians
15 ecc : float
16 eccentricity
17 tol : float, optional
18 Minimum tolerance for convergence. The default is 10**-7.
19 max_iter : int, optional
20 Maximum number of iterations to run. The default is 10**3.
22 Returns
23 -------
24 float
25 Eccentric anomaly in radians.
26 """
27 delta = 1 * 10**3
28 ecc_anom = mean_anom
29 ctr = 0
30 while np.abs(delta) > tol:
31 delta = (ecc_anom - ecc * np.sin(ecc_anom) - mean_anom) \
32 / (1 - ecc * np.cos(ecc_anom))
33 ecc_anom = ecc_anom - delta
35 ctr = ctr + 1
36 if ctr >= max_iter:
37 warn('Maximum iterations reached when finding eccentric anomaly.',
38 RuntimeWarning)
39 return ecc_anom
40 return ecc_anom
43def true_anom_from_ecc(ecc_anom, ecc):
44 """Calculate the true anomaly from the eccentric anomaly.
46 Parameters
47 ----------
48 ecc_anom : float
49 eccentric anomaly in radians
50 ecc : float
51 eccentricity
53 Returns
54 --------
55 float
56 true anomaly in radians
57 """
58 den = 1 / (1 - ecc * np.cos(ecc_anom))
59 cos_true = (np.cos(ecc_anom) - ecc) * den
60 sin_true = np.sqrt(1 - ecc**2) * np.sin(ecc_anom) * den
61 return np.arctan2(sin_true, cos_true)
64def correct_lon_ascend(lon_ascend, lon_rate, tk, toe):
65 """Correct the longitude of the ascending node.
67 Parameters
68 ----------
69 lon_ascend : float
70 Original longitude of the ascending node in radians
71 lon_rate : float
72 Rate of change of the longitude of the ascending node
73 tk : float
74 Current time of the week
75 toe : float
76 Time of ephemeris
78 Returns
79 -------
80 float
81 corrected longitude of the ascending node in radians
82 """
83 return lon_ascend + (lon_rate - wgs84.EARTH_ROT_RATE) * tk \
84 - wgs84.EARTH_ROT_RATE * toe
87def ecef_from_orbit(arg_lat, rad, lon_ascend, inc):
88 """Calculates the ECEF position from orbital parameters.
90 Parameters
91 ----------
92 arg_lat : float
93 Argument of latitude in radians
94 rad : float
95 Orbital radius in meters
96 lon_ascend : float
97 Longitude of the ascending node in radians
98 inc : float
99 Orbital inclination angle in radians
101 Returns
102 -------
103 3 x 1 numpy array
104 ECEF position in meters
105 """
106 xp = rad * np.cos(arg_lat)
107 yp = rad * np.sin(arg_lat)
109 c_lon = np.cos(lon_ascend)
110 s_lon = np.sin(lon_ascend)
111 c_inc = np.cos(inc)
113 x = xp * c_lon - yp * c_inc * s_lon
114 y = xp * s_lon + yp * c_inc * c_lon
115 z = yp * np.sin(inc)
117 return np.array([x, y, z]).reshape((3, 1))