geomag_wmm/

lib.rs

#![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);
    }
}