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