nyx_space/md/trajectory/
sc_traj.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;
20use anise::constants::orientations::J2000;
21use anise::errors::AlmanacError;
22use anise::prelude::{Almanac, Frame, Orbit};
23use arrow::array::RecordBatchReader;
24use arrow::array::{Float64Array, StringArray};
25use hifitime::TimeSeries;
26use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
27use snafu::{ensure, ResultExt};
28
29use super::TrajError;
30use super::{ExportCfg, Traj};
31use crate::cosmic::Spacecraft;
32use crate::errors::{FromAlmanacSnafu, NyxError};
33use crate::io::watermark::prj_name_ver;
34use crate::io::{InputOutputError, MissingDataSnafu, ParquetSnafu, StdIOSnafu};
35use crate::md::prelude::{Interpolatable, StateParameter};
36use crate::md::EventEvaluator;
37use crate::time::{Duration, Epoch, Format, Formatter, TimeUnits};
38use crate::State;
39use std::collections::{HashMap, HashSet};
40use std::error::Error;
41use std::fs::File;
42use std::io::{BufRead, BufReader, BufWriter, Write};
43use std::path::{Path, PathBuf};
44use std::str::FromStr;
45use std::sync::Arc;
46#[cfg(not(target_arch = "wasm32"))]
47use std::time::Instant;
48
49impl Traj<Spacecraft> {
50    /// Builds a new trajectory built from the SPICE BSP (SPK) file loaded in the provided Almanac, provided the start and stop epochs.
51    ///
52    /// If the start and stop epochs are not provided, then the full domain of the trajectory will be used.
53    pub fn from_bsp(
54        target_frame: Frame,
55        observer_frame: Frame,
56        almanac: Arc<Almanac>,
57        sc_template: Spacecraft,
58        step: Duration,
59        start_epoch: Option<Epoch>,
60        end_epoch: Option<Epoch>,
61        ab_corr: Option<Aberration>,
62        name: Option<String>,
63    ) -> Result<Self, AlmanacError> {
64        let (domain_start, domain_end) =
65            almanac
66                .spk_domain(target_frame.ephemeris_id)
67                .map_err(|e| AlmanacError::Ephemeris {
68                    action: "could not fetch domain",
69                    source: Box::new(e),
70                })?;
71
72        let start_epoch = start_epoch.unwrap_or(domain_start);
73        let end_epoch = end_epoch.unwrap_or(domain_end);
74
75        let time_series = TimeSeries::inclusive(start_epoch, end_epoch, step);
76        let mut states = Vec::with_capacity(time_series.len());
77        for epoch in time_series {
78            let orbit = almanac.transform(target_frame, observer_frame, epoch, ab_corr)?;
79
80            states.push(sc_template.with_orbit(orbit));
81        }
82
83        Ok(Self { name, states })
84    }
85    /// Allows converting the source trajectory into the (almost) equivalent trajectory in another frame
86    #[allow(clippy::map_clone)]
87    pub fn to_frame(&self, new_frame: Frame, almanac: Arc<Almanac>) -> Result<Self, NyxError> {
88        if self.states.is_empty() {
89            return Err(NyxError::Trajectory {
90                source: TrajError::CreationError {
91                    msg: "No trajectory to convert".to_string(),
92                },
93            });
94        }
95
96        #[cfg(not(target_arch = "wasm32"))]
97        let start_instant = Instant::now();
98        let mut traj = Self::new();
99        for state in &self.states {
100            let new_orbit =
101                almanac
102                    .transform_to(state.orbit, new_frame, None)
103                    .context(FromAlmanacSnafu {
104                        action: "transforming trajectory into new frame",
105                    })?;
106            traj.states.push(state.with_orbit(new_orbit));
107        }
108        traj.finalize();
109
110        #[cfg(not(target_arch = "wasm32"))]
111        info!(
112            "Converted trajectory from {} to {} in {} ms: {traj}",
113            self.first().orbit.frame,
114            new_frame,
115            (Instant::now() - start_instant).as_millis()
116        );
117
118        #[cfg(target_arch = "wasm32")]
119        info!(
120            "Converted trajectory from {} to {}: {traj}",
121            self.first().orbit.frame,
122            new_frame,
123        );
124
125        Ok(traj)
126    }
127
128    /// Exports this trajectory to the provided filename in parquet format with only the epoch, the geodetic latitude, longitude, and height at one state per minute.
129    /// Must provide a body fixed frame to correctly compute the latitude and longitude.
130    #[allow(clippy::identity_op)]
131    pub fn to_groundtrack_parquet<P: AsRef<Path>>(
132        &self,
133        path: P,
134        body_fixed_frame: Frame,
135        events: Option<Vec<&dyn EventEvaluator<Spacecraft>>>,
136        metadata: Option<HashMap<String, String>>,
137        almanac: Arc<Almanac>,
138    ) -> Result<PathBuf, Box<dyn Error>> {
139        let traj = self.to_frame(body_fixed_frame, almanac.clone())?;
140
141        let mut cfg = ExportCfg::builder()
142            .step(1.minutes())
143            .fields(vec![
144                StateParameter::Latitude,
145                StateParameter::Longitude,
146                StateParameter::Height,
147                StateParameter::Rmag,
148            ])
149            .build();
150        cfg.metadata = metadata;
151
152        traj.to_parquet(path, events, cfg, almanac)
153    }
154
155    /// Initialize a new spacecraft trajectory from the path to a CCSDS OEM file.
156    ///
157    /// CCSDS OEM only contains the orbit information but Nyx builds spacecraft trajectories.
158    /// If not spacecraft template is provided, then a default massless spacecraft will be built.
159    pub fn from_oem_file<P: AsRef<Path>>(
160        path: P,
161        tpl_option: Option<Spacecraft>,
162    ) -> Result<Self, NyxError> {
163        // Open the file
164        let file = File::open(path).map_err(|e| NyxError::CCSDS {
165            msg: format!("File opening error: {e}"),
166        })?;
167        let reader = BufReader::new(file);
168
169        let template = tpl_option.unwrap_or_default();
170
171        // Parse the Orbit Element messages
172        let mut time_system = String::new();
173
174        let ignored_tokens: HashSet<_> = [
175            "CCSDS_OMM_VERS".to_string(),
176            "CREATION_DATE".to_string(),
177            "ORIGINATOR".to_string(),
178        ]
179        .into();
180
181        let mut traj = Self::default();
182
183        let mut parse = false;
184
185        let mut center_name = None;
186        let mut orient_name = None;
187
188        'lines: for (lno, line) in reader.lines().enumerate() {
189            let line = line.map_err(|e| NyxError::CCSDS {
190                msg: format!("File read error: {e}"),
191            })?;
192            let line = line.trim();
193            if line.is_empty() {
194                continue;
195            }
196
197            if ignored_tokens.iter().any(|t| line.starts_with(t)) {
198                continue 'lines;
199            }
200            if line.starts_with("OBJECT_NAME") {
201                // Extract the object ID from the line
202                let parts: Vec<&str> = line.split('=').collect();
203                let name = parts[1].trim().to_string();
204                debug!("[line: {}] Found object {name}", lno + 1);
205                traj.name = Some(name);
206            } else if line.starts_with("CENTER_NAME") {
207                let parts: Vec<&str> = line.split('=').collect();
208                center_name = Some(parts[1].trim().to_owned());
209            } else if line.starts_with("REF_FRAME") {
210                let parts: Vec<&str> = line.split('=').collect();
211                orient_name = Some(parts[1].trim().to_owned());
212            } else if line.starts_with("TIME_SYSTEM") {
213                let parts: Vec<&str> = line.split('=').collect();
214                time_system = parts[1].trim().to_string();
215                debug!("[line: {}] Found time system `{time_system}`", lno + 1);
216            } else if line.starts_with("META_STOP") {
217                // We can start parsing now
218                parse = true;
219            } else if line.starts_with("META_START") {
220                // Stop the parsing
221                parse = false;
222            } else if line.starts_with("COVARIANCE_START") {
223                // Stop the parsing
224                warn!("[line: {}] Skipping covariance in OEM parsing", lno + 1);
225                parse = false;
226            } else if parse {
227                let frame = Frame::from_name(
228                    center_name.clone().unwrap().as_str(),
229                    orient_name.clone().unwrap().as_str(),
230                )
231                .map_err(|e| NyxError::CCSDS {
232                    msg: format!("frame error `{center_name:?} {orient_name:?}`: {e}"),
233                })?;
234                // Split the line into components
235                let parts: Vec<&str> = line.split_whitespace().collect();
236
237                if parts.len() < 7 {
238                    debug!("[line: {}] Could not understand `{parts:?}`", lno + 1);
239                } else {
240                    // Extract the values
241                    let epoch_str = format!("{} {time_system}", parts[0]);
242                    match parts[1].parse::<f64>() {
243                        Ok(x_km) => {
244                            // Look good!
245                            let y_km = parts[2].parse::<f64>().unwrap();
246                            let z_km = parts[3].parse::<f64>().unwrap();
247                            let vx_km_s = parts[4].parse::<f64>().unwrap();
248                            let vy_km_s = parts[5].parse::<f64>().unwrap();
249                            let vz_km_s = parts[6].parse::<f64>().unwrap();
250
251                            let epoch =
252                                Epoch::from_str(epoch_str.trim()).map_err(|e| NyxError::CCSDS {
253                                    msg: format!("Parsing epoch error: {e}"),
254                                })?;
255
256                            let orbit = Orbit::new(
257                                x_km, y_km, z_km, vx_km_s, vy_km_s, vz_km_s, epoch, frame,
258                            );
259
260                            traj.states.push(template.with_orbit(orbit));
261                        }
262                        Err(_) => {
263                            // Probably a comment
264                            debug!("[line: {}] Could not parse `{parts:?}`", lno + 1);
265                            continue;
266                        }
267                    };
268                }
269            }
270        }
271
272        traj.finalize();
273
274        Ok(traj)
275    }
276
277    pub fn to_oem_file<P: AsRef<Path>>(
278        &self,
279        path: P,
280        cfg: ExportCfg,
281    ) -> Result<PathBuf, NyxError> {
282        if self.states.is_empty() {
283            return Err(NyxError::CCSDS {
284                msg: "Cannot export an empty trajectory to OEM".to_string(),
285            });
286        }
287        let tick = Epoch::now().unwrap();
288        info!("Exporting trajectory to CCSDS OEM file...");
289
290        // Grab the path here before we move stuff.
291        let path_buf = cfg.actual_path(path);
292
293        let metadata = cfg.metadata.unwrap_or_default();
294
295        let file = File::create(&path_buf).map_err(|e| NyxError::CCSDS {
296            msg: format!("File creation error: {e}"),
297        })?;
298        let mut writer = BufWriter::new(file);
299
300        let err_hdlr = |e| NyxError::CCSDS {
301            msg: format!("Could not write: {e}"),
302        };
303
304        // Build the states iterator -- this does require copying the current states but I can't either get a reference or a copy of all the states.
305        let states = if cfg.start_epoch.is_some() || cfg.end_epoch.is_some() || cfg.step.is_some() {
306            // Must interpolate the data!
307            let start = cfg.start_epoch.unwrap_or_else(|| self.first().epoch());
308            let end = cfg.end_epoch.unwrap_or_else(|| self.last().epoch());
309            let step = cfg.step.unwrap_or_else(|| 1.minutes());
310            self.every_between(step, start, end).collect()
311        } else {
312            self.states.to_vec()
313        };
314
315        // Epoch formmatter.
316        let iso8601_no_ts = Format::from_str("%Y-%m-%dT%H:%M:%S.%f").unwrap();
317
318        // Write mandatory metadata
319        writeln!(writer, "CCSDS_OMM_VERS = 2.0").map_err(err_hdlr)?;
320
321        writeln!(
322            writer,
323            "COMMENT Built by {} -- https://nyxspace.com/\n",
324            prj_name_ver()
325        )
326        .map_err(err_hdlr)?;
327        writeln!(
328            writer,
329            "COMMENT Nyx Space provided under the AGPL v3 open source license -- https://nyxspace.com/pricing\n"
330        )
331        .map_err(err_hdlr)?;
332
333        writeln!(
334            writer,
335            "CREATION_DATE = {}",
336            Formatter::new(Epoch::now().unwrap(), iso8601_no_ts)
337        )
338        .map_err(err_hdlr)?;
339        writeln!(
340            writer,
341            "ORIGINATOR = {}\n",
342            metadata
343                .get("originator")
344                .unwrap_or(&"Nyx Space".to_string())
345        )
346        .map_err(err_hdlr)?;
347
348        writeln!(writer, "META_START").map_err(err_hdlr)?;
349        // Write optional metadata
350        if let Some(object_name) = metadata.get("object_name") {
351            writeln!(writer, "\tOBJECT_NAME = {}", object_name).map_err(err_hdlr)?;
352        } else if let Some(object_name) = &self.name {
353            writeln!(writer, "\tOBJECT_NAME = {}", object_name).map_err(err_hdlr)?;
354        }
355
356        let first_orbit = states[0].orbit;
357        let first_frame = first_orbit.frame;
358        let frame_str = format!(
359            "{first_frame:e} {}",
360            match first_frame.orientation_id {
361                J2000 => "ICRF".to_string(),
362                _ => format!("{first_frame:o}"),
363            }
364        );
365        let splt: Vec<&str> = frame_str.split(' ').collect();
366        let center = splt[0];
367        let ref_frame = frame_str.replace(center, " ");
368        writeln!(
369            writer,
370            "\tREF_FRAME = {}",
371            match ref_frame.trim() {
372                "J2000" => "ICRF",
373                _ => ref_frame.trim(),
374            }
375        )
376        .map_err(err_hdlr)?;
377
378        writeln!(writer, "\tCENTER_NAME = {center}",).map_err(err_hdlr)?;
379
380        writeln!(writer, "\tTIME_SYSTEM = {}", first_orbit.epoch.time_scale).map_err(err_hdlr)?;
381
382        writeln!(
383            writer,
384            "\tSTART_TIME = {}",
385            Formatter::new(states[0].epoch(), iso8601_no_ts)
386        )
387        .map_err(err_hdlr)?;
388        writeln!(
389            writer,
390            "\tUSEABLE_START_TIME = {}",
391            Formatter::new(states[0].epoch(), iso8601_no_ts)
392        )
393        .map_err(err_hdlr)?;
394        writeln!(
395            writer,
396            "\tUSEABLE_STOP_TIME = {}",
397            Formatter::new(states[states.len() - 1].epoch(), iso8601_no_ts)
398        )
399        .map_err(err_hdlr)?;
400        writeln!(
401            writer,
402            "\tSTOP_TIME = {}",
403            Formatter::new(states[states.len() - 1].epoch(), iso8601_no_ts)
404        )
405        .map_err(err_hdlr)?;
406
407        writeln!(writer, "META_STOP\n").map_err(err_hdlr)?;
408
409        for sc_state in &states {
410            let state = sc_state.orbit;
411            writeln!(
412                writer,
413                "{} {:E} {:E} {:E} {:E} {:E} {:E}",
414                Formatter::new(state.epoch, iso8601_no_ts),
415                state.radius_km.x,
416                state.radius_km.y,
417                state.radius_km.z,
418                state.velocity_km_s.x,
419                state.velocity_km_s.y,
420                state.velocity_km_s.z
421            )
422            .map_err(err_hdlr)?;
423        }
424
425        #[allow(clippy::writeln_empty_string)]
426        writeln!(writer, "").map_err(err_hdlr)?;
427
428        // Return the path this was written to
429        let tock_time = Epoch::now().unwrap() - tick;
430        info!(
431            "Trajectory written to {} in {tock_time}",
432            path_buf.display()
433        );
434        Ok(path_buf)
435    }
436
437    pub fn from_parquet<P: AsRef<Path>>(path: P) -> Result<Self, InputOutputError> {
438        let file = File::open(&path).context(StdIOSnafu {
439            action: "opening trajectory file",
440        })?;
441
442        let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
443
444        let mut metadata = HashMap::new();
445        // Build the custom metadata
446        if let Some(file_metadata) = builder.metadata().file_metadata().key_value_metadata() {
447            for key_value in file_metadata {
448                if !key_value.key.starts_with("ARROW:") {
449                    metadata.insert(
450                        key_value.key.clone(),
451                        key_value.value.clone().unwrap_or("[unset]".to_string()),
452                    );
453                }
454            }
455        }
456
457        // Check the schema
458        let mut has_epoch = false; // Required
459        let mut frame = None;
460
461        let mut found_fields = vec![
462            (StateParameter::X, false),
463            (StateParameter::Y, false),
464            (StateParameter::Z, false),
465            (StateParameter::VX, false),
466            (StateParameter::VY, false),
467            (StateParameter::VZ, false),
468            (StateParameter::PropMass, false),
469        ];
470
471        let reader = builder.build().context(ParquetSnafu {
472            action: "building output trajectory file",
473        })?;
474
475        for field in &reader.schema().fields {
476            if field.name().as_str() == "Epoch (UTC)" {
477                has_epoch = true;
478            } else {
479                for potential_field in &mut found_fields {
480                    if field.name() == potential_field.0.to_field(None).name() {
481                        potential_field.1 = true;
482                        if potential_field.0 != StateParameter::PropMass {
483                            if let Some(frame_info) = field.metadata().get("Frame") {
484                                // Frame is expected to be serialized as Dhall.
485                                match serde_dhall::from_str(frame_info).parse::<Frame>() {
486                                    Err(e) => {
487                                        return Err(InputOutputError::ParseDhall {
488                                            data: frame_info.to_string(),
489                                            err: format!("{e}"),
490                                        })
491                                    }
492                                    Ok(deser_frame) => frame = Some(deser_frame),
493                                };
494                            }
495                        }
496                        break;
497                    }
498                }
499            }
500        }
501
502        ensure!(
503            has_epoch,
504            MissingDataSnafu {
505                which: "Epoch (UTC)"
506            }
507        );
508
509        ensure!(
510            frame.is_some(),
511            MissingDataSnafu {
512                which: "Frame in metadata"
513            }
514        );
515
516        for (field, exists) in found_fields.iter().take(found_fields.len() - 1) {
517            ensure!(
518                exists,
519                MissingDataSnafu {
520                    which: format!("Missing `{}` field", field.to_field(None).name())
521                }
522            );
523        }
524
525        let sc_compat = found_fields.last().unwrap().1;
526
527        let expected_type = std::any::type_name::<Spacecraft>()
528            .split("::")
529            .last()
530            .unwrap();
531
532        if expected_type == "Spacecraft" {
533            ensure!(
534                sc_compat,
535                MissingDataSnafu {
536                    which: format!(
537                        "Missing `{}` field",
538                        found_fields.last().unwrap().0.to_field(None).name()
539                    )
540                }
541            );
542        } else if sc_compat {
543            // Not a spacecraft, remove the prop mass
544            if let Some(last_field) = found_fields.last_mut() {
545                if last_field.0 == StateParameter::PropMass && last_field.1 {
546                    last_field.1 = false;
547                }
548            }
549        }
550
551        // At this stage, we know that the measurement is valid and the conversion is supported.
552        let mut traj = Traj::default();
553
554        // Now convert each batch on the fly
555        for maybe_batch in reader {
556            let batch = maybe_batch.unwrap();
557
558            let epochs = batch
559                .column_by_name("Epoch (UTC)")
560                .unwrap()
561                .as_any()
562                .downcast_ref::<StringArray>()
563                .unwrap();
564
565            let mut shared_data = vec![];
566
567            for (field, _) in found_fields.iter().take(found_fields.len() - 1) {
568                shared_data.push(
569                    batch
570                        .column_by_name(field.to_field(None).name())
571                        .unwrap()
572                        .as_any()
573                        .downcast_ref::<Float64Array>()
574                        .unwrap(),
575                );
576            }
577
578            if expected_type == "Spacecraft" {
579                // Read the prop only if this is a spacecraft we're building
580                shared_data.push(
581                    batch
582                        .column_by_name("prop_mass (kg)")
583                        .unwrap()
584                        .as_any()
585                        .downcast_ref::<Float64Array>()
586                        .unwrap(),
587                );
588            }
589
590            // Grab the frame -- it should have been serialized with all of the data so we don't need to reload it.
591
592            // Build the states
593            for i in 0..batch.num_rows() {
594                let mut state = Spacecraft::zeros();
595                state.set_epoch(Epoch::from_gregorian_str(epochs.value(i)).map_err(|e| {
596                    InputOutputError::Inconsistency {
597                        msg: format!("{e} when parsing epoch"),
598                    }
599                })?);
600                state.set_frame(frame.unwrap()); // We checked it was set above with an ensure! call
601                state.unset_stm(); // We don't have any STM data, so let's unset this.
602
603                for (j, (param, exists)) in found_fields.iter().enumerate() {
604                    if *exists {
605                        state.set_value(*param, shared_data[j].value(i)).unwrap();
606                    }
607                }
608
609                traj.states.push(state);
610            }
611        }
612
613        // Remove any duplicates that may exist in the imported trajectory.
614        traj.finalize();
615
616        Ok(traj)
617    }
618}
619
620#[cfg(test)]
621mod ut_ccsds_oem {
622
623    use crate::md::prelude::{OrbitalDynamics, Propagator, SpacecraftDynamics};
624    use crate::time::{Epoch, TimeUnits};
625    use crate::Spacecraft;
626    use crate::{io::ExportCfg, md::prelude::Traj, Orbit};
627    use anise::almanac::Almanac;
628    use anise::constants::frames::MOON_J2000;
629    use pretty_env_logger;
630    use std::env;
631    use std::str::FromStr;
632    use std::sync::Arc;
633    use std::{collections::HashMap, path::PathBuf};
634
635    #[test]
636    fn test_load_oem_leo() {
637        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
638        let path: PathBuf = [
639            env!("CARGO_MANIFEST_DIR"),
640            "data",
641            "03_tests",
642            "ccsds",
643            "oem",
644            "LEO_10s.oem",
645        ]
646        .iter()
647        .collect();
648
649        let _ = pretty_env_logger::try_init();
650
651        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
652
653        // This trajectory has two duplicate epochs, which should be removed by the call to finalize()
654        assert_eq!(traj.states.len(), 361);
655        assert_eq!(traj.name.unwrap(), "TEST_OBJ".to_string());
656    }
657
658    #[test]
659    fn test_load_oem_meo() {
660        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
661        let path: PathBuf = [
662            env!("CARGO_MANIFEST_DIR"),
663            "data",
664            "03_tests",
665            "ccsds",
666            "oem",
667            "MEO_60s.oem",
668        ]
669        .iter()
670        .collect();
671
672        let _ = pretty_env_logger::try_init();
673
674        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
675
676        assert_eq!(traj.states.len(), 61);
677        assert_eq!(traj.name.unwrap(), "TEST_OBJ".to_string());
678    }
679
680    #[test]
681    fn test_load_oem_geo() {
682        use pretty_env_logger;
683        use std::env;
684
685        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
686        let path: PathBuf = [
687            env!("CARGO_MANIFEST_DIR"),
688            "data",
689            "03_tests",
690            "ccsds",
691            "oem",
692            "GEO_20s.oem",
693        ]
694        .iter()
695        .collect();
696
697        let _ = pretty_env_logger::try_init();
698
699        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
700
701        assert_eq!(traj.states.len(), 181);
702        assert_eq!(traj.name.as_ref().unwrap(), &"TEST_OBJ".to_string());
703
704        // Reexport this to CCSDS.
705        let cfg = ExportCfg::builder()
706            .timestamp(true)
707            .metadata(HashMap::from([
708                ("originator".to_string(), "Test suite".to_string()),
709                ("object_name".to_string(), "TEST_OBJ".to_string()),
710            ]))
711            .build();
712
713        let path: PathBuf = [
714            env!("CARGO_MANIFEST_DIR"),
715            "data",
716            "04_output",
717            "GEO_20s_rebuilt.oem",
718        ]
719        .iter()
720        .collect();
721
722        let out_path = traj.to_oem_file(path.clone(), cfg).unwrap();
723        // And reload, make sure we have the same data.
724        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
725
726        assert_eq!(traj_reloaded, traj);
727
728        // Now export after trimming one state on either end
729        let cfg = ExportCfg::builder()
730            .timestamp(true)
731            .metadata(HashMap::from([
732                ("originator".to_string(), "Test suite".to_string()),
733                ("object_name".to_string(), "TEST_OBJ".to_string()),
734            ]))
735            .step(20.seconds())
736            .start_epoch(traj.first().orbit.epoch + 1.seconds())
737            .end_epoch(traj.last().orbit.epoch - 1.seconds())
738            .build();
739        let out_path = traj.to_oem_file(path, cfg).unwrap();
740        // And reload, make sure we have the same data.
741        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
742
743        // Note that the number of states has changed because we interpolated with a step similar to the original one but
744        // we started with a different time.
745        assert_eq!(traj_reloaded.states.len(), traj.states.len() - 1);
746        assert_eq!(
747            traj_reloaded.first().orbit.epoch,
748            traj.first().orbit.epoch + 1.seconds()
749        );
750        // Note: because we used a fixed step, the last epoch is actually an offset of step size - end offset
751        // from the original trajectory
752        assert_eq!(
753            traj_reloaded.last().orbit.epoch,
754            traj.last().orbit.epoch - 19.seconds()
755        );
756    }
757
758    #[test]
759    fn test_moon_frame_long_prop() {
760        use std::path::PathBuf;
761
762        let manifest_dir =
763            PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
764
765        let almanac = Almanac::new(
766            &manifest_dir
767                .clone()
768                .join("data/01_planetary/pck08.pca")
769                .to_string_lossy(),
770        )
771        .unwrap()
772        .load(
773            &manifest_dir
774                .join("data/01_planetary/de440s.bsp")
775                .to_string_lossy(),
776        )
777        .unwrap();
778
779        let epoch = Epoch::from_str("2022-06-13T12:00:00").unwrap();
780        let orbit = Orbit::try_keplerian_altitude(
781            350.0,
782            0.02,
783            30.0,
784            45.0,
785            85.0,
786            0.0,
787            epoch,
788            almanac.frame_from_uid(MOON_J2000).unwrap(),
789        )
790        .unwrap();
791
792        let mut traj =
793            Propagator::default_dp78(SpacecraftDynamics::new(OrbitalDynamics::two_body()))
794                .with(orbit.into(), Arc::new(almanac))
795                .for_duration_with_traj(45.days())
796                .unwrap()
797                .1;
798        // Set the name of this object
799        traj.name = Some("TEST_MOON_OBJ".to_string());
800
801        // Export CCSDS OEM file
802        let path: PathBuf = [
803            env!("CARGO_MANIFEST_DIR"),
804            "data",
805            "04_output",
806            "moon_45days.oem",
807        ]
808        .iter()
809        .collect();
810
811        let out_path = traj.to_oem_file(path, ExportCfg::default()).unwrap();
812
813        // And reload
814        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
815
816        assert_eq!(traj, traj_reloaded);
817    }
818}
nyx_space/md/trajectory/sc_traj.rs

nyx_space/md/trajectory/
sc_traj.rs
