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 "{filepath} only contained (degree, order) of ({max_degree}, {max_order}) instead of requested ({degree}, {order})"
307 );
308 } else {
309 info!("{filepath} loaded with (degree, order) = ({degree}, {order})");
310 }
311 Ok(HarmonicsMem {
312 degree: max_degree,
313 order: max_order,
314 c_nm: c_nm_mat,
315 s_nm: s_nm_mat,
316 })
317 }
318
319 /// `load` handles the actual loading in memory.
320 fn load(
321 gunzipped: bool,
322 skip_first_line: bool,
323 degree: usize,
324 order: usize,
325 filepath: &str,
326 ) -> Result<HarmonicsMem, NyxError> {
327 let mut f = File::open(filepath).map_err(|_| NyxError::FileUnreadable {
328 msg: format!("File not found: {filepath}"),
329 })?;
330 let mut buffer = vec![0; 0];
331 if gunzipped {
332 let mut d = GzDecoder::new(f);
333 d.read_to_end(&mut buffer)
334 .map_err(|_| NyxError::FileUnreadable {
335 msg: "could not read file as gunzip".to_string(),
336 })?;
337 } else {
338 f.read_to_end(&mut buffer)
339 .map_err(|_| NyxError::FileUnreadable {
340 msg: "could not read file to end".to_string(),
341 })?;
342 }
343
344 let data_as_str = String::from_utf8(buffer).map_err(|_| NyxError::FileUnreadable {
345 msg: "could not decode file contents as utf8".to_string(),
346 })?;
347
348 let mut c_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
349 let mut s_nm_mat = DMatrix::from_element(degree + 1, degree + 1, 0.0);
350
351 let mut max_degree: usize = 0;
352 let mut max_order: usize = 0;
353 for (lno, line) in data_as_str.split('\n').enumerate() {
354 if lno == 0 && skip_first_line {
355 continue;
356 }
357 // These variables need to be declared as mutable because rustc does not know
358 // we won't match each ino more than once.
359 let mut cur_order: usize = 0;
360 let mut cur_degree: usize = 0;
361 let mut c_nm: f64 = 0.0;
362 let mut s_nm: f64 = 0.0;
363 for (ino, item) in line.replace(',', " ").split_whitespace().enumerate() {
364 match ino {
365 0 => match usize::from_str(item) {
366 Ok(val) => cur_degree = val,
367 Err(_) => {
368 return Err(NyxError::FileUnreadable {
369 msg: format!(
370 "Harmonics file:
371 could not parse degree on line {lno} (`{item}`)",
372 ),
373 });
374 }
375 },
376 1 => match usize::from_str(item) {
377 Ok(val) => cur_order = val,
378 Err(_) => {
379 return Err(NyxError::FileUnreadable {
380 msg: format!(
381 "Harmonics file:
382 could not parse order on line {lno} (`{item}`)"
383 ),
384 });
385 }
386 },
387 2 => match f64::from_str(&item.replace('D', "E")) {
388 Ok(val) => c_nm = val,
389 Err(_) => {
390 return Err(NyxError::FileUnreadable {
391 msg: format!(
392 "Harmonics file:
393 could not parse C_nm `{item}` on line {lno}"
394 ),
395 });
396 }
397 },
398 3 => match f64::from_str(&item.replace('D', "E")) {
399 Ok(val) => s_nm = val,
400 Err(_) => {
401 return Err(NyxError::FileUnreadable {
402 msg: format!(
403 "Harmonics file:
404 could not parse S_nm `{item}` on line {lno}"
405 ),
406 });
407 }
408 },
409 _ => break, // We aren't storing the covariance of these harmonics
410 }
411 }
412
413 if cur_degree > degree {
414 // The file is organized by degree, so once we've passed the maximum degree we want,
415 // we can safely stop reading the file.
416 break;
417 }
418
419 // Only insert this data into the hashmap if it's within the required order as well
420 if cur_order <= order {
421 c_nm_mat[(cur_degree, cur_order)] = c_nm;
422 s_nm_mat[(cur_degree, cur_order)] = s_nm;
423 }
424 // This serves as a warning.
425 max_order = if cur_order > max_order {
426 cur_order
427 } else {
428 max_order
429 };
430 max_degree = if cur_degree > max_degree {
431 cur_degree
432 } else {
433 max_degree
434 };
435 }
436 if max_degree < degree || max_order < order {
437 warn!(
438 "{filepath} only contained (degree, order) of ({max_degree}, {max_order}) instead of requested ({degree}, {order})",
439 );
440 } else {
441 info!("{filepath} loaded with (degree, order) = ({degree}, {order})");
442 }
443 Ok(HarmonicsMem {
444 order: max_order,
445 degree: max_degree,
446 c_nm: c_nm_mat,
447 s_nm: s_nm_mat,
448 })
449 }
450
451 /// Returns the maximum order of this gravity potential storage (Jnm=Jn2,Jn3...)
452 pub fn max_order_m(&self) -> usize {
453 self.order
454 }
455
456 /// Returns the maximum degree of this gravity potential storage (Jn=J2,J3...)
457 pub fn max_degree_n(&self) -> usize {
458 self.degree
459 }
460
461 /// Returns the C_nm and S_nm for the provided order and degree.
462 pub fn cs_nm(&self, degree: usize, order: usize) -> (f64, f64) {
463 (self.c_nm[(degree, order)], self.s_nm[(degree, order)])
464 }
465}
466
467#[test]
468fn test_load_harmonic_files() {
469 HarmonicsMem::from_cof("data/01_planetary/JGM3.cof.gz", 50, 50, true)
470 .expect("could not load JGM3");
471
472 HarmonicsMem::from_egm(
473 "data/01_planetary/EGM2008_to2190_TideFree.gz",
474 120,
475 120,
476 true,
477 )
478 .expect("could not load EGM2008");
479
480 HarmonicsMem::from_shadr(
481 "data/01_planetary/Luna_jggrx_1500e_sha.tab.gz",
482 1500,
483 1500,
484 true,
485 )
486 .expect("could not load jggrx");
487}