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