Syntactic Confectionery Delight PerlMonks

### calculate magnetic declination

by no_slogan (Deacon)
 on Jun 02, 2017 at 00:53 UTC ( #1191907=CUFP: print w/replies, xml ) Need Help??

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

sub magnetic_declination {
my (\$lon, \$lat, \$hgt, \$yr) = @_;
\$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);
}}

Replies are listed 'Best First'.
Re: calculate magnetic declination
by stevieb (Abbot) 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!

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!

Cool! I might suggest Geomag:: or Geomagnetism:: as a namespace. Thoughts, anyone?

Create A New User
Node Status?
node history
Node Type: CUFP [id://1191907]
Approved by kcott
Front-paged by Corion
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (8)
As of 2019-05-27 14:16 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
Do you enjoy 3D movies?

Results (156 votes). Check out past polls.

Notices?