brahe 1.4.0

Brahe is a modern satellite dynamics library for research and engineering applications designed to be easy-to-learn, high-performance, and quick-to-deploy. The north-star of the development is enabling users to solve meaningful problems and answer questions quickly, easily, and correctly.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

/*!
Module providing drag force calculations for orbit dynamics.

This module implements the computation of atmospheric drag acceleration on spacecraft.
The drag force depends on atmospheric density, which is provided by separate atmospheric
density models.

For atmospheric density models, see [`crate::orbit_dynamics::atmospheric_density_models`].

# References

1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.83-86.
 */

use nalgebra::{Vector3, Vector6};

use crate::OMEGA_EARTH;
use crate::math::SMatrix3;

const OMEGA_VECTOR: Vector3<f64> = Vector3::new(0.0, 0.0, OMEGA_EARTH);

/// Computes the perturbing, non-conservative acceleration caused by atmospheric
/// drag assuming that the ballistic properties of the spacecraft are captured by
/// the coefficient of drag.
///
/// Arguments:
///
/// - `x_object`: Satellite Cartesean state in the inertial reference frame [m; m/s]
/// - `density`: atmospheric density [kg/m^3]
/// - `mass`: Spacecraft mass [kg]
/// - `area`: Wind-facing cross-sectional area [m^2]
/// - `drag_coefficient`: coefficient of drag [dimensionless]
/// - `T`: Rotation matrix from the inertial to the true-of-date frame
///
/// Return:
///
/// - `a`: Acceleration due to drag in the X, Y, and Z inertial directions. [m/s^2]
///
/// References:
///
/// 1. O. Montenbruck, and E. Gill, _Satellite Orbits: Models, Methods and Applications_, 2012, p.83-86.
///
#[allow(non_snake_case)]
pub fn accel_drag(
    x_object: Vector6<f64>,
    density: f64,
    mass: f64,
    area: f64,
    drag_coefficient: f64,
    T: SMatrix3,
) -> Vector3<f64> {
    // Position and velocity in true-of-date system
    let r_tod: Vector3<f64> = T * x_object.fixed_rows::<3>(0);
    let v_tod: Vector3<f64> = T * x_object.fixed_rows::<3>(3);

    // Velocity relative to the Earth's atmosphere
    let v_rel = v_tod - OMEGA_VECTOR.cross(&r_tod);
    let v_abs = v_rel.norm();

    // Acceleration
    let a_tod = -0.5 * drag_coefficient * (area / mass) * density * v_abs * v_rel;

    T.transpose() * a_tod
}

#[cfg(test)]
#[cfg_attr(coverage_nightly, coverage(off))]
#[allow(clippy::too_many_arguments)]
mod tests {
    use approx::assert_abs_diff_eq;
    use rstest::rstest;

    use crate::DEGREES;
    use crate::constants::R_EARTH;
    use crate::coordinates::*;
    use crate::orbit_dynamics::atmospheric_density_models::density_harris_priester;
    use crate::time::Epoch;
    use crate::{TimeSystem, sun_position};

    use super::*;

    #[test]
    fn test_accel_drag() {
        let oe = Vector6::new(R_EARTH + 500e3, 0.01, 97.3, 15.0, 30.0, 45.0);

        let x_object = state_koe_to_eci(oe, DEGREES);

        let a = accel_drag(x_object, 1.0e-12, 1000.0, 1.0, 2.0, SMatrix3::identity());

        assert_abs_diff_eq!(a.norm(), 5.97601877277239e-8, epsilon = 1.0e-10);
    }

    #[rstest]
    #[case(60310.000, 1.000, 100.000, 2.300, 4.884992303790e+06, 4.553508537449e+06, 1.330313604797e+06, - 5.300020580854e+03, 4.925709387424e+03, 2.601867699687e+03, 6.475291922546e-07, - 5.955905823290e-07, - 3.391292982024e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 3.854234369244e+06, 5.331427650549e+06, 1.762510444136e+06, - 6.281746044740e+03, 3.815301677691e+03, 2.254442602385e+03, 6.650957772889e-07, - 3.988840098278e-07, - 2.544421167373e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 2.670937897492e+06, 5.898362795150e+06, 2.124959710177e+06, - 7.013918036984e+03, 2.555322690457e+03, 1.818275122888e+03, 6.409667332941e-07, - 2.298120374578e-07, - 1.770183557869e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 1.382105394865e+06, 6.232286239143e+06, 2.403465989137e+06, - 7.467690875201e+03, 1.197415434669e+03, 1.311370423462e+03, 5.842585088334e-07, - 9.135822408135e-08, - 1.092477881395e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 3.879697748585e+04, 6.320698885147e+06, 2.587294936269e+06, - 7.626311735293e+03, - 2.030743556792e+02, 7.545803726769e+02, 5.026456239033e-07, 1.444392491547e-08, - 5.293306533864e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 1.306071545499e+06, 6.161054495134e+06, 2.669589120058e+06, - 7.485641457809e+03, - 1.589726432986e+03, 1.705831662309e+02, 4.081446766290e-07, 8.668767563667e-08, - 9.894677441492e-09)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 2.599961458555e+06, 5.760720193579e+06, 2.647597126838e+06, - 7.053988043184e+03, - 2.907616196131e+03, - 4.171887600156e+02, 3.131156928240e-07, 1.282887478988e-07, 1.969100563219e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 3.792857897424e+06, 5.136498372385e+06, 2.522713383918e+06, - 6.351318581368e+03, - 4.105664983317e+03, - 9.855989555158e+02, 2.275007278096e-07, 1.457511811933e-07, 3.751606537078e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 4.839229618329e+06, 4.313760582551e+06, 2.300338349966e+06, - 5.407960526895e+03, - 5.138647350964e+03, - 1.512772642576e+03, 1.559048510435e-07, 1.464885073622e-07, 4.630478305728e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 5.699723056095e+06, 3.325261742120e+06, 1.989578439892e+06, - 4.262928650995e+03, - 5.968779101105e+03, - 1.978931630018e+03, 9.954914915777e-08, 1.374999502478e-07, 4.899977063788e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 6.342536886568e+06, 2.209712299398e+06, 1.602811608208e+06, - 2.962019299787e+03, - 6.566856713480e+03, - 2.367070862012e+03, 5.686726682414e-08, 1.239386286371e-07, 4.805940901794e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 6.744447614457e+06, 1.010187496333e+06, 1.155147802510e+06, - 1.555801915035e+03, - 6.912958990231e+03, - 2.663466756131e+03, 2.527983956238e-08, 1.095212143544e-07, 4.542898162750e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 6.891477821536e+06, - 2.275518102869e+05, 6.638138965866e+05, - 9.761444639723e+01, - 6.996750915333e+03, - 2.858022017467e+03, 1.713172708983e-09, 9.741641119267e-08, 4.287172531695e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 6.779215094915e+06, - 1.456761075547e+06, 1.474909195700e+05, 1.358359721155e+03, - 6.817446444227e+03, - 2.944461654893e+03, - 1.773818124483e-08, 8.957550077534e-08, 4.171240578084e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 6.412800799786e+06, - 2.631381709006e+06, - 3.743717496543e+05, 2.758884958527e+03, - 6.383491905142e+03, - 2.920400175922e+03, - 3.788560161627e-08, 8.731043787305e-08, 4.310130707932e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 5.806614110708e+06, - 3.707634060799e+06, - 8.822442375822e+05, 4.053313697470e+03, - 5.712027156556e+03, - 2.787300862703e+03, - 6.446074751998e-08, 9.011680293095e-08, 4.749509376747e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 4.983679016998e+06, - 4.645498608912e+06, - 1.357188627116e+06, 5.195069678263e+03, - 4.828170587673e+03, - 2.550345567020e+03, - 1.036692911994e-07, 9.531054328296e-08, 5.444304022724e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 3.974821657117e+06, - 5.410057126297e+06, - 1.781501428157e+06, 6.142981046602e+03, - 3.764159581142e+03, - 2.218228747299e+03, - 1.614433532925e-07, 9.757449849799e-08, 6.229798340648e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 2.817603714143e+06, - 5.972669872748e+06, - 2.139313418925e+06, 6.862473412589e+03, - 2.558363041530e+03, - 1.802883666286e+03, - 2.421864827659e-07, 8.866437983603e-08, 6.793810291675e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 1.555056064850e+06, - 6.311966594700e+06, - 2.417135091026e+06, 7.326629008929e+03, - 1.254169479678e+03, - 1.319142854603e+03, - 3.470644698157e-07, 5.766114271947e-08, 6.667697200344e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, - 2.342365879764e+05, - 6.414628848619e+06, - 2.604335853094e+06, 7.517108189447e+03, 1.012550540356e+02, - 7.843300873158e+02, - 4.725061343330e-07, - 7.931859928099e-09, 5.257236492762e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, 1.095362092957e+06, - 6.275939733234e+06, - 2.693544355988e+06, 7.424915414211e+03, 1.458348575476e+03, - 2.177780628923e+02, - 6.065981611401e-07, - 1.200154294953e-07, 1.896063186108e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, 2.383524570841e+06, - 5.900075611853e+06, - 2.680956181964e+06, 7.050975739560e+03, 2.766680692859e+03, 3.597345651486e+02, - 7.304906266784e-07, - 2.860812404277e-07, - 3.969086848276e-08)]
    #[case(60310.000, 1.000, 100.000, 2.300, 3.581106522323e+06, - 5.300114080501e+06, - 2.566534189861e+06, 6.406472303857e+03, 3.976622102432e+03, 9.266306232258e+02, - 8.157042028116e-07, - 5.034459753379e-07, - 1.255579403572e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 4.641862770238e+06, - 4.497734593543e+06, - 2.354086603151e+06, 5.512883447615e+03, 5.041134556879e+03, 1.461264246671e+03, - 8.323328365445e-07, - 7.549157872248e-07, - 2.345768393719e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 5.524250386626e+06, - 3.522593991176e+06, - 2.051209062029e+06, 4.401653269803e+03, 5.917636197580e+03, 1.942721880224e+03, - 7.566927617697e-07, - 1.006810656914e-06, - 3.546733618769e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 6.193136243056e+06, - 2.411369562038e+06, - 1.669079863560e+06, 3.113434801059e+03, 6.569869917653e+03, 2.351671474407e+03, - 5.822225786598e-07, - 1.212620773785e-06, - 4.660943067555e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 6.621329023944e+06, - 1.206477537446e+06, - 1.222102840276e+06, 1.696862974319e+03, 6.969675450702e+03, 2.671223865848e+03, - 3.295192521274e-07, - 1.328583760690e-06, - 5.470991336177e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 6.790850714079e+06, 4.550542743298e+04, - 7.273998381722e+05, 2.068459741129e+02, 7.098544211312e+03, 2.887757995344e+03, - 4.232160931968e-08, - 1.329742026701e-06, - 5.815191957405e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 6.693864225551e+06, 1.295667609670e+06, - 2.041642113971e+05, - 1.297593298229e+03, 6.948825108710e+03, 2.991652446262e+03, 2.271126461274e-07, - 1.219593187243e-06, - 5.647373599657e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 6.333183868415e+06, 2.494761038735e+06, 3.271026349663e+05, - 2.755843531779e+03, 6.524454555312e+03, 2.977862301791e+03, 4.385586274848e-07, - 1.032983373613e-06, - 5.073840303605e-07)]
    #[case(60310.000, 1.000, 100.000, 2.300, 5.722314467058e+06, 3.595279736803e+06, 8.454390492652e+05, - 4.108466612970e+03, 5.841107621980e+03, 2.846284516417e+03, 5.767674716446e-07, - 8.133252603937e-07, - 4.268119157331e-07)]
    fn validate_acceleration_drag(
        #[case] mjd_tt: f64,
        #[case] area: f64,
        #[case] mass: f64,
        #[case] cd: f64,
        #[case] r_x: f64,
        #[case] r_y: f64,
        #[case] r_z: f64,
        #[case] v_x: f64,
        #[case] v_y: f64,
        #[case] v_z: f64,
        #[case] a_x_expected: f64,
        #[case] a_y_expected: f64,
        #[case] a_z_expected: f64,
    ) {
        let x = Vector6::new(r_x, r_y, r_z, v_x, v_y, v_z);

        let epc = Epoch::from_mjd(mjd_tt, TimeSystem::TT);
        let r_sun = sun_position(epc);

        let rho = density_harris_priester(Vector3::new(r_x, r_y, r_z), r_sun);

        let a = accel_drag(x, rho, mass, area, cd, SMatrix3::identity());

        let tol = 1.0e-12;
        assert_abs_diff_eq!(a[0], a_x_expected, epsilon = tol);
        assert_abs_diff_eq!(a[1], a_y_expected, epsilon = tol);
        assert_abs_diff_eq!(a[2], a_z_expected, epsilon = tol);
    }
}