Coverage for src/gncpy/orbital_mechanics.py: 0%

1"""Utility functions for performing orbital mechanics calculations.""" 

2import numpy as np 

3from warnings import warn 

4 

5import gncpy.wgs84 as wgs84 

6 

7 

8def ecc_anom_from_mean(mean_anom, ecc, tol=10**-7, max_iter=10**3): 

9 """Calculate the eccentric anomaly from the mean anomaly. 

10 

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. 

21 

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 

34 

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 

41 

42 

43def true_anom_from_ecc(ecc_anom, ecc): 

44 """Calculate the true anomaly from the eccentric anomaly. 

45 

46 Parameters 

47 ---------- 

48 ecc_anom : float 

49 eccentric anomaly in radians 

50 ecc : float 

51 eccentricity 

52 

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) 

62 

63 

64def correct_lon_ascend(lon_ascend, lon_rate, tk, toe): 

65 """Correct the longitude of the ascending node. 

66 

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 

77 

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 

85 

86 

87def ecef_from_orbit(arg_lat, rad, lon_ascend, inc): 

88 """Calculates the ECEF position from orbital parameters. 

89 

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 

100 

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) 

108 

109 c_lon = np.cos(lon_ascend) 

110 s_lon = np.sin(lon_ascend) 

111 c_inc = np.cos(inc) 

112 

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) 

116 

117 return np.array([x, y, z]).reshape((3, 1))