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 file = File::open(path).context(StdIOSnafu {
472            action: "opening output trajectory file",
473        })?;
474
475        let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
476
477        let reader = builder.build().context(ParquetSnafu {
478            action: "building output trajectory file",
479        })?;
480
481        for field in &reader.schema().fields {
482            if field.name().as_str() == "Epoch (UTC)" {
483                has_epoch = true;
484            } else {
485                for potential_field in &mut found_fields {
486                    if field.name() == potential_field.0.to_field(None).name() {
487                        potential_field.1 = true;
488                        if potential_field.0 != StateParameter::PropMass {
489                            if let Some(frame_info) = field.metadata().get("Frame") {
490                                // Frame is expected to be serialized as Dhall.
491                                match serde_dhall::from_str(frame_info).parse::<Frame>() {
492                                    Err(e) => {
493                                        return Err(InputOutputError::ParseDhall {
494                                            data: frame_info.to_string(),
495                                            err: format!("{e}"),
496                                        })
497                                    }
498                                    Ok(deser_frame) => frame = Some(deser_frame),
499                                };
500                            }
501                        }
502                        break;
503                    }
504                }
505            }
506        }
507
508        ensure!(
509            has_epoch,
510            MissingDataSnafu {
511                which: "Epoch (UTC)"
512            }
513        );
514
515        ensure!(
516            frame.is_some(),
517            MissingDataSnafu {
518                which: "Frame in metadata"
519            }
520        );
521
522        for (field, exists) in found_fields.iter().take(found_fields.len() - 1) {
523            ensure!(
524                exists,
525                MissingDataSnafu {
526                    which: format!("Missing `{}` field", field.to_field(None).name())
527                }
528            );
529        }
530
531        let sc_compat = found_fields.last().unwrap().1;
532
533        let expected_type = std::any::type_name::<Spacecraft>()
534            .split("::")
535            .last()
536            .unwrap();
537
538        if expected_type == "Spacecraft" {
539            ensure!(
540                sc_compat,
541                MissingDataSnafu {
542                    which: format!(
543                        "Missing `{}` field",
544                        found_fields.last().unwrap().0.to_field(None).name()
545                    )
546                }
547            );
548        } else if sc_compat {
549            // Not a spacecraft, remove the prop mass
550            if let Some(last_field) = found_fields.last_mut() {
551                if last_field.0 == StateParameter::PropMass && last_field.1 {
552                    last_field.1 = false;
553                }
554            }
555        }
556
557        // At this stage, we know that the measurement is valid and the conversion is supported.
558        let mut traj = Traj::default();
559
560        // Now convert each batch on the fly
561        for maybe_batch in reader {
562            let batch = maybe_batch.unwrap();
563
564            let epochs = batch
565                .column_by_name("Epoch (UTC)")
566                .unwrap()
567                .as_any()
568                .downcast_ref::<StringArray>()
569                .unwrap();
570
571            let mut shared_data = vec![];
572
573            for (field, _) in found_fields.iter().take(found_fields.len() - 1) {
574                shared_data.push(
575                    batch
576                        .column_by_name(field.to_field(None).name())
577                        .unwrap()
578                        .as_any()
579                        .downcast_ref::<Float64Array>()
580                        .unwrap(),
581                );
582            }
583
584            if expected_type == "Spacecraft" {
585                // Read the prop only if this is a spacecraft we're building
586                shared_data.push(
587                    batch
588                        .column_by_name("prop_mass (kg)")
589                        .unwrap()
590                        .as_any()
591                        .downcast_ref::<Float64Array>()
592                        .unwrap(),
593                );
594            }
595
596            // Grab the frame -- it should have been serialized with all of the data so we don't need to reload it.
597
598            // Build the states
599            for i in 0..batch.num_rows() {
600                let mut state = Spacecraft::zeros();
601                state.set_epoch(Epoch::from_gregorian_str(epochs.value(i)).map_err(|e| {
602                    InputOutputError::Inconsistency {
603                        msg: format!("{e} when parsing epoch"),
604                    }
605                })?);
606                state.set_frame(frame.unwrap()); // We checked it was set above with an ensure! call
607                state.unset_stm(); // We don't have any STM data, so let's unset this.
608
609                for (j, (param, exists)) in found_fields.iter().enumerate() {
610                    if *exists {
611                        state.set_value(*param, shared_data[j].value(i)).unwrap();
612                    }
613                }
614
615                traj.states.push(state);
616            }
617        }
618
619        // Remove any duplicates that may exist in the imported trajectory.
620        traj.finalize();
621
622        Ok(traj)
623    }
624}
625
626#[cfg(test)]
627mod ut_ccsds_oem {
628
629    use crate::md::prelude::{OrbitalDynamics, Propagator, SpacecraftDynamics};
630    use crate::time::{Epoch, TimeUnits};
631    use crate::Spacecraft;
632    use crate::{io::ExportCfg, md::prelude::Traj, Orbit};
633    use anise::almanac::Almanac;
634    use anise::constants::frames::MOON_J2000;
635    use pretty_env_logger;
636    use std::env;
637    use std::str::FromStr;
638    use std::sync::Arc;
639    use std::{collections::HashMap, path::PathBuf};
640
641    #[test]
642    fn test_load_oem_leo() {
643        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
644        let path: PathBuf = [
645            env!("CARGO_MANIFEST_DIR"),
646            "data",
647            "tests",
648            "ccsds",
649            "oem",
650            "LEO_10s.oem",
651        ]
652        .iter()
653        .collect();
654
655        let _ = pretty_env_logger::try_init();
656
657        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
658
659        // This trajectory has two duplicate epochs, which should be removed by the call to finalize()
660        assert_eq!(traj.states.len(), 361);
661        assert_eq!(traj.name.unwrap(), "TEST_OBJ".to_string());
662    }
663
664    #[test]
665    fn test_load_oem_meo() {
666        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
667        let path: PathBuf = [
668            env!("CARGO_MANIFEST_DIR"),
669            "data",
670            "tests",
671            "ccsds",
672            "oem",
673            "MEO_60s.oem",
674        ]
675        .iter()
676        .collect();
677
678        let _ = pretty_env_logger::try_init();
679
680        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
681
682        assert_eq!(traj.states.len(), 61);
683        assert_eq!(traj.name.unwrap(), "TEST_OBJ".to_string());
684    }
685
686    #[test]
687    fn test_load_oem_geo() {
688        use pretty_env_logger;
689        use std::env;
690
691        // All three samples were taken from https://github.com/bradsease/oem/blob/main/tests/samples/real/
692        let path: PathBuf = [
693            env!("CARGO_MANIFEST_DIR"),
694            "data",
695            "tests",
696            "ccsds",
697            "oem",
698            "GEO_20s.oem",
699        ]
700        .iter()
701        .collect();
702
703        let _ = pretty_env_logger::try_init();
704
705        let traj: Traj<Spacecraft> = Traj::from_oem_file(path, None).unwrap();
706
707        assert_eq!(traj.states.len(), 181);
708        assert_eq!(traj.name.as_ref().unwrap(), &"TEST_OBJ".to_string());
709
710        // Reexport this to CCSDS.
711        let cfg = ExportCfg::builder()
712            .timestamp(true)
713            .metadata(HashMap::from([
714                ("originator".to_string(), "Test suite".to_string()),
715                ("object_name".to_string(), "TEST_OBJ".to_string()),
716            ]))
717            .build();
718
719        let path: PathBuf = [
720            env!("CARGO_MANIFEST_DIR"),
721            "output_data",
722            "GEO_20s_rebuilt.oem",
723        ]
724        .iter()
725        .collect();
726
727        let out_path = traj.to_oem_file(path.clone(), cfg).unwrap();
728        // And reload, make sure we have the same data.
729        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
730
731        assert_eq!(traj_reloaded, traj);
732
733        // Now export after trimming one state on either end
734        let cfg = ExportCfg::builder()
735            .timestamp(true)
736            .metadata(HashMap::from([
737                ("originator".to_string(), "Test suite".to_string()),
738                ("object_name".to_string(), "TEST_OBJ".to_string()),
739            ]))
740            .step(20.seconds())
741            .start_epoch(traj.first().orbit.epoch + 1.seconds())
742            .end_epoch(traj.last().orbit.epoch - 1.seconds())
743            .build();
744        let out_path = traj.to_oem_file(path, cfg).unwrap();
745        // And reload, make sure we have the same data.
746        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
747
748        // Note that the number of states has changed because we interpolated with a step similar to the original one but
749        // we started with a different time.
750        assert_eq!(traj_reloaded.states.len(), traj.states.len() - 1);
751        assert_eq!(
752            traj_reloaded.first().orbit.epoch,
753            traj.first().orbit.epoch + 1.seconds()
754        );
755        // Note: because we used a fixed step, the last epoch is actually an offset of step size - end offset
756        // from the original trajectory
757        assert_eq!(
758            traj_reloaded.last().orbit.epoch,
759            traj.last().orbit.epoch - 19.seconds()
760        );
761    }
762
763    #[test]
764    fn test_moon_frame_long_prop() {
765        use std::path::PathBuf;
766
767        let manifest_dir =
768            PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
769
770        let almanac = Almanac::new(
771            &manifest_dir
772                .clone()
773                .join("data/pck08.pca")
774                .to_string_lossy(),
775        )
776        .unwrap()
777        .load(&manifest_dir.join("data/de440s.bsp").to_string_lossy())
778        .unwrap();
779
780        let epoch = Epoch::from_str("2022-06-13T12:00:00").unwrap();
781        let orbit = Orbit::try_keplerian_altitude(
782            350.0,
783            0.02,
784            30.0,
785            45.0,
786            85.0,
787            0.0,
788            epoch,
789            almanac.frame_from_uid(MOON_J2000).unwrap(),
790        )
791        .unwrap();
792
793        let mut traj =
794            Propagator::default_dp78(SpacecraftDynamics::new(OrbitalDynamics::two_body()))
795                .with(orbit.into(), Arc::new(almanac))
796                .for_duration_with_traj(45.days())
797                .unwrap()
798                .1;
799        // Set the name of this object
800        traj.name = Some("TEST_MOON_OBJ".to_string());
801
802        // Export CCSDS OEM file
803        let path: PathBuf = [env!("CARGO_MANIFEST_DIR"), "output_data", "moon_45days.oem"]
804            .iter()
805            .collect();
806
807        let out_path = traj.to_oem_file(path, ExportCfg::default()).unwrap();
808
809        // And reload
810        let traj_reloaded: Traj<Spacecraft> = Traj::from_oem_file(out_path, None).unwrap();
811
812        assert_eq!(traj, traj_reloaded);
813    }
814}
nyx_space/md/trajectory/sc_traj.rs

nyx_space/md/trajectory/
sc_traj.rs
