ConstellationVFR/geomag-wmm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

working model

Browse source
This commit is contained in:
Louis 2024-09-28 16:55:51 +02:00
commit 7671a9abb5
6 changed files with 613 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored Normal file
View file

@ -0,0 +1 @@
/target

7
Cargo.lock generated Normal file
View file

@ -0,0 +1,7 @@
# This file is automatically @generated by Cargo.
# It is not intended for manual editing.
version = 3
[[package]]
name = "geomag-wmm"
version = "0.1.0"

6
Cargo.toml Normal file
View file

@ -0,0 +1,6 @@
[package]
name = "geomag-wmm"
version = "0.1.0"
edition = "2021"
[dependencies]

93
WMM.2020.COF Normal file
View file

@ -0,0 +1,93 @@
2020.0 WMM-2020 12/10/2019
1 0 -29404.5 0.0 6.7 0.0
1 1 -1450.7 4652.9 7.7 -25.1
2 0 -2500.0 0.0 -11.5 0.0
2 1 2982.0 -2991.6 -7.1 -30.2
2 2 1676.8 -734.8 -2.2 -23.9
3 0 1363.9 0.0 2.8 0.0
3 1 -2381.0 -82.2 -6.2 5.7
3 2 1236.2 241.8 3.4 -1.0
3 3 525.7 -542.9 -12.2 1.1
4 0 903.1 0.0 -1.1 0.0
4 1 809.4 282.0 -1.6 0.2
4 2 86.2 -158.4 -6.0 6.9
4 3 -309.4 199.8 5.4 3.7
4 4 47.9 -350.1 -5.5 -5.6
5 0 -234.4 0.0 -0.3 0.0
5 1 363.1 47.7 0.6 0.1
5 2 187.8 208.4 -0.7 2.5
5 3 -140.7 -121.3 0.1 -0.9
5 4 -151.2 32.2 1.2 3.0
5 5 13.7 99.1 1.0 0.5
6 0 65.9 0.0 -0.6 0.0
6 1 65.6 -19.1 -0.4 0.1
6 2 73.0 25.0 0.5 -1.8
6 3 -121.5 52.7 1.4 -1.4
6 4 -36.2 -64.4 -1.4 0.9
6 5 13.5 9.0 -0.0 0.1
6 6 -64.7 68.1 0.8 1.0
7 0 80.6 0.0 -0.1 0.0
7 1 -76.8 -51.4 -0.3 0.5
7 2 -8.3 -16.8 -0.1 0.6
7 3 56.5 2.3 0.7 -0.7
7 4 15.8 23.5 0.2 -0.2
7 5 6.4 -2.2 -0.5 -1.2
7 6 -7.2 -27.2 -0.8 0.2
7 7 9.8 -1.9 1.0 0.3
8 0 23.6 0.0 -0.1 0.0
8 1 9.8 8.4 0.1 -0.3
8 2 -17.5 -15.3 -0.1 0.7
8 3 -0.4 12.8 0.5 -0.2
8 4 -21.1 -11.8 -0.1 0.5
8 5 15.3 14.9 0.4 -0.3
8 6 13.7 3.6 0.5 -0.5
8 7 -16.5 -6.9 0.0 0.4
8 8 -0.3 2.8 0.4 0.1
9 0 5.0 0.0 -0.1 0.0
9 1 8.2 -23.3 -0.2 -0.3
9 2 2.9 11.1 -0.0 0.2
9 3 -1.4 9.8 0.4 -0.4
9 4 -1.1 -5.1 -0.3 0.4
9 5 -13.3 -6.2 -0.0 0.1
9 6 1.1 7.8 0.3 -0.0
9 7 8.9 0.4 -0.0 -0.2
9 8 -9.3 -1.5 -0.0 0.5
9 9 -11.9 9.7 -0.4 0.2
10 0 -1.9 0.0 0.0 0.0
10 1 -6.2 3.4 -0.0 -0.0
10 2 -0.1 -0.2 -0.0 0.1
10 3 1.7 3.5 0.2 -0.3
10 4 -0.9 4.8 -0.1 0.1
10 5 0.6 -8.6 -0.2 -0.2
10 6 -0.9 -0.1 -0.0 0.1
10 7 1.9 -4.2 -0.1 -0.0
10 8 1.4 -3.4 -0.2 -0.1
10 9 -2.4 -0.1 -0.1 0.2
10 10 -3.9 -8.8 -0.0 -0.0
11 0 3.0 0.0 -0.0 0.0
11 1 -1.4 -0.0 -0.1 -0.0
11 2 -2.5 2.6 -0.0 0.1
11 3 2.4 -0.5 0.0 0.0
11 4 -0.9 -0.4 -0.0 0.2
11 5 0.3 0.6 -0.1 -0.0
11 6 -0.7 -0.2 0.0 0.0
11 7 -0.1 -1.7 -0.0 0.1
11 8 1.4 -1.6 -0.1 -0.0
11 9 -0.6 -3.0 -0.1 -0.1
11 10 0.2 -2.0 -0.1 0.0
11 11 3.1 -2.6 -0.1 -0.0
12 0 -2.0 0.0 0.0 0.0
12 1 -0.1 -1.2 -0.0 -0.0
12 2 0.5 0.5 -0.0 0.0
12 3 1.3 1.3 0.0 -0.1
12 4 -1.2 -1.8 -0.0 0.1
12 5 0.7 0.1 -0.0 -0.0
12 6 0.3 0.7 0.0 0.0
12 7 0.5 -0.1 -0.0 -0.0
12 8 -0.2 0.6 0.0 0.1
12 9 -0.5 0.2 -0.0 -0.0
12 10 0.1 -0.9 -0.0 -0.0
12 11 -1.1 -0.0 -0.0 0.0
12 12 -0.3 0.5 -0.1 -0.1
999999999999999999999999999999999999999999999999
999999999999999999999999999999999999999999999999

429
src/lib.rs Normal file
View file

@ -0,0 +1,429 @@
#![allow(non_snake_case)]
#![allow(non_upper_case_globals)]
use std::{fs, process};
use std::io::{BufRead,BufReader};
// Defining the constants used in the model
const WMM_VALIDITY_RANGE: f64 = 5.0;
// WGS84 ellipsoid parameters
const A : f64 = 6378137.0; // Equation (8)
const f : f64 = 1.0 / 298.257223563; // Equation (8)
const e_2 : f64 = f * (2.0 - f); // Equation (8)
// Geomagnetic parameters
const a : f64 = 6371200.0; // Equation (4)
/// Calculate the factorial of a number
fn factorial(n: usize) -> f64 {
(1..=n).fold(1.0, |acc, x| acc * x as f64)
}
/// Calculate the Legendre polynomial
///
/// # Arguments
///
/// * `t` - The value to evaluate
/// * `n` - The degree of the polynomial
/// * `m` - The order of the polynomial
///
/// # Returns
///
/// * The value of the Legendre polynomial
///
/// # Reference
///
/// * Heiskanen and Moritz, 1967 (https://archive.org/details/HeiskanenMoritz1967PhysicalGeodesy/page/n33/mode/2up): equation 1-62
fn legendre_polynomial(t: f64, n: usize, m: usize) -> f64 {
let mut P = 0.0;
let mut k = 0;
while k <= (n - m) / 2 {
let num = (-1.0 as f64).powi(k.try_into().unwrap()) * factorial(2 * n - 2 * k) * t.powi(n as i32 - m as i32 - 2 * k as i32);
let den = factorial(k) * factorial(n - k) * factorial(n - m - 2 * k);
P += num / den;
k += 1;
}
P *= (2.0 as f64).powi(-(n as i32)) * (1.0 - t*t).powf(m as f64 / 2.0);
P
}
/// Calculate the Schmidt semi-normalised associated Legendre functions
/// for all degrees and orders up to 14
///
/// # Arguments
///
/// * `mu` - The value to evaluate
///
/// # Reference
///
/// * Equation (5) in "The US/UK Workd Magnetic Model for 2020-2025"
fn schmidt_semi_normalised_associated_legendre(mu: f64) -> [[f64; 14]; 14] {
let mut Psn = [[0.0; 14]; 14];
for n in 0..=13 {
for m in 0..=n {
if m == 0 {
Psn[n][m] = legendre_polynomial(mu, n, m);
} else {
Psn[n][m] = f64::sqrt(2.0 * factorial(n - m) / factorial(n + m)) * legendre_polynomial(mu, n, m);
}
}
}
Psn
}
pub enum Error {
DateOutsideOfValidityRange
}
impl std::fmt::Display for Error {
fn fmt(&self, formatter: &mut std::fmt::Formatter) -> std::fmt::Result {
match self {
Error::DateOutsideOfValidityRange => write!(formatter, "Date outside of validity range!")
}
}
}
pub struct MagneticModel {
from_year: f64,
until_year: f64,
g: [[f64; 13]; 13],
h: [[f64; 13]; 13],
g_sv: [[f64; 13]; 13],
h_sv: [[f64; 13]; 13]
}
impl MagneticModel {
/// Calculate the geomagnetic field components at a given location and time
///
/// # Arguments
///
/// * `h` - The height above the WGS84 ellipsoid [metres]
/// * `phi` - The latitude [degrees]
/// * `lambda` - The longitude [degrees]
/// * `year` - The year [decimal] e.g. 2022.5
///
/// # Returns
///
/// * `X` - The northward intensity [nT]
/// * `Y` - The eastward intensity [nT]
/// * `Z` - The downward intensity [nT]
/// * `H` - The horizontal intensity [nT]
/// * `F` - The total intensity [nT]
/// * `I` - The inclination [degrees]
/// * `D` - The declination [degrees]
/// * `X_dot` - The rate of change of the northward intensity [nT/year]
/// * `Y_dot` - The rate of change of the eastward intensity [nT/year]
/// * `Z_dot` - The rate of change of the downward intensity [nT/year]
/// * `H_dot` - The rate of change of the horizontal intensity [nT/year]
/// * `F_dot` - The rate of change of the total intensity [nT/year]
/// * `I_dot` - The rate of change of the inclination [degrees/year]
/// * `D_dot` - The rate of change of the declination [degrees/year]
///
/// # Note
///
/// The conversion from height above MSL to height above the WGS 84 ellipsoid has a very small effect on the
/// resulting magnetic field values (of the order of 1 nT or less).
pub fn calculate(&self, h: f64, phi: f64, lambda: f64, year: f64) -> Result<(f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64), Error> {
// Translation of coordinates from geodetic to geocentric
let lambda = lambda.to_radians();
let phi = phi.to_radians();
let _R_c = A / f64::sqrt(1.0 - e_2 * f64::sin(phi)*f64::sin(phi)); // Equation (8)
let _p = (_R_c + h) * f64::cos(phi); // Equation (7)
let _z = (_R_c * (1.0 - e_2) + h) * f64::sin(phi); // Equation (7)
let r = f64::sqrt(_p*_p + _z*_z); // Equation (7)
let phi_prime = f64::asin(_z / r); // Equation (7)
let Psn_sinphi = schmidt_semi_normalised_associated_legendre(f64::sin(phi_prime));
let mut g_t = [[0.0; 13]; 13];
let mut h_t = [[0.0; 13]; 13];
let mut dPsn_sinphi_dphi = [[0.0; 13]; 13];
let time_delta = year - self.from_year;
if time_delta < 0.0 || time_delta > WMM_VALIDITY_RANGE {
return Result::Err(Error::DateOutsideOfValidityRange);
}
// Prepare matrices
for n in 0..=12 {
for m in 0..=n {
g_t[n][m] = self.g[n][m] + time_delta * self.g_sv[n][m]; // Equation (9): calculating Gauss coefficients for the desired time
h_t[n][m] = self.h[n][m] + time_delta * self.h_sv[n][m]; // Equation (9): calculating Gauss coefficients for the desired time
dPsn_sinphi_dphi[n][m] = (n as f64 + 1.0) * f64::tan(phi_prime) * Psn_sinphi[n][m] - // Equation (16): used to calculate the field vector components
f64::sqrt((n as f64 + 1.0).powi(2) - (m as f64).powi(2)) / f64::cos(phi_prime) * Psn_sinphi[n+1][m];
}
}
// Field vector components
let mut X_prime = 0.0; // will be computed using Equation (10)
let mut Y_prime = 0.0; // will be computed using Equation (11)
let mut Z_prime = 0.0; // will be computed using Equation (12)
let mut X_dot_prime = 0.0; // will be computed using Equation (13)
let mut Y_dot_prime = 0.0; // will be computed using Equation (14)
let mut Z_dot_prime = 0.0; // will be computed using Equation (15)
for n in 1..=12 {
let mut X_prime_tmp = 0.0;
let mut Y_prime_tmp = 0.0;
let mut Z_prime_tmp = 0.0;
let mut X_dot_prime_tmp = 0.0;
let mut Y_dot_prime_tmp = 0.0;
let mut Z_dot_prime_tmp = 0.0;
for m in 0..=n {
X_prime_tmp += (g_t[n][m] * f64::cos(m as f64 * lambda) + h_t[n][m] * f64::sin(m as f64 * lambda)) * dPsn_sinphi_dphi[n][m];
Y_prime_tmp += m as f64 * (g_t[n][m] * f64::sin(m as f64 * lambda) - h_t[n][m] * f64::cos(m as f64 * lambda)) * Psn_sinphi[n][m];
Z_prime_tmp += (g_t[n][m] * f64::cos(m as f64 * lambda) + h_t[n][m] * f64::sin(m as f64 * lambda)) * Psn_sinphi[n][m];
X_dot_prime_tmp += (self.g_sv[n][m] * f64::cos(m as f64 * lambda) + self.h_sv[n][m] * f64::sin(m as f64 * lambda)) * dPsn_sinphi_dphi[n][m];
Y_dot_prime_tmp += m as f64 * (self.g_sv[n][m] * f64::sin(m as f64 * lambda) - self.h_sv[n][m] * f64::cos(m as f64 * lambda)) * Psn_sinphi[n][m];
Z_dot_prime_tmp += (self.g_sv[n][m] * f64::cos(m as f64 * lambda) + self.h_sv[n][m] * f64::sin(m as f64 * lambda)) * Psn_sinphi[n][m];
}
let k_temp = (a / r).powi(n as i32 + 2);
X_prime += k_temp * X_prime_tmp;
Y_prime += k_temp * Y_prime_tmp;
Z_prime += (n as f64 + 1.0) * k_temp * Z_prime_tmp;
X_dot_prime += k_temp * X_dot_prime_tmp;
Y_dot_prime += k_temp * Y_dot_prime_tmp;
Z_dot_prime += (n as f64 + 1.0) * k_temp * Z_dot_prime_tmp;
}
X_prime *= -1.0;
Y_prime *= 1.0 / f64::cos(phi_prime);
Z_prime *= -1.0;
X_dot_prime *= -1.0;
Y_dot_prime *= 1.0 / f64::cos(phi_prime);
Z_dot_prime *= -1.0;
// Rotate the field vector components to the ellipsoidal reference frame
let X = X_prime * f64::cos(phi_prime - phi) - Z_prime * f64::sin(phi_prime - phi);
let Y = Y_prime;
let Z = X_prime * f64::sin(phi_prime - phi) + Z_prime * f64::cos(phi_prime - phi);
let X_dot = X_dot_prime * f64::cos(phi_prime - phi) - Z_dot_prime * f64::sin(phi_prime - phi);
let Y_dot = Y_dot_prime;
let Z_dot = X_dot_prime * f64::sin(phi_prime - phi) + Z_dot_prime * f64::cos(phi_prime - phi);
// Compute the magnetic elements
let H = f64::sqrt(X*X + Y*Y);
let F = f64::sqrt(H*H + Z*Z);
let I = f64::atan2(Z, H);
let D = f64::atan2(Y, X);
let H_dot = (X*X_dot + Y*Y_dot) / H;
let F_dot = (X*X_dot + Y*Y_dot + Z*Z_dot) / F;
let I_dot = (Z_dot*H - Z*H_dot) / (F*F);
let D_dot = (Y_dot*X - Y*X_dot) / (H*H);
return Result::Ok (
(X, Y, Z, H, F, I.to_degrees(), D.to_degrees(),
X_dot, Y_dot, Z_dot, H_dot, F_dot, I_dot.to_degrees(), D_dot.to_degrees())
);
}
}
/// Function to initialise the magnetic model
///
/// # Arguments
///
/// * `path` - The path to the WMM model file
pub fn initialise_magnetic_model(path: &str) -> MagneticModel {
// Read WMM model file
let model_file = fs::File::open(path).expect("Model file not found!");
let mut model_file = BufReader::new(model_file).lines();
// Parse header line
let from_year = model_file
.next().expect("Model file is empty!")
.expect("Error reading model file!")
.split_whitespace()
.nth(0).expect("Error parsing model file!")
.parse::<f64>().expect("Error parsing model file!");
// Determine the validity range of the model
let until_year = from_year + WMM_VALIDITY_RANGE;
// Read the rest of the file
let mut g = [[0.0; 13]; 13];
let mut h = [[0.0; 13]; 13];
let mut g_sv = [[0.0; 13]; 13];
let mut h_sv = [[0.0; 13]; 13];
for line in model_file {
match line {
Ok(line) => {
if line.starts_with("9999") {
break;
}
let mut line = line.split_whitespace();
let n_line: i32 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
let m_line: i32 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
let g_line: f64 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
let h_line: f64 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
let g_sv_line: f64 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
let h_sv_line: f64 = line.next().expect("Error parsing model file!").parse().expect("Error parsing model file!");
if n_line > 12 {
break;
}
if m_line > n_line || m_line < 0 {
eprintln!("Corrupt record in model file!");
process::exit(1);
}
g[n_line as usize][m_line as usize] = g_line;
g_sv[n_line as usize][m_line as usize] = g_sv_line;
if m_line != 0 {
h[n_line as usize][m_line as usize] = h_line;
h_sv[n_line as usize][m_line as usize] = h_sv_line;
}
},
Err(e) => {
eprintln!("Error reading model file: {:?}", e);
}
}
}
println!("Model is valid from {:04}-01-01 to {:04}-01-01", from_year, until_year);
MagneticModel {
from_year,
until_year,// geomg1
g,
h,
g_sv,
h_sv
}
}
#[cfg(test)]
mod tests {
use super::*;
fn round_to_decimal_places(value: f64, decimal_places: i32) -> f64 {
let multiplier = 10.0_f64.powi(decimal_places);
(value * multiplier).round() / multiplier
}
/// Test function
///
/// # Arguments
///
/// * `year` - The year
/// * `h` - The height above the WGS84 ellipsoid [km]
/// * `phi` - The latitude [degrees]
/// * `l` - The longitude [degrees]
/// * `expected_X` - The expected northward intensity [nT]
/// * `expected_Y` - The expected eastward intensity [nT]
/// * `expected_Z` - The expected downward intensity [nT]
/// * `expected_H` - The expected horizontal intensity [nT]
/// * `expected_F` - The expected total intensity [nT]
/// * `expected_I` - The expected inclination [degrees]
/// * `expected_D` - The expected declination [degrees]
/// * `expected_GV` - The expected grid variation [degrees]
/// * `expected_X_dot` - The expected rate of change of the northward intensity [nT/year]
/// * `expected_Y_dot` - The expected rate of change of the eastward intensity [nT/year]
/// * `expected_Z_dot` - The expected rate of change of the downward intensity [nT/year]
/// * `expected_H_dot` - The expected rate of change of the horizontal intensity [nT/year]
/// * `expected_F_dot` - The expected rate of change of the total intensity [nT/year]
/// * `expected_I_dot` - The expected rate of change of the inclination [degrees/year]
/// * `expected_D_dot` - The expected rate of change of the declination [degrees/year]
///
/// # Note
///
/// The grid variation is not tested as it is not calculated in the model (yet)
fn test(year: f64, h: f64, phi: f64, l: f64,
expected_X: f64, expected_Y: f64, expected_Z: f64, expected_H: f64, expected_F: f64, expected_I: f64, expected_D: f64, expected_GV: f64,
expected_X_dot: f64, expected_Y_dot: f64, expected_Z_dot: f64, expected_H_dot: f64, expected_F_dot: f64, expected_I_dot: f64, expected_D_dot: f64) {
let magnetic_model = initialise_magnetic_model("WMM.2020.COF");
let h = h * 1000.0;
let result = magnetic_model.calculate(h, phi, l, year);
match result {
Ok(results) => {
assert_eq!(round_to_decimal_places(results.0, 1), expected_X);
assert_eq!(round_to_decimal_places(results.1, 1), expected_Y);
assert_eq!(round_to_decimal_places(results.2, 1), expected_Z);
assert_eq!(round_to_decimal_places(results.3, 1), expected_H);
assert_eq!(round_to_decimal_places(results.4, 1), expected_F);
assert_eq!(round_to_decimal_places(results.5, 2), expected_I);
assert_eq!(round_to_decimal_places(results.6, 2), expected_D);
assert_eq!(round_to_decimal_places(results.7, 1), expected_X_dot);
assert_eq!(round_to_decimal_places(results.8, 1), expected_Y_dot);
assert_eq!(round_to_decimal_places(results.9, 1), expected_Z_dot);
assert_eq!(round_to_decimal_places(results.10, 1), expected_H_dot);
assert_eq!(round_to_decimal_places(results.11, 1), expected_F_dot);
assert!(f64::abs(results.12 - expected_I_dot) <= 0.01);
assert!(f64::abs(results.13 - expected_D_dot) <= 0.01);
},
Err(e) => {
eprintln!("Error: {}", e);
process::exit(1);
}
}
}
#[test]
fn wmm2020_case1() {
test(2020.0, 0.0, 80.0, 0.0,
6570.4, -146.3, 54606.0, 6572.0, 55000.1, 83.14, -1.28, -1.28,
-16.2, 59.0, 42.9, -17.5, 40.5, 0.02, 0.51
);
}
#[test]
fn wmm2020_case2() {
test(2020.0, 0.0, 0.0, 120.0, 39624.3, 109.9, -10932.5, 39624.4, 41104.9, -15.42, 0.16, 020.0, 24.2, -60.8, 49.2, 24.0, 10.1, 0.08, -0.09);
}
#[test]
fn wmm2020_case3() {
test(2020.0, 0.0, -80.0, 240.0, 5940.6, 15772.1, -52480.8, 16853.8, 55120.6, -72.20, 69.36, -50.64, 30.4, 1.8, 91.7, 12.4, -83.5, 0.04, -0.10);
}
#[test]
fn wmm2020_case4() {
test(2020.0, 100.0, 80.0, 0.0, 6261.8, -185.5, 52429.1, 6264.5, 52802.0, 83.19, -1.70, -1.70, -15.1, 56.4, 39.2, -16.8, 36.9, 0.02, 0.51);
}
#[test]
fn wmm2020_case5() {
test(2020.0, 100.0, 0.0, 120.0, 37636.7, 104.9, -10474.8, 37636.9, 39067.3, -15.55, 0.16, 0.0, 22.9, -56.1, 45.1, 22.8, 9.8, 0.07, -0.09);
}
#[test]
fn wmm2020_case6() {
test(2020.0, 100.0, -80.0, 240.0, 5744.9, 14799.5, -49969.4, 15875.4, 52430.6, -72.37, 68.78, -51.22, 28.0, 1.4, 85.6, 11.4, -78.1, 0.04, -0.09);
}
#[test]
fn wmm2020_case7() {
test(2022.5, 0.0, 80.0, 0.0, 6529.9, 1.1, 54713.4, 6529.9, 55101.7, 83.19, 0.01, 0.01, -16.2, 59.0, 42.9, -16.2, 40.7, 0.02, 0.52);
}
#[test]
fn wmm2020_case8() {
test(2022.5, 0.0, 0.0, 120.0, 39684.7, -42.2, -10809.5, 39684.7, 41130.5, -15.24, -0.06, 0.0, 24.2, -60.8, 49.2, 24.2, 10.5, 0.08, -0.09);
}
#[test]
fn wmm2020_case9() {
test(2022.5, 0.0, -80.0, 240.0, 6016.5, 15776.7, -52251.6, 16885.0, 54912.1, -72.09, 69.13, -50.87, 30.4, 1.8, 91.7, 12.6, -83.4, 0.04, -0.09);
}
#[test]
fn wmm2020_case10() {
test(2022.5, 100.0, 80.0, 0.0, 6224.0, -44.5, 52527.0, 6224.2, 52894.5, 83.24, -0.41, -0.41, -15.1, 56.4, 39.2, -15.5, 37.1, 0.02, 0.52);
}
#[test]
fn wmm2020_case11() {
test(2022.5, 100.0, 0.0, 120.0, 37694.0, -35.3, -10362.0, 37694.1, 39092.4, -15.37, -0.05, 0.0, 22.9, -56.1, 45.1, 23.0, 10.2, 0.07, -0.09);
}
#[test]
fn wmm2020_case12() {
test(2022.5, 100.0, -80.0, 240.0, 5815.0, 14803.0, -49755.3, 15904.1, 52235.4, -72.27, 68.55, -51.47, 28.0, 1.4, 85.6, 11.6, -78.0, 0.04, -0.09);
}
}

77
src/main.rs Normal file
View file

@ -0,0 +1,77 @@
use geomag_wmm::initialise_magnetic_model;
fn main() {
// Read path from command-line first argument
let args: Vec<String> = std::env::args().collect();
if args.len() != 2 {
eprintln!("Usage: {} <path>", args[0]);
std::process::exit(1);
}
let path = &args[1];
let magnetic_model = initialise_magnetic_model(path);
// Collecting user input
println!("Height above WGS84 [metres]:");
let mut height_above_wgs84 = String::new();
std::io::stdin().read_line(&mut height_above_wgs84).expect("Failed to read line");
let height_above_wgs84: f64 = height_above_wgs84.trim().parse().expect("Please type a number!");
println!("Latitude [degrees]:");
let mut latitude = String::new();
std::io::stdin().read_line(&mut latitude).expect("Failed to read line");
let latitude: f64 = latitude.trim().parse().expect("Please type a number!");
println!("Longitude [degrees]:");
let mut longitude = String::new();
std::io::stdin().read_line(&mut longitude).expect("Failed to read line");
let longitude: f64 = longitude.trim().parse().expect("Please type a number!");
println!("Date [YYYY-MM-DD|YYYY.decimal]:");
let mut date = String::new();
let mut decimal_year = 0.0;
std::io::stdin().read_line(&mut date).expect("Failed to read line");
if date.contains('-') {
let date: Vec<&str> = date.trim().split("-").collect();
let year: f64 = date[0].parse().expect("Please type a number!");
let month: f64 = date[1].parse().expect("Please type a number!");
let day: f64 = date[2].parse().expect("Please type a number!");
decimal_year = year + (month - 1.0) / 12.0 + (day - 1.0) / 365.25;
}
else if date.contains('.') {
decimal_year = date.trim().parse().expect("Please type a number!");
}
else {
eprintln!("Please type a valid date format! (YYYY-MM-DD or YYYY.decimal)");
std::process::exit(1);
}
// Calculating
let result = magnetic_model.calculate(height_above_wgs84, latitude, longitude, decimal_year);
match result {
Ok(results) => {
println!("Results:");
println!("X: {:?} nT", results.0);
println!("Y: {:?} nT", results.1);
println!("Z: {:?} nT", results.2);
println!("H: {:?} nT", results.3);
println!("F: {:?} nT", results.4);
println!("I: {:?} degrees", results.5);
println!("D: {:?} degrees", results.6);
println!("X_dot: {:?} nT/year", results.7);
println!("Y_dot: {:?} nT/year", results.8);
println!("Z_dot: {:?} nT/year", results.9);
println!("H_dot: {:?} nT/year", results.10);
println!("F_dot: {:?} nT/year", results.11);
println!("I_dot: {:?} degrees/year", results.12);
println!("D_dot: {:?} degrees/year", results.13);
},
Err(e) => {
eprintln!("Error: {}", e);
std::process::exit(1);
}
}
}