module Data.Geo.Jord.Earth
    ( Earth(..)
    , Ellipsoid(..)
    , eccentricity
    , meanRadius
    , polarRadius
    , spherical
    
    , wgs84
    , grs80
    , wgs72
    
    , s84
    , s80
    , s72
    , r84
    , r80
    , r72
    ) where
import Data.Geo.Jord.Length
data Earth
    = Ellipsoidal Ellipsoid
    | Spherical Length
    deriving (Eq, Show)
data Ellipsoid = Ellipsoid
    { equatorialRadius :: Length 
    , inverseFlattening :: Double 
    } deriving (Eq, Show)
eccentricity :: Earth -> Double
eccentricity (Ellipsoidal e) = sqrt (1.0 - (b * b) / (a * a))
  where
    a = semiMajorAxis e
    b = semiMinorAxis a (inverseFlattening e)
eccentricity (Spherical _) = 0
meanRadius :: Earth -> Length
meanRadius (Ellipsoidal e) = metres ((2.0 * a + b) / 3.0)
  where
    a = semiMajorAxis e
    b = semiMinorAxis a (inverseFlattening e)
meanRadius (Spherical r) = r
polarRadius :: Earth -> Length
polarRadius (Ellipsoidal e) = metres (semiMinorAxis a f)
  where
    a = semiMajorAxis e
    f = inverseFlattening e
polarRadius (Spherical r) = r
spherical :: Earth -> Earth
spherical e = Spherical (meanRadius e)
wgs84 :: Earth
wgs84 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257223563))
grs80 :: Earth
grs80 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257222101))
wgs72 :: Earth
wgs72 = Ellipsoidal (Ellipsoid (metres 6378135.0) (1.0 / 298.26))
s84 :: Earth
s84 = spherical wgs84
s80 :: Earth
s80 = spherical grs80
s72 :: Earth
s72 = spherical wgs72
r84 :: Length
r84 = meanRadius s84
r80 :: Length
r80 = meanRadius s80
r72 :: Length
r72 = meanRadius s72
semiMajorAxis :: Ellipsoid -> Double
semiMajorAxis = toMetres . equatorialRadius
semiMinorAxis :: Double -> Double -> Double
semiMinorAxis a f = a * (1.0 - f)