commit 7671a9abb567145928babd120797962d2b85ca57 Author: Louis Date: Sat Sep 28 16:55:51 2024 +0200 working model diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ea8c4bf --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +/target diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..6c25a36 --- /dev/null +++ b/Cargo.lock @@ -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" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..6a740fb --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,6 @@ +[package] +name = "geomag-wmm" +version = "0.1.0" +edition = "2021" + +[dependencies] diff --git a/WMM.2020.COF b/WMM.2020.COF new file mode 100644 index 0000000..7940efb --- /dev/null +++ b/WMM.2020.COF @@ -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 diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..ed7fc79 --- /dev/null +++ b/src/lib.rs @@ -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::().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); + } +} \ No newline at end of file diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..b6723f4 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,77 @@ +use geomag_wmm::initialise_magnetic_model; + +fn main() { + // Read path from command-line first argument + let args: Vec = std::env::args().collect(); + if args.len() != 2 { + eprintln!("Usage: {} ", 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); + } + } +} \ No newline at end of file