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

nyx_space/od/ground_station/
mod.rs
