nyx_space/od/ground_station/
trk_device.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use super::{ODAlmanacSnafu, ODError, ODTrajSnafu, TrackingDevice};
20use crate::md::prelude::{Interpolatable, Traj};
21use crate::od::msr::measurement::Measurement;
22use crate::od::msr::MeasurementType;
23use crate::time::Epoch;
24use crate::Spacecraft;
25use anise::errors::AlmanacResult;
26use anise::frames::Frame;
27use anise::prelude::{Almanac, Orbit};
28use hifitime::TimeUnits;
29use indexmap::IndexSet;
30use rand_pcg::Pcg64Mcg;
31use snafu::ResultExt;
32use std::sync::Arc;
33
34use super::GroundStation;
35
36impl TrackingDevice<Spacecraft> for GroundStation {
37    fn measurement_types(&self) -> &IndexSet<MeasurementType> {
38        &self.measurement_types
39    }
40
41    /// Perform a measurement from the ground station to the receiver (rx).
42    fn measure(
43        &mut self,
44        epoch: Epoch,
45        traj: &Traj<Spacecraft>,
46        rng: Option<&mut Pcg64Mcg>,
47        almanac: Arc<Almanac>,
48    ) -> Result<Option<Measurement>, ODError> {
49        match self.integration_time {
50            Some(integration_time) => {
51                // If out of traj bounds, return None, else the whole strand is rejected.
52                let rx_0 = match traj.at(epoch - integration_time) {
53                    Ok(rx) => rx,
54                    Err(_) => return Ok(None),
55                };
56
57                let rx_1 = match traj.at(epoch).context(ODTrajSnafu) {
58                    Ok(rx) => rx,
59                    Err(_) => return Ok(None),
60                };
61
62                let obstructing_body = if !self.frame.ephem_origin_match(rx_0.frame()) {
63                    Some(rx_0.frame())
64                } else {
65                    None
66                };
67
68                let aer_t0 = self
69                    .azimuth_elevation_of(rx_0.orbit, obstructing_body, &almanac)
70                    .context(ODAlmanacSnafu {
71                        action: "computing AER",
72                    })?;
73                let aer_t1 = self
74                    .azimuth_elevation_of(rx_1.orbit, obstructing_body, &almanac)
75                    .context(ODAlmanacSnafu {
76                        action: "computing AER",
77                    })?;
78
79                if aer_t0.elevation_deg < self.elevation_mask_deg
80                    || aer_t1.elevation_deg < self.elevation_mask_deg
81                {
82                    debug!(
83                        "{} (el. mask {:.3} deg) but object moves from {:.3} to {:.3} deg -- no measurement",
84                        self.name, self.elevation_mask_deg, aer_t0.elevation_deg, aer_t1.elevation_deg
85                    );
86                    return Ok(None);
87                } else if aer_t0.is_obstructed() || aer_t1.is_obstructed() {
88                    debug!(
89                        "{} obstruction at t0={}, t1={} -- no measurement",
90                        self.name,
91                        aer_t0.is_obstructed(),
92                        aer_t1.is_obstructed()
93                    );
94                    return Ok(None);
95                }
96
97                // Noises are computed at the midpoint of the integration time.
98                let noises = self.noises(epoch - integration_time * 0.5, rng)?;
99
100                let mut msr = Measurement::new(self.name.clone(), epoch + noises[0].seconds());
101
102                for (ii, msr_type) in self.measurement_types.iter().enumerate() {
103                    let msr_value = msr_type.compute_two_way(aer_t0, aer_t1, noises[ii + 1])?;
104                    msr.push(*msr_type, msr_value);
105                }
106
107                Ok(Some(msr))
108            }
109            None => self.measure_instantaneous(traj.at(epoch).context(ODTrajSnafu)?, rng, almanac),
110        }
111    }
112
113    fn name(&self) -> String {
114        self.name.clone()
115    }
116
117    fn location(&self, epoch: Epoch, frame: Frame, almanac: Arc<Almanac>) -> AlmanacResult<Orbit> {
118        almanac.transform_to(self.to_orbit(epoch, &almanac).unwrap(), frame, None)
119    }
120
121    fn measure_instantaneous(
122        &mut self,
123        rx: Spacecraft,
124        rng: Option<&mut Pcg64Mcg>,
125        almanac: Arc<Almanac>,
126    ) -> Result<Option<Measurement>, ODError> {
127        let obstructing_body = if !self.frame.ephem_origin_match(rx.frame()) {
128            Some(rx.frame())
129        } else {
130            None
131        };
132
133        let aer = self
134            .azimuth_elevation_of(rx.orbit, obstructing_body, &almanac)
135            .context(ODAlmanacSnafu {
136                action: "computing AER",
137            })?;
138
139        if aer.elevation_deg >= self.elevation_mask_deg && !aer.is_obstructed() {
140            // Only update the noises if the measurement is valid.
141            let noises = self.noises(rx.orbit.epoch, rng)?;
142
143            let mut msr = Measurement::new(self.name.clone(), rx.orbit.epoch + noises[0].seconds());
144
145            for (ii, msr_type) in self.measurement_types.iter().enumerate() {
146                let msr_value = msr_type.compute_one_way(aer, noises[ii + 1])?;
147                msr.push(*msr_type, msr_value);
148            }
149
150            Ok(Some(msr))
151        } else {
152            debug!(
153                "{} {} (el. mask {:.3} deg), object at {:.3} deg -- no measurement",
154                self.name, rx.orbit.epoch, self.elevation_mask_deg, aer.elevation_deg
155            );
156            Ok(None)
157        }
158    }
159
160    /// Returns the measurement noise of this ground station.
161    ///
162    /// # Methodology
163    /// Noises are modeled using a [StochasticNoise] process, defined by the sigma on the turn-on bias and on the steady state noise.
164    /// The measurement noise is computed assuming that all measurements are independent variables, i.e. the measurement matrix is
165    /// a diagonal matrix. The first item in the diagonal is the range noise (in km), set to the square of the steady state sigma. The
166    /// second item is the Doppler noise (in km/s), set to the square of the steady state sigma of that Gauss Markov process.
167    fn measurement_covar(&self, msr_type: MeasurementType, epoch: Epoch) -> Result<f64, ODError> {
168        let stochastics = self.stochastic_noises.as_ref().unwrap();
169
170        Ok(stochastics
171            .get(&msr_type)
172            .ok_or(ODError::NoiseNotConfigured {
173                kind: format!("{msr_type:?}"),
174            })?
175            .covariance(epoch))
176    }
177
178    fn measurement_bias(&self, msr_type: MeasurementType, _epoch: Epoch) -> Result<f64, ODError> {
179        let stochastics = self.stochastic_noises.as_ref().unwrap();
180
181        if let Some(gm) = stochastics
182            .get(&msr_type)
183            .ok_or(ODError::NoiseNotConfigured {
184                kind: format!("{msr_type:?}"),
185            })?
186            .bias
187        {
188            Ok(gm.constant.unwrap_or(0.0))
189        } else {
190            Ok(0.0)
191        }
192    }
193}
nyx_space/od/ground_station/trk_device.rs

nyx_space/od/ground_station/
trk_device.rs
