calculate magnetic declination

Inspired by a recent thread, 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.

# Magnetic declination calculation based on WMM2015 earth magnetism mo
+del.
# See https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
use strict;
use warnings;

$, = " ";
$\ = "\n";
print magnetic_declination(240, -80, 100e3, 2017.5); # WMM sample data
print magnetic_declination(-122-4/60, 37+23/60, 32, 2017); # Mountain 
+View

# magnetic_declination($lon, $lat, $hgt, $yr)
# $lon: degrees longitude (east is positive)
# $lat: degrees latitude (north is positive)
# $hgt: elevation from sea level in meters, default=0
# $yr: year, default=2015
# Returns magnetic declination for the given location in degrees.
# In array context, also returns magnetic inclination (dip).
[download]

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);
}}
[download]

Comment on calculate magnetic declination Select or Download Code

Replies are listed 'Best First'.
Re: calculate magnetic declination by stevieb (Canon) on Jun 02, 2017 at 11:18 UTC
This is absolutely fantastic! I'm the author of the thread that inspired this, and it sure looks like it alleviates me of having to deal with this in the future :) I'll gladly adopt this and put it on the CPAN! Give me the weekend and I'll put it all together and post back here when it's complete. Thanks!	[reply]
Re^2: calculate magnetic declination by stevieb (Canon) on Jun 02, 2017 at 18:56 UTC
Extremely rough first cut after distribution-ing the code during a 10 minute break at work, here on my Github. I'll easily have a full test suite and it all ready for the CPAN by tomorrow. Thanks again!	[reply]
Re^3: calculate magnetic declination by no_slogan (Deacon) on Jun 02, 2017 at 19:14 UTC
Cool! I might suggest Geomag:: or Geomagnetism:: as a namespace. Thoughts, anyone?	[reply]
Re^4: calculate magnetic declination by stevieb (Canon) on Jun 02, 2017 at 19:54 UTC
Re^5: calculate magnetic declination by zentara (Archbishop) on Jun 02, 2017 at 20:03 UTC
Some notes below your chosen depth have not been shown here

Back to Cool Uses for Perl