igrf (empty) → 0.2.0.0
raw patch · 7 files changed
+283/−0 lines, 7 filesdep +addep +basedep +polynomialsetup-changed
Dependencies added: ad, base, polynomial
Files
- LICENSE +30/−0
- README.md +4/−0
- Setup.hs +2/−0
- igrf.cabal +37/−0
- src/IGRF.hs +70/−0
- src/Math/SphericalHarmonics.hs +101/−0
- src/Math/SphericalHarmonics/AssociatedLegendre.hs +39/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, J. Douglas McClean + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of J. Douglas McClean nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,4 @@+igrf +==== + +International Geomagnetic Reference Field
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple +main = defaultMain
+ igrf.cabal view
@@ -0,0 +1,37 @@+-- Initial igrf.cabal generated by cabal init. For further documentation, +-- see http://haskell.org/cabal/users-guide/ + +name: igrf +version: 0.2.0.0 +synopsis: International Geomagnetic Reference Field +description: + An implemetation of the Internation Geomagnetic Reference Field, a model of the Earth's magnetic field. + . + Includes model coefficients from the 11th Edition of the IGRF. +homepage: https://github.com/dmcclean/igrf +license: BSD3 +license-file: LICENSE +author: J. Douglas McClean +maintainer: douglas.mcclean@gmail.com +-- copyright: +category: Science +stability: Experimental +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 + +source-repository head + type: git + location: git://github.com/dmcclean/igrf.git + +library + exposed-modules: IGRF, + Math.SphericalHarmonics, + Math.SphericalHarmonics.AssociatedLegendre + -- other-modules: + -- other-extensions: + build-depends: base >=4.7 && <4.8, + ad >= 4.2, + polynomial >= 0.7.1 + hs-source-dirs: src + default-language: Haskell2010
+ src/IGRF.hs view
@@ -0,0 +1,70 @@+-- | An implementation of the International Geomagnetic Reference Field, as defined at <http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html>. +module IGRF +( + MagneticModel(..) +, igrf11 +, fieldAtTime +, evaluateModelGradientInLocalTangentPlane +) +where + +import Math.SphericalHarmonics + +-- | Represents a spherical harmonic model of a magnetic field. +data MagneticModel a = MagneticModel + { + fieldAtEpoch :: SphericalHarmonicModel a -- ^ Field at model epoch in nT, reference radius in km + , secularVariation :: SphericalHarmonicModel a -- ^ Secular variation in nT / yr, reference radius in km + } + +-- | Gets a spherical harmonic model of a magnetic field at a specified time offset from the model epoch. +fieldAtTime :: (Num a, Eq a) => MagneticModel a -- ^ Magnetic field model + -> a -- ^ Time since model epoch (year) + -> SphericalHarmonicModel a -- ^ Spherical harmonic model of magnetic field at specified time. Field in nT, reference radius in km +fieldAtTime m t = combine (fieldAtEpoch m) (scale t $ secularVariation m) + +-- | The International Geomagnetic Reference Field model, 11th edition. +-- Model epoch is January 1st, 2010. +igrf11 :: (Floating a) => MagneticModel a +igrf11 = MagneticModel + { + fieldAtEpoch = f + , secularVariation = s + } + where + f = SphericalHarmonicModel + { + modelDegree = 13 + , referenceRadius = r + , coefficients = [(0, 0), + (-29496.5, 0), (-1585.9, 4945.1), + (-2396.6, 0), (3026.0, -2707.7), (1668.6, -575.4), + (1339.7, 0), (-2326.3, -160.5), (1231.7, 251.7), (634.2, -536.8), + (912.6, 0), (809.0, 286.5), (166.6, -211.2), (-357.1, 164.4), (89.7, -309.2), + (-231.1, 0), (357.2, 44.7), (200.3, 188.9), (-141.2, -118.1), (-163.1, 0.1), (-7.7, 100.9), + (72.8, 0), (68.6, -20.8), (76.0, 44.2), (-141.4, 61.5), (-22.9, -66.3), (13.1, 3.1), (-77.9, 54.9), + (80.4, 0), (-75.0, -57.8), (-4.7, -21.2), (45.3, 6.6), (14.0, 24.9), (10.4, 7.0), (1.6, -27.7), (4.9, -3.4), + (24.3, 0), (8.2, 10.9), (-14.5, -20.0), (-5.7, 11.9), (-19.3, -17.4), (11.6, 16.7), (10.9, 7.1), (-14.1, -10.8), (-3.7, 1.7), + (5.4, 0), (9.4, -20.5), (3.4, 11.6), (-5.3, 12.8), (3.1, -7.2), (-12.4, -7.4), (-0.8, 8.0), (8.4, 2.2), (-8.4, -6.1), (-10.1, 7.0), + (-2.0, 0), (-6.3, 2.8), (0.9, -0.1), (-1.1, 4.7), (-0.2, 4.4), (2.5, -7.2), (-0.3, -1.0), (2.2, -4.0), (3.1, -2.0), (-1.0, -2.0), (-2.8, -8.3), + (3.0, 0), (-1.5, 0.1), (-2.1, 1.7), (1.6, -0.6), (-0.5, -1.8), (0.5, 0.9), (-0.8, -0.4), (0.4, -2.5), (1.8, -1.3), (0.2, -2.1), (0.8, -1.9), (3.8, -1.8), + (-2.1, 0), (-0.2, -0.8), (0.3, 0.3), (1.0, 2.2), (-0.7, -2.5), (0.9, 0.5), (-0.1, 0.6), (0.5, 0.0), (-0.4, 0.1), (-0.4, 0.3), (0.2, -0.9), (-0.8, -0.2), (0.0, 0.8), + (-0.2, 0), (-0.9, -0.8), (0.3, 0.3), (0.4, 1.7), (-0.4, -0.6), (1.1, -1.2), (-0.3, -0.1), (0.8, 0.5), (-0.2, 0.1), (0.4, 0.5), (0.0, 0.4), (0.4, -0.2), (-0.3, -0.5), (-0.3, -0.8) + ] + } + s = SphericalHarmonicModel + { + modelDegree = 8 + , referenceRadius = r + , coefficients = [(0, 0), + (11.4, 0), (16.7, -28.8), + (-11.3, 0), (-3.9, -23.0), (2.7, -12.9), + (1.3, 0), (-3.9, 8.6), (-2.9, -2.9), (-8.1, -2.1), + (-1.4, 0), (2.0, 0.4), (-8.9, 3.2), (4.4, 3.6), (-2.3, -0.8), + (-0.5, 0), (0.5, 0.5), (-1.5, 1.5), (-0.7, 0.9), (1.3, 3.7), (1.4, -0.6), + (-0.3, 0), (-0.3, -0.1), (-0.3, -2.1), (1.9, -0.4), (-1.6, -0.5), (-0.2, 0.8), (1.8, 0.5), + (0.2, 0), (-0.1, 0.6), (-0.6, 0.3), (1.4, -0.2), (0.3, -0.1), (0.1, -0.8), (-0.8, -0.3), (0.4, 0.2), + (-0.1, 0), (0.1, 0.0), (-0.5, 0.2), (0.3, 0.5), (-0.3, 0.4), (0.3, 0.1), (0.2, -0.1), (-0.5, 0.4), (0.2, 0.4) + ] + } + r = 6371.2
+ src/Math/SphericalHarmonics.hs view
@@ -0,0 +1,101 @@+{-# LANGUAGE DeriveFunctor #-} + +-- | Provides spherical harmonic models of scalar-valued functions. +module Math.SphericalHarmonics +( + SphericalHarmonicModel(..) +, combine +, scale +, evaluateModel +, evaluateModelGradient +, evaluateModelGradientInLocalTangentPlane +) +where + +import Math.SphericalHarmonics.AssociatedLegendre +import Numeric.AD + +-- | Represents a spherical harmonic model of a scalar-valued function. +data SphericalHarmonicModel a = SphericalHarmonicModel + { + modelDegree :: Int -- ^ The maximum degree of the model. Must be >= 0. + , referenceRadius :: a -- ^ The reference radius used to define the model. + , coefficients :: [(a, a)] -- ^ G and H coefficients of the model and their secular variations. + -- These coefficients are stored in the order [(g_0_0, h_0_0), (g_1_0, h1_0_), 1_1, 2_0, 2_1, 2_2, 3_0, 3_1, 3_2, 3_3, ...] + -- There must be Triangle('modelDegree' + 1) coefficients. + } + deriving (Functor) + +-- TODO: consider how to relax the reference radius error condition +-- TODO: make SphericalHarmonicModel an instance of additive typeclass +-- | Adds two compatible spherical harmonic models. +combine :: (Num a, Eq a) => SphericalHarmonicModel a -> SphericalHarmonicModel a -> SphericalHarmonicModel a +combine m1 m2 | (referenceRadius m1 /= referenceRadius m2) = error "Incompatible model reference radii." + | otherwise = SphericalHarmonicModel + { + modelDegree = max (modelDegree m1) (modelDegree m2) + , referenceRadius = referenceRadius m1 + , coefficients = combineCoefficients (coefficients m1) (coefficients m2) + } + where + combineCoefficients c1 c2 = take (max (length c1) (length c2)) $ zipWith addPairs (c1 ++ repeat (0,0)) (c2 ++ repeat (0,0)) + addPairs (g1, h1) (g2, h2) = (g1 + g2, h1 + h2) + +-- | Linearly scales a spherical harmonic model. +scale :: (Num a) => a -> SphericalHarmonicModel a -> SphericalHarmonicModel a +scale x m = m { coefficients = fmap scalePair (coefficients m) } + where + scalePair (g, h) = (x * g, x * h) + +-- | Computes the scalar value of the spherical harmonic model at a specified spherical position. +evaluateModel :: (Floating a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model + -> a -- ^ Spherical radius + -> a -- ^ Spherical colatitude (radian) + -> a -- ^ Spherical longitude (radian) + -> a -- ^ Model value +evaluateModel model r colat lon = refR * sumOverDegree + where + refR = referenceRadius model + deg = modelDegree model + gs = map fst $ coefficients model + hs = map snd $ coefficients model + sumOverDegree = sum $ fmap degreeTerm [0..deg] + degreeTerm n = ((refR / r) ^ (n + 1)) * (sum $ fmap (orderTerm n) [0..n]) + orderTerm n m = lonFactor * (p (cos colat)) + where + scaledLon = lon * fromIntegral m + lonFactor = (g * cos scaledLon) + (h * sin scaledLon) + p = schmidtSemiNormalizedAssociatedLegendreFunction n m + g = gs !! computeIndex n m + h = hs !! computeIndex n m + +-- | Computes the gradient of the scalar value of the spherical harmonic model, in spherical coordinates, at a specified location. +evaluateModelGradient :: (Floating a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model + -> a -- ^ Spherical radius + -> a -- ^ Spherical colatitude (radian) + -> a -- ^ Spherical longitude (radian) + -> (a, a, a) -- ^ Radial, colatitudinal, and longitudinal components of gradient +evaluateModelGradient model r colat lon = makeTuple . fmap negate $ modelGrad [r, colat, lon] + where + modelGrad = grad (\[r', c', l'] -> evaluateModel (fmap auto model) r' c' l') + makeTuple [x, y, z] = (x, y, z) + +-- | Computes the gradient of the scalar value of the spherical harmonic model at a specified location, in Cartesian coordinates. +-- The result is expressed in a reference frame locally tangent to the specified location. +evaluateModelGradientInLocalTangentPlane :: (Floating a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model + -> a -- ^ Spherical radius + -> a -- ^ Spherical colatitude (radian) + -> a -- ^ Spherical longitude (radian) + -> (a, a, a) -- ^ East, North, and up components of gradient +evaluateModelGradientInLocalTangentPlane model r colat lon = (e, n, u) + where + (r', colat', lon') = evaluateModelGradient model r colat lon + e = lon' / (r * sin colat) + n = -colat' / r -- negated because the colatitude increase southward + u = r' + +computeIndex :: Int -> Int -> Int +computeIndex n m = triangle n + m + +triangle :: Int -> Int +triangle n = (n * (n + 1)) `div` 2
+ src/Math/SphericalHarmonics/AssociatedLegendre.hs view
@@ -0,0 +1,39 @@+-- | Provides definitions of associated Legendre functions used in spherical harmonic models. +module Math.SphericalHarmonics.AssociatedLegendre +( + associatedLegendreFunction +, schmidtSemiNormalizedAssociatedLegendreFunction +) +where + +import Math.Polynomial hiding (x) +import Math.Polynomial.Legendre + +-- definition from http://www.mathworks.com/help/matlab/ref/legendre.html#f89-998354 +-- | Computes the associated Legendre function of degree 'n' and order 'm'. +-- Note that the resulting function may not be a polynomial, as when `m` is odd it involves a fractional power of `x`. +-- As used in the geodesy and magnetics literature, these functions do not include the Condon-Shortley phase. +associatedLegendreFunction :: (Floating a, Ord a) => Int -- ^ Degree 'n' of the desired associated Legendre function. + -> Int -- ^ Order 'm' of the desired associated Legendre function. + -> a -> a +associatedLegendreFunction n m = f + where + f x = (nonPolyTerm x) * (evalPoly p' x) + nonPolyTerm x = (1 - (x * x)) ** (fromIntegral m / 2) + p' = polyDerivs p !! m + p = legendre n + +-- definition from http://www.mathworks.com/help/matlab/ref/legendre.html#f89-998354 +-- | Computes the Schmidt semi-normalized associated Legendre function of degree 'n' and order 'm'. +-- As used in the geodesy and magnetics literature, these functions do not include the Condon-Shortley phase. +schmidtSemiNormalizedAssociatedLegendreFunction :: (Floating a, Ord a) => Int -- ^ Degree 'n' of the desired function. + -> Int -- ^ Order 'm' of the desired function. + -> a -> a +schmidtSemiNormalizedAssociatedLegendreFunction n 0 = associatedLegendreFunction n 0 +schmidtSemiNormalizedAssociatedLegendreFunction n m = (* factor) . associatedLegendreFunction n m + where + factor = (sqrt $ 2 / rawFactor) + rawFactor = fromIntegral $ rawFactor' (fromIntegral n) (fromIntegral m) + +rawFactor' :: Integer -> Integer -> Integer +rawFactor' n m = product . map (max 1) $ enumFromTo (n - m + 1) (n + m)