nyx_space/od/ground_station/
mod.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 anise::astro::{Aberration, AzElRange, Location};
20use anise::errors::{AlmanacError, AlmanacResult};
21use anise::prelude::{Almanac, Frame, Orbit};
22use der::{Decode, Encode, Reader};
23use indexmap::{IndexMap, IndexSet};
24use snafu::ensure;
25
26use super::msr::MeasurementType;
27use super::noise::{GaussMarkov, StochasticNoise};
28use super::{ODAlmanacSnafu, ODError, ODTrajSnafu, TrackingDevice};
29use crate::io::ConfigRepr;
30use crate::od::NoiseNotConfiguredSnafu;
31use crate::time::Epoch;
32use hifitime::Duration;
33use rand_pcg::Pcg64Mcg;
34use serde::{Deserialize, Serialize};
35use std::fmt::{self, Debug};
36
37pub mod builtin;
38pub mod trk_device;
39
40#[cfg(feature = "python")]
41use pyo3::exceptions::PyValueError;
42#[cfg(feature = "python")]
43use pyo3::prelude::*;
44#[cfg(feature = "python")]
45use pyo3::types::{PyBytes, PyType};
46
47/// GroundStation defines a two-way ranging and doppler station.
48#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
49#[cfg_attr(feature = "python", pyclass)]
50pub struct GroundStation {
51    pub name: String,
52    pub location: Location,
53    pub measurement_types: IndexSet<MeasurementType>,
54    /// Duration needed to generate a measurement (if unset, it is assumed to be instantaneous)
55    pub integration_time: Option<Duration>,
56    /// Whether to correct for light travel time
57    pub light_time_correction: bool,
58    /// Noise on the timestamp of the measurement
59    pub timestamp_noise_s: Option<StochasticNoise>,
60    pub stochastic_noises: Option<IndexMap<MeasurementType, StochasticNoise>>,
61}
62
63#[cfg_attr(feature = "python", pymethods)]
64impl GroundStation {
65    /// Computes the azimuth and elevation of the provided object seen from this ground station, both in degrees.
66    /// This is a shortcut to almanac.azimuth_elevation_range_sez.
67    pub fn azimuth_elevation_of(
68        &self,
69        rx: Orbit,
70        obstructing_body: Option<Frame>,
71        almanac: &Almanac,
72    ) -> AlmanacResult<AzElRange> {
73        let ab_corr = if self.light_time_correction {
74            Aberration::LT
75        } else {
76            Aberration::NONE
77        };
78        almanac.azimuth_elevation_range_sez(
79            rx,
80            self.to_orbit(rx.epoch, almanac)?,
81            obstructing_body,
82            ab_corr,
83        )
84    }
85
86    /// Return this ground station as an orbit in its current frame
87    pub fn to_orbit(&self, epoch: Epoch, almanac: &Almanac) -> AlmanacResult<Orbit> {
88        Orbit::try_latlongalt(
89            self.location.latitude_deg,
90            self.location.longitude_deg,
91            self.location.height_km,
92            epoch,
93            almanac.frame_info(self.location.frame).map_err(|source| {
94                AlmanacError::GenericError {
95                    err: source.to_string(),
96                }
97            })?,
98        )
99        .map_err(|source| AlmanacError::AlmanacPhysics {
100            action: "building ground station location",
101            source: Box::new(source),
102        })
103    }
104}
105
106impl GroundStation {
107    /// Initializes a point on the surface of a celestial object.
108    /// This is meant for analysis, not for spacecraft navigation.
109    pub fn from_point(
110        name: String,
111        latitude_deg: f64,
112        longitude_deg: f64,
113        height_km: f64,
114        frame: Frame,
115    ) -> Self {
116        Self {
117            name,
118            location: Location {
119                latitude_deg,
120                longitude_deg,
121                height_km,
122                frame: frame.into(),
123                terrain_mask: vec![],
124                terrain_mask_ignored: true,
125            },
126            measurement_types: IndexSet::new(),
127            integration_time: None,
128            light_time_correction: false,
129            timestamp_noise_s: None,
130            stochastic_noises: None,
131        }
132    }
133
134    /// Returns a copy of this ground station with the new measurement type added (or replaced)
135    pub fn with_msr_type(mut self, msr_type: MeasurementType, noise: StochasticNoise) -> Self {
136        if self.stochastic_noises.is_none() {
137            self.stochastic_noises = Some(IndexMap::new());
138        }
139
140        self.stochastic_noises
141            .as_mut()
142            .unwrap()
143            .insert(msr_type, noise);
144
145        self.measurement_types.insert(msr_type);
146
147        self
148    }
149
150    /// Returns a copy of this ground station without the provided measurement type (if defined, else no error)
151    pub fn without_msr_type(mut self, msr_type: MeasurementType) -> Self {
152        if let Some(noises) = self.stochastic_noises.as_mut() {
153            noises.swap_remove(&msr_type);
154        }
155
156        self.measurement_types.swap_remove(&msr_type);
157
158        self
159    }
160
161    pub fn with_integration_time(mut self, integration_time: Option<Duration>) -> Self {
162        self.integration_time = integration_time;
163
164        self
165    }
166
167    /// Returns a copy of this ground station with the measurement type noises' constant bias set to the provided value.
168    pub fn with_msr_bias_constant(
169        mut self,
170        msr_type: MeasurementType,
171        bias_constant: f64,
172    ) -> Result<Self, ODError> {
173        if self.stochastic_noises.is_none() {
174            self.stochastic_noises = Some(IndexMap::new());
175        }
176
177        let stochastics = self.stochastic_noises.as_mut().unwrap();
178
179        let this_noise = stochastics
180            .get_mut(&msr_type)
181            .ok_or(ODError::NoiseNotConfigured {
182                kind: format!("{msr_type:?}"),
183            })
184            .unwrap();
185
186        if this_noise.bias.is_none() {
187            this_noise.bias = Some(GaussMarkov::ZERO);
188        }
189
190        this_noise.bias.unwrap().constant = Some(bias_constant);
191
192        Ok(self)
193    }
194
195    /// Returns the noises for all measurement types configured for this ground station at the provided epoch, timestamp noise is the first entry.
196    fn noises(&mut self, epoch: Epoch, rng: Option<&mut Pcg64Mcg>) -> Result<Vec<f64>, ODError> {
197        let mut noises = vec![0.0; self.measurement_types.len() + 1];
198
199        if let Some(rng) = rng {
200            ensure!(
201                self.stochastic_noises.is_some(),
202                NoiseNotConfiguredSnafu {
203                    kind: "ground station stochastics".to_string(),
204                }
205            );
206            // Add the timestamp noise first
207
208            if let Some(mut timestamp_noise) = self.timestamp_noise_s {
209                noises[0] = timestamp_noise.sample(epoch, rng);
210            }
211
212            let stochastics = self.stochastic_noises.as_mut().unwrap();
213
214            for (ii, msr_type) in self.measurement_types.iter().enumerate() {
215                noises[ii + 1] = stochastics
216                    .get_mut(msr_type)
217                    .ok_or(ODError::NoiseNotConfigured {
218                        kind: format!("{msr_type:?}"),
219                    })?
220                    .sample(epoch, rng);
221            }
222        }
223
224        Ok(noises)
225    }
226
227    fn available_data(&self) -> u8 {
228        let mut bits: u8 = 0;
229
230        if self.integration_time.is_some() {
231            bits |= 1 << 0;
232        }
233        if self.timestamp_noise_s.is_some() {
234            bits |= 1 << 1;
235        }
236        if self.stochastic_noises.is_some() {
237            bits |= 1 << 2;
238        }
239        bits
240    }
241}
242
243#[cfg(feature = "python")]
244#[cfg_attr(feature = "python", pymethods)]
245impl GroundStation {
246    /// Decodes an ASN.1 DER encoded byte array into a GroundStation object.
247    ///
248    /// :type data: bytes
249    /// :rtype: GroundStation
250    #[classmethod]
251    pub fn from_asn1(_cls: &Bound<'_, PyType>, data: &[u8]) -> PyResult<Self> {
252        match Self::from_der(data) {
253            Ok(obj) => Ok(obj),
254            Err(e) => Err(PyValueError::new_err(format!("ASN.1 decoding error: {e}"))),
255        }
256    }
257
258    /// Encodes this GroundStation object into an ASN.1 DER encoded byte array.
259    ///
260    /// :rtype: bytes
261    pub fn to_asn1<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
262        let mut buf = Vec::new();
263        match self.encode_to_vec(&mut buf) {
264            Ok(_) => Ok(PyBytes::new(py, &buf)),
265            Err(e) => Err(PyValueError::new_err(format!("ASN.1 encoding error: {e}"))),
266        }
267    }
268}
269
270impl Default for GroundStation {
271    fn default() -> Self {
272        let mut measurement_types = IndexSet::new();
273        measurement_types.insert(MeasurementType::Range);
274        measurement_types.insert(MeasurementType::Doppler);
275        Self {
276            name: "UNDEFINED".to_string(),
277            measurement_types,
278            location: Location::default(),
279            integration_time: None,
280            light_time_correction: false,
281            timestamp_noise_s: None,
282            stochastic_noises: None,
283        }
284    }
285}
286
287impl ConfigRepr for GroundStation {}
288
289#[derive(der::Sequence)]
290struct MsrNoisePair {
291    msr_type: MeasurementType,
292    noise: StochasticNoise,
293}
294
295impl<'a> Decode<'a> for GroundStation {
296    fn decode<R: Reader<'a>>(decoder: &mut R) -> der::Result<Self> {
297        let name: String = decoder.decode()?;
298        let location = decoder.decode()?;
299        // Measurement types are stored as a sequence of measurement types
300        let msr_types_vec: Vec<MeasurementType> = decoder.decode()?;
301        let measurement_types = IndexSet::from_iter(msr_types_vec);
302
303        let light_time_correction = decoder.decode()?;
304
305        // The flags tell us what happens next
306        let flags: u8 = decoder.decode()?;
307
308        let integration_time = if flags & (1 << 0) != 0 {
309            Some(Duration::from_total_nanoseconds(decoder.decode()?))
310        } else {
311            None
312        };
313
314        let timestamp_noise_s = if flags & (1 << 1) != 0 {
315            Some(decoder.decode()?)
316        } else {
317            None
318        };
319
320        let stochastic_noises = if flags & (1 << 2) != 0 {
321            // Stochastic noises are stored as a sequence of (MeasurementType, StochasticNoise) tuples (SEQUENCE of SEQUENCE)
322            // We define a helper struct for decoding
323
324            let stochastics_vec: Vec<MsrNoisePair> = decoder.decode()?;
325            let mut map = IndexMap::new();
326            for pair in stochastics_vec {
327                map.insert(pair.msr_type, pair.noise);
328            }
329            Some(map)
330        } else {
331            None
332        };
333
334        Ok(GroundStation {
335            name,
336            location,
337            measurement_types,
338            integration_time,
339            light_time_correction,
340            timestamp_noise_s,
341            stochastic_noises,
342        })
343    }
344}
345
346impl Encode for GroundStation {
347    fn encoded_len(&self) -> der::Result<der::Length> {
348        let msr_types_vec: Vec<MeasurementType> = self.measurement_types.iter().copied().collect();
349
350        let integration_time_ns = self.integration_time.map(|d| d.total_nanoseconds());
351
352        let stochastics_vec = self.stochastic_noises.as_ref().map(|map| {
353            map.iter()
354                .map(|(k, v)| MsrNoisePair {
355                    msr_type: *k,
356                    noise: *v,
357                })
358                .collect::<Vec<MsrNoisePair>>()
359        });
360
361        self.name.encoded_len()?
362            + self.location.encoded_len()?
363            + msr_types_vec.encoded_len()?
364            + self.light_time_correction.encoded_len()?
365            + self.available_data().encoded_len()?
366            + integration_time_ns.encoded_len()?
367            + self.timestamp_noise_s.encoded_len()?
368            + stochastics_vec.encoded_len()?
369    }
370
371    fn encode(&self, encoder: &mut impl der::Writer) -> der::Result<()> {
372        self.name.encode(encoder)?;
373        self.location.encode(encoder)?;
374
375        let msr_types_vec: Vec<MeasurementType> = self.measurement_types.iter().copied().collect();
376        msr_types_vec.encode(encoder)?;
377
378        self.light_time_correction.encode(encoder)?;
379        self.available_data().encode(encoder)?;
380
381        let integration_time_ns = self.integration_time.map(|d| d.total_nanoseconds());
382        integration_time_ns.encode(encoder)?;
383
384        self.timestamp_noise_s.encode(encoder)?;
385
386        let stochastics_vec = self.stochastic_noises.as_ref().map(|map| {
387            map.iter()
388                .map(|(k, v)| MsrNoisePair {
389                    msr_type: *k,
390                    noise: *v,
391                })
392                .collect::<Vec<MsrNoisePair>>()
393        });
394        stochastics_vec.encode(encoder)?;
395
396        Ok(())
397    }
398}
399
400impl fmt::Display for GroundStation {
401    // Prints the Keplerian orbital elements with units
402    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
403        write!(f, "{} ({})", self.name, self.location,)
404    }
405}
406
407#[cfg(test)]
408mod gs_ut {
409
410    use anise::astro::{Location, TerrainMask};
411    use anise::constants::frames::IAU_EARTH_FRAME;
412    use indexmap::{IndexMap, IndexSet};
413
414    use crate::io::ConfigRepr;
415    use crate::od::prelude::*;
416
417    #[test]
418    fn test_load_single() {
419        use std::env;
420        use std::path::PathBuf;
421
422        use hifitime::TimeUnits;
423
424        let test_data: PathBuf = [
425            env::var("CARGO_MANIFEST_DIR").unwrap(),
426            "data".to_string(),
427            "03_tests".to_string(),
428            "config".to_string(),
429            "one_ground_station.yaml".to_string(),
430        ]
431        .iter()
432        .collect();
433
434        assert!(test_data.exists(), "Could not find the test data");
435
436        let gs = GroundStation::load(test_data).unwrap();
437
438        dbg!(&gs);
439
440        let mut measurement_types = IndexSet::new();
441        measurement_types.insert(MeasurementType::Range);
442        measurement_types.insert(MeasurementType::Doppler);
443
444        let mut stochastics = IndexMap::new();
445        stochastics.insert(
446            MeasurementType::Range,
447            StochasticNoise {
448                bias: Some(GaussMarkov::new(1.days(), 5e-3).unwrap()),
449                ..Default::default()
450            },
451        );
452        stochastics.insert(
453            MeasurementType::Doppler,
454            StochasticNoise {
455                bias: Some(GaussMarkov::new(1.days(), 5e-5).unwrap()),
456                ..Default::default()
457            },
458        );
459
460        let expected_gs = GroundStation {
461            name: "Demo ground station".to_string(),
462            location: Location {
463                latitude_deg: 2.3522,
464                longitude_deg: 48.8566,
465                height_km: 0.4,
466                frame: IAU_EARTH_FRAME.into(),
467                terrain_mask: TerrainMask::from_flat_terrain(5.0),
468                terrain_mask_ignored: false,
469            },
470            measurement_types,
471            stochastic_noises: Some(stochastics),
472
473            light_time_correction: false,
474            timestamp_noise_s: None,
475            integration_time: Some(60 * Unit::Second),
476        };
477
478        println!("{}", serde_yml::to_string(&expected_gs).unwrap());
479
480        assert_eq!(expected_gs, gs);
481    }
482
483    #[test]
484    fn test_load_many() {
485        use hifitime::TimeUnits;
486        use std::env;
487        use std::path::PathBuf;
488
489        let test_file: PathBuf = [
490            env::var("CARGO_MANIFEST_DIR").unwrap(),
491            "data".to_string(),
492            "03_tests".to_string(),
493            "config".to_string(),
494            "many_ground_stations.yaml".to_string(),
495        ]
496        .iter()
497        .collect();
498
499        let stations = GroundStation::load_many(test_file).unwrap();
500
501        dbg!(&stations);
502
503        let mut measurement_types = IndexSet::new();
504        measurement_types.insert(MeasurementType::Range);
505        measurement_types.insert(MeasurementType::Doppler);
506
507        let mut stochastics = IndexMap::new();
508        stochastics.insert(
509            MeasurementType::Range,
510            StochasticNoise {
511                bias: Some(GaussMarkov::new(1.days(), 5e-3).unwrap()),
512                ..Default::default()
513            },
514        );
515        stochastics.insert(
516            MeasurementType::Doppler,
517            StochasticNoise {
518                bias: Some(GaussMarkov::new(1.days(), 5e-5).unwrap()),
519                ..Default::default()
520            },
521        );
522
523        let expected = vec![
524            GroundStation {
525                name: "Demo ground station".to_string(),
526                location: Location {
527                    latitude_deg: 2.3522,
528                    longitude_deg: 48.8566,
529                    height_km: 0.4,
530                    frame: IAU_EARTH_FRAME.into(),
531                    terrain_mask: TerrainMask::from_flat_terrain(5.0),
532                    terrain_mask_ignored: false,
533                },
534                measurement_types: measurement_types.clone(),
535                stochastic_noises: Some(stochastics.clone()),
536                light_time_correction: false,
537                timestamp_noise_s: None,
538                integration_time: None,
539            },
540            GroundStation {
541                name: "Canberra".to_string(),
542                location: Location {
543                    latitude_deg: -35.398333,
544                    longitude_deg: 148.981944,
545                    height_km: 0.691750,
546                    frame: IAU_EARTH_FRAME.into(),
547                    terrain_mask: TerrainMask::from_flat_terrain(5.0),
548                    terrain_mask_ignored: false,
549                },
550                measurement_types,
551                stochastic_noises: Some(stochastics),
552                light_time_correction: false,
553                timestamp_noise_s: None,
554                integration_time: None,
555            },
556        ];
557
558        assert_eq!(expected, stations);
559
560        // Serialize back
561        let reser = serde_yml::to_string(&expected).unwrap();
562        dbg!(reser);
563    }
564}
nyx_space/od/ground_station/mod.rs

nyx_space/od/ground_station/
mod.rs
