Beefy Boxes and Bandwidth Generously Provided by pair Networks
good chemistry is complicated,
and a little bit messy -LW
 
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] ], ); 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); }}

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?

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: CUFP [id://1191907]
Approved by kcott
Front-paged by Corion
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (9)
As of 2018-06-24 17:44 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    Should cpanminus be part of the standard Perl release?



    Results (126 votes). Check out past polls.

    Notices?