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}