, 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);
}}
`

Comment oncalculate magnetic declinationSelectorDownloadCode