nyx_space/io/gravity.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 crate::linalg::DMatrix;
20use crate::NyxError;
21use flate2::read::GzDecoder;
22use log::{info, warn};
23use serde::{Deserialize, Serialize};
24use serde_dhall::{SimpleType, StaticType};
25use std::collections::HashMap;
26use std::fmt::Debug;
27use std::fs::File;
28use std::io::prelude::*;
29use std::path::{Path, PathBuf};
30use std::str::FromStr;
31
32/// Configuration holder for gravity field.
33///
34/// Data is first loaded as a SHADR, if that fails, Nyx will try to load it as a COF file.
35#[derive(Clone, Serialize, Deserialize, Debug)]
36pub struct GravityFieldConfig {
37 /// Path to the file, relative to the current working director
38 pub filepath: PathBuf,
39 /// Set to true if the data is gunzipped
40 pub gunzipped: bool,
41 /// Desired degree
42 pub degree: usize,
43 /// Desired order
44 pub order: usize,
45}
46
47/// `HarmonicsMem` loads the requested gravity potential files and stores them in memory (in a HashMap).
48///
49/// WARNING: This memory backend may require a lot of RAM (e.g. EMG2008 2190x2190 requires nearly 400 MB of RAM).
50#[derive(Clone)]
51pub struct GravityFieldData {
52 degree: usize,
53 order: usize,
54 c_nm: DMatrix<f64>,
55 s_nm: DMatrix<f64>,
56}
57
58impl GravityFieldData {
59 pub fn from_config(cfg: GravityFieldConfig) -> Result<Self, NyxError> {
60 if !cfg.gunzipped && cfg.filepath.ends_with(".cof")
61 || cfg.gunzipped && cfg.filepath.ends_with(".cof.gz")
62 {
63 Self::from_cof(cfg.filepath, cfg.degree, cfg.order, cfg.gunzipped)
64 } else {
65 Self::from_shadr(cfg.filepath, cfg.degree, cfg.order, cfg.gunzipped)
66 }
67 }
68
69 /// Initialize `HarmonicsMem` with a custom J2 value
70 pub fn from_j2(j2: f64) -> GravityFieldData {
71 let mut c_nm = DMatrix::from_element(3, 3, 0.0);
72 c_nm[(2, 0)] = j2;
73
74 GravityFieldData {
75 degree: 2,
76 order: 0,
77 c_nm,
78 s_nm: DMatrix::from_element(3, 3, 0.0),
79 }
80 }
81
82 /// Initialize `HarmonicsMem` as an EARTH J<sub>2</sub> only using the JGM3 model (available in GMAT)
83 ///
84 /// Use the embedded Earth parameter. If others are needed, load from `from_shadr` or `from_egm`.
85 /// *WARNING:* This is an EARTH gravity model, and _should not_ be used around any other body.
86 pub fn j2_jgm3() -> GravityFieldData {
87 Self::from_j2(-4.841_653_748_864_70e-04)
88 }
89
90 /// Initialize `HarmonicsMem` as an EARTH J<sub>2</sub> only using the JGM2 model (available in GMAT)
91 ///
92 /// Use the embedded Earth parameter. If others are needed, load from `from_shadr` or `from_egm`.
93 /// *WARNING:* This is an EARTH gravity model, and _should not_ be used around any other body.
94 pub fn j2_jgm2() -> GravityFieldData {
95 Self::from_j2(-4.841_653_9e-04)
96 }
97
98 /// Initialize `HarmonicsMem` as J<sub>2</sub> only using the EGM2008 model (from the GRACE mission, best model as of 2018)
99 ///
100 /// *WARNING:* This is an EARTH gravity model, and _should not_ be used around any other body.
101 pub fn j2_egm2008() -> GravityFieldData {
102 Self::from_j2(-0.484_165_143_790_815e-03)
103 }
104
105 /// Initialize `HarmonicsMem` from the file path (must be a gunzipped file)
106 ///
107 /// Gravity models provided by `nyx`:
108 /// + EMG2008 to 2190 for Earth (tide free)
109 /// + Moon to 1500 (from SHADR file)
110 /// + Mars to 120 (from SHADR file)
111 /// + Venus to 150 (from SHADR file)
112 pub fn from_shadr<P: AsRef<Path> + Debug>(
113 filepath: P,
114 degree: usize,
115 order: usize,
116 gunzipped: bool,
117 ) -> Result<GravityFieldData, NyxError> {
118 Self::load(
119 filepath, gunzipped, true, //SHADR has a header which we ignore
120 degree, order,
121 )
122 }
123
124 pub fn from_cof<P: AsRef<Path> + Debug>(
125 filepath: P,
126 degree: usize,
127 order: usize,
128 gunzipped: bool,
129 ) -> Result<GravityFieldData, NyxError> {
130 let mut f = File::open(&filepath).map_err(|_| NyxError::FileUnreadable {
131 msg: format!("File not found: {filepath:?}"),
132 })?;
133 let mut buffer = vec![0; 0];
134 if gunzipped {
135 let mut d = GzDecoder::new(f);
136 d.read_to_end(&mut buffer)
137 .map_err(|_| NyxError::FileUnreadable {
138 msg: "could not read file as gunzip".to_string(),
139 })?;
140 } else {
141 f.read_to_end(&mut buffer)
142 .map_err(|_| NyxError::FileUnreadable {
143 msg: "could not read file to end".to_string(),
144 })?;
145 }
146
147 let data_as_str = String::from_utf8(buffer).map_err(|_| NyxError::FileUnreadable {
148 msg: "could not decode file contents as utf8".to_string(),
149 })?;
150
151 // Since the COF files are so specific, we just code everything up in here.
152
153 let mut c_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
154 let mut s_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
155 let mut max_order: usize = 0;
156 let mut max_degree: usize = 0;
157 for (lno, line) in data_as_str.split('\n').enumerate() {
158 if line.is_empty() || !line.starts_with('R') {
159 continue; // This is either a comment, a header or "END"
160 }
161 // These variables need to be declared as mutable because rustc does not know
162 // we nwon't match each ino more than once.
163 let mut cur_degree: usize = 0;
164 let mut cur_order: usize = 0;
165 let mut c_nm: f64 = 0.0;
166 let mut s_nm: f64 = 0.0;
167 for (ino, item) in line.split_whitespace().enumerate() {
168 match ino {
169 0 => continue, // We need this so we don't break at every first item
170 1 => match usize::from_str(item) {
171 Ok(val) => cur_degree = val,
172 Err(_) => {
173 return Err(NyxError::FileUnreadable {
174 msg: format!(
175 "Harmonics file:
176 could not parse degree `{item}` on line {lno}"
177 ),
178 });
179 }
180 },
181 2 => match usize::from_str(item) {
182 Ok(val) => cur_order = val,
183 Err(_) => {
184 return Err(NyxError::FileUnreadable {
185 msg: format!(
186 "Harmonics file:
187 could not parse order `{item}` on line {lno}"
188 ),
189 });
190 }
191 },
192 3 => {
193 // If we are at degree zero, then there is only one item, so we can parse that and
194 // set the S_nm to zero.
195 if degree == 0 {
196 s_nm = 0.0;
197 match f64::from_str(item) {
198 Ok(val) => c_nm = val,
199 Err(_) => {
200 return Err(NyxError::FileUnreadable {
201 msg: format!(
202 "Harmonics file:
203 could not parse C_nm `{item}` on line {lno}"
204 ),
205 });
206 }
207 }
208 } else {
209 // There is a space as a delimiting character between the C_nm and S_nm only if the S_nm
210 // is a positive number, otherwise, they are continuous (what a great format).
211 if (item.matches('-').count() == 3 && !item.starts_with('-'))
212 || item.matches('-').count() == 4
213 {
214 // Now we have two items concatenated into one... great
215 let parts: Vec<&str> = item.split('-').collect();
216 if parts.len() == 5 {
217 // That mean we have five minus signs, so both the C and S are negative.
218 let c_nm_str = "-".to_owned() + parts[1] + "-" + parts[2];
219 match f64::from_str(&c_nm_str) {
220 Ok(val) => c_nm = val,
221 Err(_) => {
222 return Err(NyxError::FileUnreadable {
223 msg: format!(
224 "Harmonics file:
225 could not parse C_nm `{item}` on line {lno}"
226 ),
227 });
228 }
229 }
230 // That mean we have five minus signs, so both the C and S are negative.
231 let s_nm_str = "-".to_owned() + parts[3] + "-" + parts[4];
232 match f64::from_str(&s_nm_str) {
233 Ok(val) => s_nm = val,
234 Err(_) => {
235 return Err(NyxError::FileUnreadable {
236 msg: format!(
237 "Harmonics file:
238 could not parse S_nm `{item}` on line {lno}"
239 ),
240 });
241 }
242 }
243 } else {
244 // That mean we have fouve minus signs, and since both values are concatenated, C_nm is positive and S_nm is negative
245 let c_nm_str = parts[0].to_owned() + "-" + parts[1];
246 match f64::from_str(&c_nm_str) {
247 Ok(val) => c_nm = val,
248 Err(_) => {
249 return Err(NyxError::FileUnreadable {
250 msg: format!(
251 "Harmonics file:
252 could not parse C_nm `{item}` on line {lno}"
253 ),
254 });
255 }
256 }
257 // That mean we have five minus signs, so both the C and S are negative.
258 let s_nm_str = "-".to_owned() + parts[2] + "-" + parts[3];
259 match f64::from_str(&s_nm_str) {
260 Ok(val) => s_nm = val,
261 Err(_) => {
262 return Err(NyxError::FileUnreadable {
263 msg: format!(
264 "Harmonics file:
265 could not parse S_nm `{item}` on line {lno}"
266 ),
267 });
268 }
269 }
270 }
271 } else {
272 // We only have the first item, and that's the C_nm
273 match f64::from_str(item) {
274 Ok(val) => c_nm = val,
275 Err(_) => {
276 return Err(NyxError::FileUnreadable {
277 msg: format!(
278 "Harmonics file:
279 could not parse C_nm `{item}` on line {lno}"
280 ),
281 });
282 }
283 }
284 }
285 }
286 }
287 4 => match f64::from_str(item) {
288 // If this exists, then the S_nm is positive.
289 Ok(val) => s_nm = val,
290 Err(_) => {
291 return Err(NyxError::FileUnreadable {
292 msg: format!(
293 "Harmonics file:
294 could not parse S_nm `{item}` on line {lno}"
295 ),
296 });
297 }
298 },
299 _ => break, // We aren't storing the covariance of these harmonics
300 }
301 }
302
303 if cur_degree > degree {
304 // The file is organized by degree, so once we've passed the maximum degree we want,
305 // we can safely stop reading the file.
306 break;
307 }
308
309 // Only insert this data into the hashmap if it's within the required order as well
310 if cur_order <= order {
311 c_nm_mat[(cur_degree, cur_order)] = c_nm;
312 s_nm_mat[(cur_degree, cur_order)] = s_nm;
313 }
314 // This serves as a warning.
315 max_order = if cur_order > max_order {
316 cur_order
317 } else {
318 max_order
319 };
320 max_degree = if cur_degree > max_degree {
321 cur_degree
322 } else {
323 max_degree
324 };
325 }
326 if max_degree < degree || max_order < order {
327 warn!(
328 "{filepath:?} only contained (degree, order) of ({max_degree}, {max_order}) instead of requested ({degree}, {order})"
329 );
330 } else {
331 info!("{filepath:?} loaded with (degree, order) = ({degree}, {order})");
332 }
333 Ok(GravityFieldData {
334 degree: max_degree,
335 order: max_order,
336 c_nm: c_nm_mat,
337 s_nm: s_nm_mat,
338 })
339 }
340
341 /// `load` handles the actual loading in memory.
342 fn load<P: AsRef<Path> + Debug>(
343 filepath: P,
344 gunzipped: bool,
345 skip_first_line: bool,
346 degree: usize,
347 order: usize,
348 ) -> Result<GravityFieldData, NyxError> {
349 let mut f = File::open(&filepath).map_err(|_| NyxError::FileUnreadable {
350 msg: format!("File not found: {filepath:?}"),
351 })?;
352 let mut buffer = vec![0; 0];
353 if gunzipped {
354 let mut d = GzDecoder::new(f);
355 d.read_to_end(&mut buffer)
356 .map_err(|_| NyxError::FileUnreadable {
357 msg: "could not read file as gunzip".to_string(),
358 })?;
359 } else {
360 f.read_to_end(&mut buffer)
361 .map_err(|_| NyxError::FileUnreadable {
362 msg: "could not read file to end".to_string(),
363 })?;
364 }
365
366 let data_as_str = String::from_utf8(buffer).map_err(|_| NyxError::FileUnreadable {
367 msg: "could not decode file contents as utf8".to_string(),
368 })?;
369
370 let mut c_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
371 let mut s_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
372
373 let mut max_degree: usize = 0;
374 let mut max_order: usize = 0;
375 for (lno, line) in data_as_str.split('\n').enumerate() {
376 if lno == 0 && skip_first_line {
377 continue;
378 }
379 // These variables need to be declared as mutable because rustc does not know
380 // we won't match each ino more than once.
381 let mut cur_order: usize = 0;
382 let mut cur_degree: usize = 0;
383 let mut c_nm: f64 = 0.0;
384 let mut s_nm: f64 = 0.0;
385 for (ino, item) in line.replace(',', " ").split_whitespace().enumerate() {
386 match ino {
387 0 => match usize::from_str(item) {
388 Ok(val) => cur_degree = val,
389 Err(_) => {
390 return Err(NyxError::FileUnreadable {
391 msg: format!(
392 "Harmonics file:
393 could not parse degree on line {lno} (`{item}`)",
394 ),
395 });
396 }
397 },
398 1 => match usize::from_str(item) {
399 Ok(val) => cur_order = val,
400 Err(_) => {
401 return Err(NyxError::FileUnreadable {
402 msg: format!(
403 "Harmonics file:
404 could not parse order on line {lno} (`{item}`)"
405 ),
406 });
407 }
408 },
409 2 => match f64::from_str(&item.replace('D', "E")) {
410 Ok(val) => c_nm = val,
411 Err(_) => {
412 return Err(NyxError::FileUnreadable {
413 msg: format!(
414 "Harmonics file:
415 could not parse C_nm `{item}` on line {lno}"
416 ),
417 });
418 }
419 },
420 3 => match f64::from_str(&item.replace('D', "E")) {
421 Ok(val) => s_nm = val,
422 Err(_) => {
423 return Err(NyxError::FileUnreadable {
424 msg: format!(
425 "Harmonics file:
426 could not parse S_nm `{item}` on line {lno}"
427 ),
428 });
429 }
430 },
431 _ => break, // We aren't storing the covariance of these harmonics
432 }
433 }
434
435 if cur_degree > degree {
436 // The file is organized by degree, so once we've passed the maximum degree we want,
437 // we can safely stop reading the file.
438 break;
439 }
440
441 // Only insert this data into the hashmap if it's within the required order as well
442 if cur_order <= order {
443 c_nm_mat[(cur_degree, cur_order)] = c_nm;
444 s_nm_mat[(cur_degree, cur_order)] = s_nm;
445 }
446 // This serves as a warning.
447 max_order = if cur_order > max_order {
448 cur_order
449 } else {
450 max_order
451 };
452 max_degree = if cur_degree > max_degree {
453 cur_degree
454 } else {
455 max_degree
456 };
457 }
458 if max_degree < degree || max_order < order {
459 warn!(
460 "{filepath:?} only contained (degree, order) of ({max_degree}, {max_order}) instead of requested ({degree}, {order})",
461 );
462 } else {
463 info!("{filepath:?} loaded with (degree, order) = ({degree}, {order})");
464 }
465 Ok(GravityFieldData {
466 order: max_order,
467 degree: max_degree,
468 c_nm: c_nm_mat,
469 s_nm: s_nm_mat,
470 })
471 }
472
473 /// Returns the maximum order of this gravity potential storage (Jnm=Jn2,Jn3...)
474 pub fn max_order_m(&self) -> usize {
475 self.order
476 }
477
478 /// Returns the maximum degree of this gravity potential storage (Jn=J2,J3...)
479 pub fn max_degree_n(&self) -> usize {
480 self.degree
481 }
482
483 /// Returns the C_nm and S_nm for the provided order and degree.
484 pub fn cs_nm(&self, degree: usize, order: usize) -> (f64, f64) {
485 (self.c_nm[(degree, order)], self.s_nm[(degree, order)])
486 }
487}
488
489impl StaticType for GravityFieldConfig {
490 fn static_type() -> SimpleType {
491 let mut fields = HashMap::new();
492
493 fields.insert("filepath".to_string(), String::static_type());
494 fields.insert("gunzipped".to_string(), bool::static_type());
495 fields.insert("degree".to_string(), usize::static_type());
496 fields.insert("order".to_string(), usize::static_type());
497
498 SimpleType::Record(fields)
499 }
500}
501
502#[test]
503fn test_load_harmonic_files() {
504 GravityFieldData::from_cof("data/01_planetary/JGM3.cof.gz", 50, 50, true)
505 .expect("could not load JGM3");
506
507 GravityFieldData::from_shadr(
508 "data/01_planetary/EGM2008_to2190_TideFree.gz",
509 120,
510 120,
511 true,
512 )
513 .expect("could not load EGM2008");
514
515 GravityFieldData::from_shadr(
516 "data/01_planetary/Luna_jggrx_1500e_sha.tab.gz",
517 1500,
518 1500,
519 true,
520 )
521 .expect("could not load jggrx");
522}