igrf 0.2.0.0 → 0.4.0.0
raw patch · 6 files changed
+266/−102 lines, 6 filesdep +polydep +semiringsdep +textdep −polynomialdep ~base
Dependencies added: poly, semirings, text, vector-space
Dependencies removed: polynomial
Dependency ranges changed: base
Files
- README.md +2/−0
- igrf.cabal +9/−5
- src/IGRF.hs +108/−37
- src/IGRF/Parser.hs +45/−0
- src/Math/SphericalHarmonics.hs +91/−51
- src/Math/SphericalHarmonics/AssociatedLegendre.hs +11/−9
README.md view
@@ -2,3 +2,5 @@ ==== International Geomagnetic Reference Field + +[](http://hackage.haskell.org/package/igrf)
igrf.cabal view
@@ -2,10 +2,10 @@ -- see http://haskell.org/cabal/users-guide/ name: igrf -version: 0.2.0.0 +version: 0.4.0.0 synopsis: International Geomagnetic Reference Field description: - An implemetation of the Internation Geomagnetic Reference Field, a model of the Earth's magnetic field. + An implemetation of the International 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 @@ -26,12 +26,16 @@ library exposed-modules: IGRF, + IGRF.Parser, Math.SphericalHarmonics, Math.SphericalHarmonics.AssociatedLegendre -- other-modules: -- other-extensions: - build-depends: base >=4.7 && <4.8, + build-depends: base >=4.7 && <5, ad >= 4.2, - polynomial >= 0.7.1 + poly >= 0.4, + semirings >= 0.5, + text >= 1.2, + vector-space >= 0.8.7 hs-source-dirs: src - default-language: Haskell2010+ default-language: Haskell2010
src/IGRF.hs view
@@ -3,68 +3,139 @@ ( MagneticModel(..) , igrf11 +, igrf12 +, igrf13 , fieldAtTime , evaluateModelGradientInLocalTangentPlane ) where +import Data.VectorSpace import Math.SphericalHarmonics -- | Represents a spherical harmonic model of a magnetic field. -data MagneticModel a = MagneticModel +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 +fieldAtTime :: (Fractional 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) +fieldAtTime m t = (fieldAtEpoch m) ^+^ (t *^ secularVariation m) -- | The International Geomagnetic Reference Field model, 11th edition. -- Model epoch is January 1st, 2010. -igrf11 :: (Floating a) => MagneticModel a +igrf11 :: (Fractional 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) + f = scaledSphericalHarmonicModel r fcs + fcs = [[(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) + s = scaledSphericalHarmonicModel r scs + scs = [[(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 + +-- | The International Geomagnetic Reference Field model, 12th edition. +-- Model epoch is January 1st, 2015. +igrf12 :: (Fractional a) => MagneticModel a +igrf12 = MagneticModel + { + fieldAtEpoch = f, + secularVariation = s + } + where + f = scaledSphericalHarmonicModel r fcs + fcs = [[(0, 0)], + [(-29442.0, 0), (-1501.0, 4797.1)], + [(-2445.1, 0), (3012.9, -2845.6), (1676.7, -641.9)], + [(1350.7, 0), (-2352.3, -115.3), (1225.6, 244.9), (582.0, -538.4)], + [(907.6, 0), (813.7, 283.3), (120.4, -188.7), (-334.9, 180.9), (70.4, -329.5)], + [(-232.6, 0), (360.1, 47.3), (192.4, 197.0), (-140.9, -119.3), (-157.5, 16.0), (4.1, 100.2)], + [(70.0, 0), (67.7, -20.8), (72.7, 33.2), (-129.9, 58.9), (-28.9, -66.7), (13.2, 7.3), (-70.9, 62.6)], + [(81.6, 0), (-76.1, -54.1), (-6.8, -19.5), (51.8, 5.7), (15.0, 24.4), (9.4, 3.4), (-2.8, -27.4), (6.8, -2.2)], + [(24.2, 0), (8.8, 10.1), (-16.9, -18.3), (-3.2, 13.3), (-20.6, -14.6), (13.4, 16.2), (11.7, 5.7), (-15.9, -9.1), (-2.0, 2.1)], + [(5.4, 0), (8.8, -21.6), (3.1, 10.8), (-3.3, 11.8), (0.7, -6.8), (-13.3, -6.9), (-0.1, 7.8), (8.7, 1.0), (-9.1, -4.0), (-10.5, 8.4)], + [(-1.9, 0), (-6.3, 3.2), (0.1, -0.4), (0.5, 4.6), (-0.5, 4.4), (1.8, -7.9), (-0.7, -0.6), (2.1, -4.2), (2.4, -2.8), (-1.8, -1.2), (-3.6, -8.7)], + [(3.1, 0), (-1.5, -0.1), (-2.3, 2.0), (2.0, -0.7), (-0.8, -1.1), (0.6, 0.8), (-0.7, -0.2), (0.2, -2.2), (1.7, -1.4), (-0.2, -2.5), (0.4, -2.0), (3.5, -2.4)], + [(-1.9, 0), (-0.2, -1.1), (0.4, 0.4), (1.2, 1.9), (-0.8, -2.2), (0.9, 0.3), (0.1, 0.7), (0.5, -0.1), (-0.3, 0.3), (-0.4, 0.2), (0.2, -0.9), (-0.9, -0.1), (0.0, 0.7)], + [(0.0, 0), (-0.9, -0.9), (0.4, 0.4), (0.5, 1.6), (-0.5, -0.5), (1.0, -1.2), (-0.2, -0.1), (0.8, 0.4), (-0.1, -0.1), (0.3, 0.4), (0.1, 0.5), (0.5, -0.3), (-0.4, -0.4), (-0.3, -0.8)] + ] + s = scaledSphericalHarmonicModel r scs + scs = [[(0,0)], + [(10.3, 0), (18.1, -26.6)], + [(-8.7, 0), (-3.3, -27.4), (2.1, -14.1)], + [(3.4, 0), (-5.5, 8.2), (-0.7, -0.4), (-10.1, 1.8)], + [(-0.7, 0), (0.2, -1.3), (-9.1, 5.3), (4.1, 2.9), (-4.3, -5.2)], + [(-0.2, 0), (0.5, 0.6), (-1.3, 1.7), (-0.1, -1.2), (1.4, 3.4), (3.9, 0)], + [(-0.3, 0), (-0.1, 0), (-0.7, -2.1), (2.1, -0.7), (-1.2, 0.2), (0.3, 0.9), (1.6, 1)], + [(0.3, 0), (-0.2, 0.8), (-0.5, 0.4), (1.3, -0.2), (0.1, -0.3), (-0.6, -0.6), (-0.8, 0.1), (0.2, -0.2)], + [(0.2, 0), (0, -0.3), (-0.6, 0.3), (0.5, 0.1), (-0.2, 0.5), (0.4, -0.2), (0.1, -0.3), (-0.4, 0.3), (0.3, 0)] + ] + r = 6371.2 + +-- | The International Geomagnetic Reference Field model, 13th edition. +-- Model epoch is January 1st, 2020. +igrf13 :: (Fractional a) => MagneticModel a +igrf13 = MagneticModel + { + fieldAtEpoch = f, + secularVariation = s + } + where + f = scaledSphericalHarmonicModel r fcs + fcs = [[(0.0,0.0)], + [(-29404.8,0.0),(-1450.9,4652.5)], + [(-2499.6,0.0),(2982.0,-2991.6),(1677.0,-734.6)], + [(1363.2,0.0),(-2381.2,-82.1),(1236.2,241.9),(525.7,-543.4)], + [(903.0,0.0),(809.5,281.9),(86.3,-158.4),(-309.4,199.7),(48.0,-349.7)], + [(-234.3,0.0),(363.2,47.7),(187.8,208.3),(-140.7,-121.2),(-151.2,32.3),(13.5,98.9)], + [(66.0,0.0),(65.5,-19.1),(72.9,25.1),(-121.5,52.8),(-36.2,-64.5),(13.5,8.9),(-64.7,68.1)], + [(80.6,0.0),(-76.7,-51.5),(-8.2,-16.9),(56.5,2.2),(15.8,23.5),(6.4,-2.2),(-7.2,-27.2),(9.8,-1.8)], + [(23.7,0.0),(9.7,8.4),(-17.6,-15.3),(-0.5,12.8),(-21.1,-11.7),(15.3,14.9),(13.7,3.6),(-16.5,-6.9),(-0.3,2.8)], + [(5.0,0.0),(8.4,-23.4),(2.9,11.0),(-1.5,9.8),(-1.1,-5.1),(-13.2,-6.3),(1.1,7.8),(8.8,0.4),(-9.3,-1.4),(-11.9,9.6)], + [(-1.9,0.0),(-6.2,3.4),(-0.1,-0.2),(1.7,3.6),(-0.9,4.8),(0.7,-8.6),(-0.9,-0.1),(1.9,-4.3),(1.4,-3.4),(-2.4,-0.1),(-3.8,-8.8)], + [(3.0,0.0),(-1.4,0.0),(-2.5,2.5),(2.3,-0.6),(-0.9,-0.4),(0.3,0.6),(-0.7,-0.2),(-0.1,-1.7),(1.4,-1.6),(-0.6,-3.0),(0.2,-2.0),(3.1,-2.6)], + [(-2.0,0.0),(-0.1,-1.2),(0.5,0.5),(1.3,1.4),(-1.2,-1.8),(0.7,0.1),(0.3,0.8),(0.5,-0.2),(-0.3,0.6),(-0.5,0.2),(0.1,-0.9),(-1.1,0.0),(-0.3,0.5)], + [(0.1,0.0),(-0.9,-0.9),(0.5,0.6),(0.7,1.4),(-0.3,-0.4),(0.8,-1.3),(0.0,-0.1),(0.8,0.3),(0.0,-0.1),(0.4,0.5),(0.1,0.5),(0.5,-0.4),(-0.5,-0.4),(-0.4,-0.6)] + ] + s = scaledSphericalHarmonicModel r scs + scs = [[(0.0,0.0)], + [(5.7,0.0),(7.4,-25.9)], + [(-11.0,0.0),(-7.0,-30.2),(-2.1,-22.4)], + [(2.2,0.0),(-5.9,6.0),(3.1,-1.1),(-12.0,0.5)], + [(-1.2,0.0),(-1.6,-0.1),(-5.9,6.5),(5.2,3.6),(-5.1,-5.0)], + [(-0.3,0.0),(0.5,0.0),(-0.6,2.5),(0.2,-0.6),(1.3,3.0),(0.9,0.3)], + [(-0.5,0.0),(-0.3,0.0),(0.4,-1.6),(1.3,-1.3),(-1.4,0.8),(0.0,0.0),(0.9,1.0)], + [(-0.1,0.0),(-0.2,0.6),(0.0,0.6),(0.7,-0.8),(0.1,-0.2),(-0.5,-1.1),(-0.8,0.1),(0.8,0.3)], + [(0.0,0.0),(0.1,-0.2),(-0.1,0.6),(0.4,-0.2),(-0.1,0.5),(0.4,-0.3),(0.3,-0.4),(-0.1,0.5),(0.4,0.0)] + ] r = 6371.2
+ src/IGRF/Parser.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE OverloadedStrings #-} + +-- | Parse <https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt> +module IGRF.Parser (parseModels) where + +import Control.Arrow +import Math.SphericalHarmonics +import Data.Text (Text) +import qualified Data.List as L +import qualified Data.Text as T + +-- | Parse <https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt> +-- and return a list of models. +parseModels :: Text -> [(Text, Text, SphericalHarmonicModel Double)] +parseModels file = map (parseModel . selectColumn) [3..length (head nonComments) - 1] + where + nonComments = fmap T.words $ filter (not . T.isPrefixOf "#") $ T.lines file + selectColumn i = fmap (\xs -> (xs !! 0, xs !! 1, xs !! 2, xs !! i)) nonComments + +parseModel :: [(Text, Text, Text, Text)] -> (Text, Text, SphericalHarmonicModel Double) +parseModel ((_, _, _, header1) : (_, _, _, header2) : values) = (header1, header2, sphericalHarmonicModel model) + where + n :: Int + n = maximum $ map (\(_, x, _, _) -> read (T.unpack x)) values + + zeroModel :: [[(Double, Double)]] + zeroModel = map (\i -> replicate (i + 1) (0, 0)) [0..n] + + model :: [[(Double, Double)]] + model = L.foldl' (flip go) zeroModel values + + go :: (Text, Text, Text, Text) -> [[(Double, Double)]] -> [[(Double, Double)]] + go (gh, i, j, x) = modify + ((if gh == "g" then first else second) $ const $ read $ T.unpack x) + (read $ T.unpack i) + (read $ T.unpack j) + +modify :: (a -> a) -> Int -> Int -> [[a]] -> [[a]] +modify f i j xss = xss' + where + xs = xss !! i + x = xs !! j + x' = f x + xs' = take j xs <> [x'] <> drop (j + 1) xs + xss' = take i xss <> [xs'] <> drop (i + 1) xss
src/Math/SphericalHarmonics.hs view
@@ -1,76 +1,99 @@ {-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE TypeFamilies #-} -- | Provides spherical harmonic models of scalar-valued functions. module Math.SphericalHarmonics ( - SphericalHarmonicModel(..) -, combine -, scale + SphericalHarmonicModel +, sphericalHarmonicModel +, scaledSphericalHarmonicModel , evaluateModel +, evaluateModelCartesian , evaluateModelGradient +, evaluateModelGradientCartesian , evaluateModelGradientInLocalTangentPlane ) where +import Data.Complex +import Data.VectorSpace hiding (magnitude) 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. - } +data SphericalHarmonicModel a = SphericalHarmonicModel [[(a, a)]] 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) - } +-- | Creates a spherical harmonic model. +-- Result in an error if the length of the list is not a triangular number. +sphericalHarmonicModel :: (Fractional a) => [[(a, a)]] -- ^ A list of g and h coefficients for the model + -> SphericalHarmonicModel a -- ^ The spherical harmonic model +sphericalHarmonicModel cs | valid = SphericalHarmonicModel cs + | otherwise = error "The number of coefficients is not a triangular number." 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) + valid = and $ zipWith (==) (fmap length cs) [1..length cs] --- | Linearly scales a spherical harmonic model. -scale :: (Num a) => a -> SphericalHarmonicModel a -> SphericalHarmonicModel a -scale x m = m { coefficients = fmap scalePair (coefficients m) } +-- | Creates a spherical harmonic model, scaling coefficients for the supplied reference radius. +-- Result in an error if the length of the list is not a triangular number. +scaledSphericalHarmonicModel :: (Fractional a) => a -- ^ The reference radius + -> [[(a, a)]] -- ^ A list of g and h coefficients for the model + -> SphericalHarmonicModel a -- ^ The spherical harmonic model +scaledSphericalHarmonicModel r cs = sphericalHarmonicModel cs' where - scalePair (g, h) = (x * g, x * h) + cs' = normalizeReferenceRadius r cs +instance(Fractional a, Eq a) => AdditiveGroup (SphericalHarmonicModel a) where + zeroV = SphericalHarmonicModel [[(0,0)]] + negateV = fmap negate + (SphericalHarmonicModel m1) ^+^ (SphericalHarmonicModel m2) = SphericalHarmonicModel (combineCoefficients m1 m2) + where + combineCoefficients [] cs = cs + combineCoefficients cs [] = cs + combineCoefficients (c1:cs1) (c2:cs2) = zipWith addPairs c1 c2 : combineCoefficients cs1 cs2 + addPairs (g1, h1) (g2, h2) = (g1 + g2, h1 + h2) + +instance (Fractional a, Eq a) => VectorSpace (SphericalHarmonicModel a) where + type Scalar (SphericalHarmonicModel a) = a + x *^ m = fmap (* x) m + +normalizeReferenceRadius :: (Fractional a) => a -> [[(a, a)]] -> [[(a, a)]] +normalizeReferenceRadius r = zipWith (fmap . mapWholePair . transform) [0 :: Int ..] + where + transform n = (* (r ^ (2 + n))) + -- | Computes the scalar value of the spherical harmonic model at a specified spherical position. -evaluateModel :: (Floating a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model +evaluateModel :: (RealFloat 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 +evaluateModel m r colat lon = evaluateModel' m r (cos colat) (cis lon) + +-- | Computes the scalar value of the spherical harmonic model at a specified Cartesian position. +evaluateModelCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model + -> a -- ^ X position + -> a -- ^ Y position + -> a -- ^ Z position + -> a -- ^ Model value +evaluateModelCartesian m x y z = evaluateModel' m r cosColat cisLon 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 + r = sqrt $ (x*x) + (y*y) + (z*z) + cosColat = z / r + cisLon = normalize $ mkPolar x y +evaluateModel' :: (RealFloat a, Ord a) => SphericalHarmonicModel a + -> a -- r + -> a -- cosColat + -> Complex a -- cisLon + -> a +evaluateModel' (SphericalHarmonicModel cs) r cosColat cisLon = sum $ zipWith (*) (iterate (/ r) (recip r)) (zipWith evaluateDegree [0..] cs) + where + sines = 1 : iterate (* cisLon) cisLon + evaluateDegree n cs' = sum $ zipWith3 evaluateOrder (fmap (schmidtSemiNormalizedAssociatedLegendreFunction n) [0..n]) cs' sines + evaluateOrder p (g, h) cisMLon = ((g * realPart cisMLon) + (h * imagPart cisMLon)) * (p (cosColat)) + -- | 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 +evaluateModelGradient :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model -> a -- ^ Spherical radius -> a -- ^ Spherical colatitude (radian) -> a -- ^ Spherical longitude (radian) @@ -78,11 +101,22 @@ 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 +-- The result is expressed in right-handed coordinates centered at the origin of the sphere, with the positive Z-axis piercing the +-- north pole and the positive x-axis piercing the reference meridian. +evaluateModelGradientCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model + -> a -- ^ X position + -> a -- ^ Y position + -> a -- ^ Z position + -> (a, a, a) -- X, Y, and Z components of gradient +evaluateModelGradientCartesian model x y z = makeTuple . fmap negate $ modelGrad [x, y, z] + where + modelGrad = grad (\[x', y', z'] -> evaluateModelCartesian (fmap auto model) 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 sphere at the specified location. +evaluateModelGradientInLocalTangentPlane :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model -> a -- ^ Spherical radius -> a -- ^ Spherical colatitude (radian) -> a -- ^ Spherical longitude (radian) @@ -94,8 +128,14 @@ n = -colat' / r -- negated because the colatitude increase southward u = r' -computeIndex :: Int -> Int -> Int -computeIndex n m = triangle n + m +normalize :: (RealFloat a) => Complex a -> Complex a +normalize r@(x :+ y) | isInfinite m' = 0 + | otherwise = (x * m') :+ (y * m') + where + m' = recip . magnitude $ r -triangle :: Int -> Int -triangle n = (n * (n + 1)) `div` 2+mapWholePair :: (a -> b) -> (a, a) -> (b, b) +mapWholePair f (a, b) = (f a, f b) + +makeTuple :: [a] -> (a, a, a) +makeTuple [x, y, z] = (x, y, z)
src/Math/SphericalHarmonics/AssociatedLegendre.hs view
@@ -6,8 +6,9 @@ ) where -import Math.Polynomial hiding (x) -import Math.Polynomial.Legendre +import Data.Poly (VPoly, eval, deriv) +import Data.Poly.Orthogonal (legendre) +import Data.Euclidean (Field, WrappedFractional(..)) -- definition from http://www.mathworks.com/help/matlab/ref/legendre.html#f89-998354 -- | Computes the associated Legendre function of degree 'n' and order 'm'. @@ -18,10 +19,11 @@ -> 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 + f x = nonPolyTerm x * unwrapFractional (eval p' (WrapFractional x)) + nonPolyTerm x = (1 - x * x) ** (fromIntegral m / 2) + p' = iterate deriv p !! m + p :: (Eq t, Field t) => VPoly t + 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'. @@ -32,8 +34,8 @@ 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) - + 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)