, calculate the magnetic declination (angle between true north and magnetic north) for a given longitude/latitude. If someone wants to put this on CPAN, be my guest, but then you're on the hook to update it to use the WMM2020 model when it's released in late 2019.
BEGIN {
my @WMM2015 = (
[],
[ [-29438.5, 0, 10.7, 0],
[-1501.1, 4796.2, 17.9, -26.8] ],
[ [-2445.3, 0, -8.6, 0],
[3012.5, -2845.6, -3.3, -27.1],
[1676.6, -642, 2.4, -13.3] ],
[ [1351.1, 0, 3.1, 0],
[-2352.3, -115.3, -6.2, 8.4],
[1225.6, 245, -0.4, -0.4],
[581.9, -538.3, -10.4, 2.3] ],
[ [907.2, 0, -0.4, 0],
[813.7, 283.4, 0.8, -0.6],
[120.3, -188.6, -9.2, 5.3],
[-335, 180.9, 4, 3],
[70.3, -329.5, -4.2, -5.3] ],
[ [-232.6, 0, -0.2, 0],
[360.1, 47.4, 0.1, 0.4],
[192.4, 196.9, -1.4, 1.6],
[-141, -119.4, 0, -1.1],
[-157.4, 16.1, 1.3, 3.3],
[4.3, 100.1, 3.8, 0.1] ],
[ [69.5, 0, -0.5, 0],
[67.4, -20.7, -0.2, 0],
[72.8, 33.2, -0.6, -2.2],
[-129.8, 58.8, 2.4, -0.7],
[-29, -66.5, -1.1, 0.1],
[13.2, 7.3, 0.3, 1],
[-70.9, 62.5, 1.5, 1.3] ],
[ [81.6, 0, 0.2, 0],
[-76.1, -54.1, -0.2, 0.7],
[-6.8, -19.4, -0.4, 0.5],
[51.9, 5.6, 1.3, -0.2],
[15, 24.4, 0.2, -0.1],
[9.3, 3.3, -0.4, -0.7],
[-2.8, -27.5, -0.9, 0.1],
[6.7, -2.3, 0.3, 0.1] ],
[ [24, 0, 0, 0],
[8.6, 10.2, 0.1, -0.3],
[-16.9, -18.1, -0.5, 0.3],
[-3.2, 13.2, 0.5, 0.3],
[-20.6, -14.6, -0.2, 0.6],
[13.3, 16.2, 0.4, -0.1],
[11.7, 5.7, 0.2, -0.2],
[-16, -9.1, -0.4, 0.3],
[-2, 2.2, 0.3, 0] ],
[ [5.4, 0, 0, 0],
[8.8, -21.6, -0.1, -0.2],
[3.1, 10.8, -0.1, -0.1],
[-3.1, 11.7, 0.4, -0.2],
[0.6, -6.8, -0.5, 0.1],
[-13.3, -6.9, -0.2, 0.1],
[-0.1, 7.8, 0.1, 0],
[8.7, 1, 0, -0.2],
[-9.1, -3.9, -0.2, 0.4],
[-10.5, 8.5, -0.1, 0.3] ],
[ [-1.9, 0, 0, 0],
[-6.5, 3.3, 0, 0.1],
[0.2, -0.3, -0.1, -0.1],
[0.6, 4.6, 0.3, 0],
[-0.6, 4.4, -0.1, 0],
[1.7, -7.9, -0.1, -0.2],
[-0.7, -0.6, -0.1, 0.1],
[2.1, -4.1, 0, -0.1],
[2.3, -2.8, -0.2, -0.2],
[-1.8, -1.1, -0.1, 0.1],
[-3.6, -8.7, -0.2, -0.1] ],
[ [3.1, 0, 0, 0],
[-1.5, -0.1, 0, 0],
[-2.3, 2.1, -0.1, 0.1],
[2.1, -0.7, 0.1, 0],
[-0.9, -1.1, 0, 0.1],
[0.6, 0.7, 0, 0],
[-0.7, -0.2, 0, 0],
[0.2, -2.1, 0, 0.1],
[1.7, -1.5, 0, 0],
[-0.2, -2.5, 0, -0.1],
[0.4, -2, -0.1, 0],
[3.5, -2.3, -0.1, -0.1] ],
[ [-2, 0, 0.1, 0],
[-0.3, -1, 0, 0],
[0.4, 0.5, 0, 0],
[1.3, 1.8, 0.1, -0.1],
[-0.9, -2.2, -0.1, 0],
[0.9, 0.3, 0, 0],
[0.1, 0.7, 0.1, 0],
[0.5, -0.1, 0, 0],
[-0.4, 0.3, 0, 0],
[-0.4, 0.2, 0, 0],
[0.2, -0.9, 0, 0],
[-0.9, -0.2, 0, 0],
[0, 0.7, 0, 0] ],
);
my $DEG2RAD = atan2(1,1)/45;
sub magnetic_declination {
my ($lon, $lat, $hgt, $yr) = @_;
$lon *= $DEG2RAD;
$lat *= $DEG2RAD;
$hgt //= 0;
$yr //= 2015;
warn "Model is valid only from 2015 to 2020" if $yr < 2015 || $yr >
+ 2020;
my ($geo_r, $geo_lat) = do { # geocentric coordinates
my $A = 6378137; # reference ellipsoid semimajor axis
my $f = 1 / 298.257223563; # flattening
my $e2 = $f * (2 - $f); # eccentricity
my $Rc = $A / sqrt(1 - $e2 * sin($lat)**2); # radius of curvatur
+e
my $p = ($Rc + $hgt) * cos($lat); # radius in x-y plane
my $z = ($Rc * (1 - $e2) + $hgt) * sin($lat);
(sqrt($p*$p + $z*$z), atan2($z, $p))
};
my $s = sin($geo_lat);
my $c = cos($geo_lat);
# associated Legendre polynomials (P) and derivatives (dP)
my @P = ([1],[$s,$c]);
my @dP = ([0],[$c,-$s]);
for my $n (2 .. $#WMM2015) {
my $k = 2*$n-1;
for my $m (0 .. $n-2) {
my $k1 = $k / ($n - $m);
my $k2 = ($n + $m - 1) / ($n - $m);
$P[$n][$m] = $k1 * $s * $P[$n-1][$m] - $k2 * $P[$n-2][$m];
$dP[$n][$m] = $k1 * ($s * $dP[$n-1][$m] + $c * $P[$n-1][$m])
- $k2 * $dP[$n-2][$m];
}
my $y = $k * $P[$n-1][$n-1];
my $dy = $k * $dP[$n-1][$n-1];
$P[$n][$n-1] = $s * $y;
$dP[$n][$n-1] = $s * $dy + $c * $y;
$P[$n][$n] = $c * $y;
$dP[$n][$n] = $c * $dy - $s * $y;
}
# Schmidt quasi-normalization
for my $n (1 .. $#WMM2015) {
my $f = sqrt(2);
for my $m (1 .. $n) {
$f /= sqrt(($n-$m+1) * ($n+$m));
$P[$n][$m] *= $f;
$dP[$n][$m] *= $f;
}
}
my $X = 0; # magnetic field north component in nT
my $Y = 0; # east component
my $Z = 0; # vertical component
my $t = $yr - 2015;
my $r = 6371200 / $geo_r; # radius relative to geomagnetic referenc
+e
my $R = $r * $r;
my @c = map cos($_*$lon), 0 .. $#WMM2015;
my @s = map sin($_*$lon), 0 .. $#WMM2015;
for my $n (1 .. $#WMM2015) {
my $x = my $y = my $z = 0;
for my $m (0 .. $n) {
my $row = $WMM2015[$n][$m];
my $g = $row->[0] + $t * $row->[2];
my $h = $row->[1] + $t * $row->[3];
$x += ($g * $c[$m] + $h * $s[$m]) * $dP[$n][$m];
$y += ($g * $s[$m] - $h * $c[$m]) * $P[$n][$m] * $m;
$z += ($g * $c[$m] + $h * $s[$m]) * $P[$n][$m];
}
$R *= $r;
$X -= $x * $R;
$Y += $y * $R;
$Z -= $z * $R * ($n+1);
}
$Y /= $c;
$c = cos($geo_lat - $lat); # transform back to geodetic coords
$s = sin($geo_lat - $lat);
($X, $Z) = ($X * $c - $Z * $s, $X * $s + $Z * $c);
my $decl = atan2($Y, $X) / $DEG2RAD;
return $decl unless wantarray;
my $incl = atan2($Z, sqrt($X*$X + $Y*$Y)) / $DEG2RAD;
return ($decl, $incl);
}}