ferromagnetic/igrf/
mod.rs

1use crate::{MagneticComponents, OrthogonalStrength};
2
3mod coeffs;
4mod math;
5
6pub struct IGRFresults {
7    pub result: MagneticComponents,
8    // Annual changes
9    pub sv: MagneticComponents,
10}
11
12pub struct IGRF {
13    coeffs: coeffs::IGRFCoeffs,
14}
15impl Default for IGRF {
16    fn default() -> IGRF {
17        IGRF {
18            coeffs: coeffs::igrf_data(),
19        }
20    }
21}
22impl IGRF {
23    pub fn calc(&self, lat: f64, lon: f64, alt: f64, date: f64) -> IGRFresults {
24        //input validation
25        let (start_coeffs, end_coeffs, nmax) = self.coeffs.coeffs(date);
26        let (a, b) = math::shval3(lat, lon, alt, nmax as usize, &start_coeffs, &end_coeffs);
27        let dif_a = math::Difh::from_orthognal_strength(&a);
28        let dif_b = math::Difh::from_orthognal_strength(&b);
29
30        let result = MagneticComponents {
31            declination: dif_a.declination.to_degrees(),
32            inclination: dif_a.inclination.to_degrees(),
33            horizontal_intensity: dif_a.horizontal_intensity,
34            orthogonal_strength: a,
35            total_intensity: dif_a.total_intensity,
36        };
37
38        let sv = MagneticComponents {
39            declination: {
40                let mut ddot = (dif_b.declination - dif_a.declination).to_degrees();
41                if ddot > 180.0 {
42                    ddot -= 360.0
43                }
44                if ddot <= -180.0 {
45                    ddot += 360.0
46                }
47                ddot * 60.0
48            },
49            inclination: (dif_b.inclination - dif_a.inclination).to_degrees() * 60.,
50            horizontal_intensity: dif_b.horizontal_intensity - dif_a.horizontal_intensity,
51            orthogonal_strength: OrthogonalStrength {
52                north: b.north - result.orthogonal_strength.north,
53                east: b.east - result.orthogonal_strength.east,
54                down: b.down - result.orthogonal_strength.down,
55            },
56            total_intensity: dif_b.total_intensity - dif_a.total_intensity,
57        };
58
59        IGRFresults { result, sv }
60    }
61}