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