{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TypeFamilies #-}
module Math.SphericalHarmonics
(
SphericalHarmonicModel
, sphericalHarmonicModel
, scaledSphericalHarmonicModel
, evaluateModel
, evaluateModelCartesian
, evaluateModelGradient
, evaluateModelGradientCartesian
, evaluateModelGradientInLocalTangentPlane
)
where
import Data.Complex
import Data.VectorSpace hiding (magnitude)
import Math.SphericalHarmonics.AssociatedLegendre
import Numeric.AD
data SphericalHarmonicModel a = SphericalHarmonicModel [[(a, a)]]
deriving (Functor)
sphericalHarmonicModel :: (Fractional a) => [[(a, a)]]
-> SphericalHarmonicModel a
sphericalHarmonicModel cs | valid = SphericalHarmonicModel cs
| otherwise = error "The number of coefficients is not a triangular number."
where
valid = and $ zipWith (==) (fmap length cs) [1..length cs]
scaledSphericalHarmonicModel :: (Fractional a) => a
-> [[(a, a)]]
-> SphericalHarmonicModel a
scaledSphericalHarmonicModel r cs = sphericalHarmonicModel cs'
where
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)))
evaluateModel :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> a
-> a
evaluateModel m r colat lon = evaluateModel' m r (cos colat) (cis lon)
evaluateModelCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> a
-> a
evaluateModelCartesian m x y z = evaluateModel' m r cosColat cisLon
where
r = sqrt $ (x*x) + (y*y) + (z*z)
cosColat = z / r
cisLon = normalize $ mkPolar x y
evaluateModel' :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> Complex a
-> 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))
evaluateModelGradient :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> a
-> (a, a, a)
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')
evaluateModelGradientCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> a
-> (a, a, a)
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')
evaluateModelGradientInLocalTangentPlane :: (RealFloat a, Ord a) => SphericalHarmonicModel a
-> a
-> a
-> a
-> (a, a, a)
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
u = r'
normalize :: (RealFloat a) => Complex a -> Complex a
normalize r@(x :+ y) | isInfinite m' = 0
| otherwise = (x * m') :+ (y * m')
where
m' = recip . magnitude $ r
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)