From 76e5e864040ec2c0270f1f80b28d7b5353adb899 Mon Sep 17 00:00:00 2001 From: Trevor Campbell Date: Wed, 11 Oct 2023 16:39:46 +1100 Subject: [PATCH] Convert Geomagnetism.java to rust --- src/geomagnetism.rs | 564 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 564 insertions(+) create mode 100644 src/geomagnetism.rs diff --git a/src/geomagnetism.rs b/src/geomagnetism.rs new file mode 100644 index 0000000..5a135ed --- /dev/null +++ b/src/geomagnetism.rs @@ -0,0 +1,564 @@ +// +// Module to calculate magnetic declination, magnetic field strength, +// inclination etc. for any point on the earth.

+//

Adapted from the geomagc software and World Magnetic Model of the NOAA +// Satellite and Information Service, National Geophysical Data Center

+// http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml +//

© Deep Pradhan, 2017

*/ +// Converted from original Java code + +use std::str::SplitWhitespace; + +use chrono::{Datelike, DateTime, Utc}; + +/** The input string array which contains each line of input for the wmm.cof input file. +* The columns in this file are as follows: n, m, gnm, hnm, dgnm, dhnm*/ +const WMM_COF: [&str; 91] = [" 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", ]; + +/** Mean radius of IAU-66 ellipsoid, in km*/ +const IAU66_RADIUS: f64 = 6371.2; + +/** Semi-major axis of WGS-1984 ellipsoid, in km*/ +const WGS84_A: f64 = 6378.137; + +/** Semi-minor axis of WGS-1984 ellipsoid, in km*/ +const WGS84_B: f64 = 6356.7523142; + +/** The maximum number of degrees of the spherical harmonic model*/ +const MAX_DEG: usize = 12; + +pub struct Geomagnetism { + /** Geomagnetic declination (decimal degrees) [opposite of variation, positive Eastward/negative Westward]*/ + declination: f64, + + /** Geomagnetic inclination/dip angle (degrees) [positive downward]*/ + inclination: f64, + + /** Geomagnetic field intensity/strength (nano Teslas)*/ + intensity: f64, + + /** Geomagnetic horizontal field intensity/strength (nano Teslas)*/ + bh: f64, + + /** Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward]*/ + bz: f64, + + /** Geomagnetic North South (northerly component) field intensity/strength (nano Tesla)*/ + bx: f64, + + /** Geomagnetic East West (easterly component) field intensity/strength (nano Teslas)*/ + by: f64, + + /** The maximum order of spherical harmonic model*/ + maxord: usize, + + /** The Gauss coefficients of main geomagnetic model (nt)*/ + c: [[f64; 13]; 13], + + /** The Gauss coefficients of secular geomagnetic model (nt/yr)*/ + cd: [[f64; 13]; 13], + + /** The time adjusted geomagnetic gauss coefficients (nt)*/ + tc: [[f64; 13]; 13], + + /** The theta derivative of p(n,m) (unnormalized)*/ + dp: [[f64; 13]; 13], + + /** The Schmidt normalization factors*/ + snorm: [f64; 169], + + /** The sine of (m*spherical coordinate longitude)*/ + sp: [f64; 13], + + /** The cosine of (m*spherical coordinate longitude)*/ + cp: [f64; 13], + + fn_: [f64; 13], + + fm: [f64; 13], + + /** The associated Legendre polynomials for m = 1 (unnormalized)*/ + pp: [f64; 13], + + k: [[f64; 13]; 13], + + /** The variables otime (old time), oalt (old altitude), + * olat (old latitude), olon (old longitude), are used to + * store the values used from the previous calculation to + * save on calculation time if some inputs don't change*/ + otime: f64, + oalt: f64, + olat: f64, + olon: f64, + + /** The date in years, for the start of the valid time of the fit coefficients*/ + epoch: f64, + + r: f64, + d: f64, + ca: f64, + sa: f64, + ct: f64, + st: f64, +} + +impl Geomagnetism { + /** Initialise the instance without calculations*/ + fn create() -> Geomagnetism { + let maxord = MAX_DEG; + let sp = [0.0; 13]; + let mut cp = [1.0; 13]; + cp[0] = 1.0; + let mut snorm = [0.0; 169]; + snorm[0] = 1.0; + let mut pp = [0.0; 13]; + pp[0] = 1.0; + let dp = [[0.0; 13]; 13]; + let mut c = [[0.0; 13]; 13]; + let mut cd = [[0.0; 13]; 13]; + + let mut k = [[0.0; 13]; 13]; + let mut fn_ = [0.0; 13]; + let mut fm = [0.0; 13]; + let tc = [[0.0; 13]; 13]; + + + let epoch = WMM_COF[0].trim().split_whitespace().next().unwrap().parse::().unwrap_or(0.0); + let mut tokens: SplitWhitespace; + let mut gnm: f64; + let mut hnm: f64; + let mut dgnm: f64; + let mut dhnm: f64; + { + let mut i: usize = 1; + let mut m: usize; + let mut n: usize; + while i < WMM_COF.len() { + { + tokens = WMM_COF[i].trim().split_whitespace(); + n = tokens.next().unwrap().parse::().unwrap_or(0); + m = tokens.next().unwrap().parse::().unwrap_or(0); + gnm = tokens.next().unwrap().parse::().unwrap_or(0.0); + hnm = tokens.next().unwrap().parse::().unwrap_or(0.0); + dgnm = tokens.next().unwrap().parse::().unwrap_or(0.0); + dhnm = tokens.next().unwrap().parse::().unwrap_or(0.0); + if m <= n { + c[m][n] = gnm; + cd[m][n] = dgnm; + if m != 0 { + c[n][m - 1] = hnm; + cd[n][m - 1] = dhnm; + } + } + } + i += 1; + } + } + + // Convert schmidt normalized gauss coefficients to unnormalized + snorm[0] = 1.0; + let mut flnmj: f64; + { + let mut j: usize; + let mut n: usize = 1; + while n <= maxord { + { + snorm[n] = snorm[n - 1] * (2.0 * n as f64 - 1.0) / n as f64; + j = 2; + { + let mut m: usize = 0; + let d1: usize = 1; + let mut d2: usize = (n - m + d1) / d1; + while d2 > 0 { + { + k[m][n] = (((n as i32 - 1) * (n as i32 - 1)) - (m * m) as i32) as f64 / ((2 * n as i32 - 1) * (2 * n as i32 - 3)) as f64; + if m > 0 { + flnmj = ((n - m + 1) * j) as f64 / (n + m) as f64; + snorm[n + m * 13] = snorm[n + (m - 1) * 13] * flnmj.sqrt(); + j = 1; + c[n][m - 1] = snorm[n + m * 13] * c[n][m - 1]; + cd[n][m - 1] = snorm[n + m * 13] * cd[n][m - 1]; + } + c[m][n] = snorm[n + m * 13] * c[m][n]; + cd[m][n] = snorm[n + m * 13] * cd[m][n]; + } + d2 -= 1; + m += d1; + } + } + + fn_ [n] = (n + 1) as f64; + fm[n] = n as f64; + } + n += 1; + } + } + + k[1][1] = 0.0; + fm[0] = 0.0; + Self { + declination: 0.0, + inclination: 0.0, + intensity: 0.0, + bh: 0.0, + bz: 0.0, + bx: 0.0, + by: 0.0, + maxord, + c, + cd, + tc, + dp, + snorm, + sp, + cp, + fn_, + fm, + pp, + k, + otime: -1000.0, + oalt: -1000.0, + olat: -1000.0, + olon: -1000.0, + epoch, + r: 0.0, + d: 0.0, + ca: 0.0, + sa: 0.0, + ct: 0.0, + st: 0.0, + } + } + + /** Initialise the instance and calculate for given location, altitude and date + * @param longitude Longitude in decimal degrees + * @param latitude Latitude in decimal degrees + * @param altitude Altitude in metres (with respect to WGS-1984 ellipsoid) + * @param calendar Calendar for date of calculation*/ + pub fn new(longitude: f64, latitude: f64, altitude: Option, calendar: Option>) -> Geomagnetism { + let mut geo = Self::create(); + geo.calculate(longitude, latitude, altitude, calendar); + geo + } + + /** Calculate for given location, altitude and date + * @param longitude Longitude in decimal degrees + * @param latitude Latitude in decimal degrees + * @param altitude Altitude in metres (with respect to WGS-1984 ellipsoid) + * @param calendar Calendar for date of calculation*/ + fn calculate(&mut self, longitude: f64, latitude: f64, altitude: Option, calendar: Option>) { + + let current_date = chrono::Utc::now(); + + let date = calendar.unwrap_or(current_date); + + let max_days = if date.naive_utc().date().leap_year() {366.0} else {365.0}; + + let altitude = altitude.unwrap_or(0.0); + let rlon: f64 = longitude.to_radians(); + let rlat: f64 = latitude.to_radians(); + let altitude_km: f64 = altitude / 1000.0; + let year_fraction: f64 = date.year() as f64 + date.ordinal() as f64 / max_days; + let dt: f64 = year_fraction - self.epoch; + let srlon: f64 = rlon.sin(); + let srlat: f64 = rlat.sin(); + let crlon: f64 = rlon.cos(); + let crlat: f64 = rlat.cos(); + let srlat2: f64 = srlat * srlat; + let crlat2: f64 = crlat * crlat; + let a2: f64 = WGS84_A * WGS84_A; + let b2: f64 = WGS84_B * WGS84_B; + let c2: f64 = a2 - b2; + let a4: f64 = a2 * a2; + let b4: f64 = b2 * b2; + let c4: f64 = a4 - b4; + self.sp[1] = srlon; + self.cp[1] = crlon; + // Convert from geodetic coords. to spherical coords. + if altitude_km != self.oalt || latitude != self.olat { + let q: f64 = (a2 - c2 * srlat2).sqrt(); + let q1: f64 = altitude_km * q; + let q2: f64 = ((q1 + a2) / (q1 + b2)) * ((q1 + a2) / (q1 + b2)); + let r2: f64 = (altitude_km * altitude_km) + 2.0 * q1 + (a4 - c4 * srlat2) / (q * q); + self.ct = srlat / (q2 * crlat2 + srlat2).sqrt(); + self.st = (1.0 - (self.ct * self.ct)).sqrt(); + self.r = r2.sqrt(); + self.d = (a2 * crlat2 + b2 * srlat2).sqrt(); + self.ca = (altitude_km + self.d) / self.r; + self.sa = c2 * crlat * srlat / (self.r * self.d); + } + if longitude != self.olon { + { + let mut m: usize = 2; + while m <= self.maxord { + { + self.sp[m] = self.sp[1] * self.cp[m - 1] + self.cp[1] * self.sp[m - 1]; + self.cp[m] = self.cp[1] * self.cp[m - 1] - self.sp[1] * self.sp[m - 1]; + } + m += 1; + } + } + } + let aor: f64 = IAU66_RADIUS / self.r; + let mut ar: f64 = aor * aor; + let mut br: f64 = 0.0; + let mut bt: f64 = 0.0; + let mut bp: f64 = 0.0; + let mut bpp: f64 = 0.0; + let mut par: f64; + let mut parp: f64; + let mut temp1: f64; + let mut temp2: f64; + { + let mut n: usize = 1; + while n <= self.maxord { + { + ar = ar * aor; + { + let mut m: usize = 0; + let d3: usize = 1; + let mut d4: usize = (n + m + d3) / d3; + while d4 > 0 { + { + // Compute unnormalized associated legendre polynomials and derivatives via recursion relations + if altitude_km != self.oalt || latitude != self.olat { + if n == m { + self.snorm[n + m * 13] = self.st * self.snorm[n - 1 + (m - 1) * 13]; + self.dp[m][n] = self.st * self.dp[m - 1][n - 1] + self.ct * self.snorm[n - 1 + (m - 1) * 13]; + } + if n == 1 && m == 0 { + self.snorm[n + m * 13] = self.ct * self.snorm[n - 1 + m * 13]; + self.dp[m][n] = self.ct * self.dp[m][n - 1] - self.st * self.snorm[n - 1 + m * 13]; + } + if n > 1 && n != m { + if m > n - 2 { + self.snorm[n - 2 + m * 13] = 0.0; + } + + if m > n - 2 { + self.dp[m][n - 2] = 0.0; + } + + self.snorm[n + m * 13] = self.ct * self.snorm[n - 1 + m * 13] - self.k[m][n] * self.snorm[n - 2 + m * 13]; + self.dp[m][n] = self.ct * self.dp[m][n - 1] - self.st * self.snorm[n - 1 + m * 13] - self.k[m][n] * self.dp[m][n - 2]; + } + } + // Time adjust the gauss coefficients + if year_fraction != self.otime { + self.tc[m][n] = self.c[m][n] + dt * self.cd[m][n]; + if m != 0 { + self.tc[n][m - 1] = self.c[n][m - 1] + dt * self.cd[n][m - 1]; + } + } + // Accumulate terms of the spherical harmonic expansions + par = ar * self.snorm[n + m * 13]; + if m == 0 { + temp1 = self.tc[m][n] * self.cp[m]; + temp2 = self.tc[m][n] * self.sp[m]; + } else { + temp1 = self.tc[m][n] * self.cp[m] + self.tc[n][m - 1] * self.sp[m]; + temp2 = self.tc[m][n] * self.sp[m] - self.tc[n][m - 1] * self.cp[m]; + } + bt = bt - ar * temp1 * self.dp[m][n]; + bp += self.fm[m] * temp2 * par; + br += self. fn_[n] * temp1 * par; + // Special case: north/south geographic poles + if self.st == 0.0 && m == 1 { + if n == 1 { + self.pp[n] = self.pp[n - 1]; + } else { + self.pp[n] = self.ct * self.pp[n - 1] - self.k[m][n] * self.pp[n - 2]; + } + + parp = ar * self.pp[n]; + bpp += self.fm[m] * temp2 * parp; + } + } + d4 -= 1; + m += d3; + } + } + } + n += 1; + } + } + + if self.st == 0.0 { + bp = bpp; + } else { + bp /= self.st; + } + + // Rotate magnetic vector components from spherical to geodetic coordinates + // bx must be the east-west field component + // by must be the north-south field component + // bz must be the vertical field component. + self.bx = -bt * self.ca - br * self.sa; + self.by = bp; + self.bz = bt * self.sa - br * self.ca; + // Compute declination (dec), inclination (dip) and total intensity (ti) + self.bh = ((self.bx * self.bx) + (self.by * self.by)).sqrt(); + self.intensity = ((self.bh * self.bh) + (self.bz * self.bz)).sqrt(); + // Calculate the declination. + self.declination = self.by.atan2(self.bx).to_degrees(); + self.inclination = self.bz.atan2(self.bh).to_degrees(); + self.otime = year_fraction; + self.oalt = altitude_km; + self.olat = latitude; + self.olon = longitude; + } + + /** @return Geomagnetic declination (degrees) [opposite of variation, positive Eastward/negative Westward]*/ + pub fn get_declination(&self) -> f64 { + return self.declination; + } + + /** @return Geomagnetic inclination/dip angle (degrees) [positive downward]*/ + pub fn get_inclination(&self) -> f64 { + return self.inclination; + } + + /** @return Geomagnetic field intensity/strength (nano Teslas)*/ + pub fn get_intensity(&self) -> f64 { + return self.intensity; + } + + /** @return Geomagnetic horizontal field intensity/strength (nano Teslas)*/ + pub fn get_horizontal_intensity(&self) -> f64 { + return self.bh; + } + + /** @return Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward]*/ + pub fn get_vertical_intensity(&self) -> f64 { + return self.bz; + } + + /** @return Geomagnetic North South (northerly component) field intensity/strength (nano Tesla)*/ + pub fn get_north_intensity(&self) -> f64 { + return self.bx; + } + + /** @return Geomagnetic East West (easterly component) field intensity/strength (nano Teslas)*/ + pub fn get_east_intensity(&self) -> f64 { + return self.by; + } +} + +#[cfg(test)] +mod tests { + use super::Geomagnetism; + + #[test] + fn test_declination() { + let geo = Geomagnetism::new(150.9231, -34.1356, None, None); + assert_between(geo.declination, 12.7, 12.8); + let geo = Geomagnetism::new(0.5, 51.48, None, None); + assert_between(geo.declination, 0.5, 1.0); + let geo = Geomagnetism::new(139.76, 35.55, None, None); + assert_between(geo.declination, -8.0, -7.5); + let geo = Geomagnetism::new(-16.65, 13.24, None, None); + assert_between(geo.declination, -7.0, -6.0); + let geo = Geomagnetism::new(-73.8, 40.6, None, None); + assert_between(geo.declination, -13.0, -12.0); + } + + fn assert_between(variable: f64, bottom: f64, top: f64) -> bool { + let result = variable >= bottom && variable <= top; + if !result { + assert!(result, "Variable {} not between {} and {}", variable, bottom, top); + } + result + } + +} +