{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RoleAnnotations #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}

{- | An Ellipsoid is a reasonable best fit for the surface of the
Earth over some defined area. WGS84 is the standard used for the whole
of the Earth. Other Ellipsoids are considered a best fit for some
specific area.
-}

module Geodetics.Ellipsoids (
   -- * Useful constants
   degree,
   arcminute,
   arcsecond,
   kilometer,
   _2, _3, _4, _5, _6, _7,
   -- * Helmert transform between geodetic reference systems
   Helmert (..),
   inverseHelmert,
   ECEF,
   applyHelmert,
   -- * Ellipsoid models of the Geoid
   Ellipsoid (..),
   WGS84 (..),
   LocalEllipsoid (..),
   flattening,
   minorRadius,
   eccentricity2,
   eccentricity'2,
   -- * Auxiliary latitudes and related Values
   normal,
   latitudeRadius,
   meridianRadius,
   primeVerticalRadius,
   isometricLatitude,
   -- * Tiny linear algebra library for 3D vectors
   Vec3,
   Matrix3,
   add3,
   scale3,
   negate3,
   transform3,
   invert3,
   trans3,
   dot3,
   cross3
) where


-- | All angles in this library are in radians. This is one degree in radians.
degree :: Double
degree :: Double
degree = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
180

-- | One arc-minute in radians.
arcminute :: Double
arcminute :: Double
arcminute = Double
degree Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
60

-- | One arc-second in radians.
arcsecond :: Double
arcsecond :: Double
arcsecond = Double
arcminute Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
60


-- | All distances in this library are in meters. This is one kilometer in meters.
kilometer :: Double
kilometer :: Double
kilometer = Double
1000

-- | Lots of geodetic calculations involve integer powers. Writing e.g. @x ^ 2@ causes
-- GHC to complain that the @2@ has ambiguous type. @x ** 2@ doesn't complain
-- but is much slower. So for convenience, here are small integers with type @Int@.
_2, _3, _4, _5, _6, _7 :: Int
_2 :: Int
_2 = Int
2
_3 :: Int
_3 = Int
3
_4 :: Int
_4 = Int
4
_5 :: Int
_5 = Int
5
_6 :: Int
_6 = Int
6
_7 :: Int
_7 = Int
7

-- | 3d vector as @(X,Y,Z)@.
type Vec3 a = (a,a,a)

-- | 3x3 transform matrix for Vec3.
type Matrix3 a = Vec3 (Vec3 a)


-- | Multiply a vector by a scalar.
scale3 :: (Num a) =>  Vec3 a -> a -> Vec3 a
scale3 :: forall a. Num a => Vec3 a -> a -> Vec3 a
scale3 (a
x,a
y,a
z) a
s = (a
xa -> a -> a
forall a. Num a => a -> a -> a
*a
s, a
ya -> a -> a
forall a. Num a => a -> a -> a
*a
s, a
za -> a -> a
forall a. Num a => a -> a -> a
*a
s)


-- | Negation of a vector.
negate3 :: (Num a) => Vec3 a -> Vec3 a
negate3 :: forall a. Num a => Vec3 a -> Vec3 a
negate3 (a
x,a
y,a
z) = (a -> a
forall a. Num a => a -> a
negate a
x, a -> a
forall a. Num a => a -> a
negate a
y, a -> a
forall a. Num a => a -> a
negate a
z)

-- | Add two vectors
add3 :: (Num a) => Vec3 a -> Vec3 a -> Vec3 a
add3 :: forall a. Num a => Vec3 a -> Vec3 a -> Vec3 a
add3 (a
x1,a
y1,a
z1) (a
x2,a
y2,a
z2) = (a
x1a -> a -> a
forall a. Num a => a -> a -> a
+a
x2, a
y1a -> a -> a
forall a. Num a => a -> a -> a
+a
y2, a
z1a -> a -> a
forall a. Num a => a -> a -> a
+a
z2)


-- | Multiply a matrix by a vector in the Dimensional type system.
transform3 :: (Num a) =>
   Matrix3 a -> Vec3 a -> Vec3 a
transform3 :: forall a. Num a => Matrix3 a -> Vec3 a -> Vec3 a
transform3 (Vec3 a
tx,Vec3 a
ty,Vec3 a
tz) Vec3 a
v = (Vec3 a -> Vec3 a -> a
forall {a}. Num a => (a, a, a) -> (a, a, a) -> a
t Vec3 a
tx Vec3 a
v, Vec3 a -> Vec3 a -> a
forall {a}. Num a => (a, a, a) -> (a, a, a) -> a
t Vec3 a
ty Vec3 a
v, Vec3 a -> Vec3 a -> a
forall {a}. Num a => (a, a, a) -> (a, a, a) -> a
t Vec3 a
tz Vec3 a
v)
   where
      t :: (a, a, a) -> (a, a, a) -> a
t (a
x1,a
y1,a
z1) (a
x2,a
y2,a
z2) = a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2


-- | Inverse of a 3x3 matrix.
invert3 :: (Fractional a) => Matrix3 a -> Matrix3 a
invert3 :: forall a. Fractional a => Matrix3 a -> Matrix3 a
invert3 ((a
x1,a
y1,a
z1),
         (a
x2,a
y2,a
z2),
         (a
x3,a
y3,a
z3)) =
      ((a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
y2 a
z2 a
y3 a
z3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
z1 a
y1 a
z3 a
y3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
y1 a
z1 a
y2 a
z2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det),
       (a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
z2 a
x2 a
z3 a
x3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
x1 a
z1 a
x3 a
z3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
z1 a
x1 a
z2 a
x2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det),
       (a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
x2 a
y2 a
x3 a
y3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
y1 a
x1 a
y3 a
x3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det, a -> a -> a -> a -> a
forall {a}. Num a => a -> a -> a -> a -> a
det2 a
x1 a
y1 a
x2 a
y2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
det))
   where
      det :: a
det = (a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2a -> a -> a
forall a. Num a => a -> a -> a
*a
z3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2a -> a -> a
forall a. Num a => a -> a -> a
*a
x3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2a -> a -> a
forall a. Num a => a -> a -> a
*a
y3) a -> a -> a
forall a. Num a => a -> a -> a
- (a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2a -> a -> a
forall a. Num a => a -> a -> a
*a
x3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2a -> a -> a
forall a. Num a => a -> a -> a
*a
z3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2a -> a -> a
forall a. Num a => a -> a -> a
*a
y3)
      det2 :: a -> a -> a -> a -> a
det2 a
a a
b a
c a
d = a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
ba -> a -> a
forall a. Num a => a -> a -> a
*a
c

-- | Transpose of a 3x3 matrix.
trans3 :: Matrix3 a -> Matrix3 a
trans3 :: forall a. Matrix3 a -> Matrix3 a
trans3 ((a
x1,a
y1,a
z1),(a
x2,a
y2,a
z2),(a
x3,a
y3,a
z3)) = ((a
x1,a
x2,a
x3),(a
y1,a
y2,a
y3),(a
z1,a
z2,a
z3))


-- | Dot product of two vectors
dot3 :: (Num a) => Vec3 a -> Vec3 a -> a
dot3 :: forall {a}. Num a => (a, a, a) -> (a, a, a) -> a
dot3 (a
x1,a
y1,a
z1) (a
x2,a
y2,a
z2) = a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2

-- | Cross product of two vectors
cross3 :: (Num a) => Vec3 a -> Vec3 a -> Vec3 a
cross3 :: forall a. Num a => Vec3 a -> Vec3 a -> Vec3 a
cross3 (a
x1,a
y1,a
z1) (a
x2,a
y2,a
z2) = (a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2 a -> a -> a
forall a. Num a => a -> a -> a
- a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2, a
z1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2 a -> a -> a
forall a. Num a => a -> a -> a
- a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
z2, a
x1a -> a -> a
forall a. Num a => a -> a -> a
*a
y2 a -> a -> a
forall a. Num a => a -> a -> a
- a
y1a -> a -> a
forall a. Num a => a -> a -> a
*a
x2)


-- | The 7 parameter Helmert transformation. The monoid instance allows composition but
-- is only accurate for the small values used in practical ellipsoids.
data Helmert = Helmert {
   Helmert -> Double
cX, Helmert -> Double
cY, Helmert -> Double
cZ :: Double,  -- ^ Offset in meters
   Helmert -> Double
helmertScale :: Double,  -- ^ Parts per million
   Helmert -> Double
rX, Helmert -> Double
rY, Helmert -> Double
rZ :: Double  -- ^ Rotation around each axis in radians.
} deriving (Helmert -> Helmert -> Bool
(Helmert -> Helmert -> Bool)
-> (Helmert -> Helmert -> Bool) -> Eq Helmert
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Helmert -> Helmert -> Bool
== :: Helmert -> Helmert -> Bool
$c/= :: Helmert -> Helmert -> Bool
/= :: Helmert -> Helmert -> Bool
Eq, Int -> Helmert -> ShowS
[Helmert] -> ShowS
Helmert -> String
(Int -> Helmert -> ShowS)
-> (Helmert -> String) -> ([Helmert] -> ShowS) -> Show Helmert
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Helmert -> ShowS
showsPrec :: Int -> Helmert -> ShowS
$cshow :: Helmert -> String
show :: Helmert -> String
$cshowList :: [Helmert] -> ShowS
showList :: [Helmert] -> ShowS
Show)

instance Semigroup Helmert where
    Helmert
h1 <> :: Helmert -> Helmert -> Helmert
<> Helmert
h2 = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Helmert
Helmert (Helmert -> Double
cX Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
cX Helmert
h2) (Helmert -> Double
cY Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
cY Helmert
h2) (Helmert -> Double
cZ Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
cZ Helmert
h2)
                       (Helmert -> Double
helmertScale Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
helmertScale Helmert
h2)
                       (Helmert -> Double
rX Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
rX Helmert
h2) (Helmert -> Double
rY Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
rY Helmert
h2) (Helmert -> Double
rZ Helmert
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
rZ Helmert
h2)

instance Monoid Helmert where
   mempty :: Helmert
mempty = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Helmert
Helmert Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0
   mappend :: Helmert -> Helmert -> Helmert
mappend = Helmert -> Helmert -> Helmert
forall a. Semigroup a => a -> a -> a
(<>)

-- | The inverse of a Helmert transformation.
inverseHelmert :: Helmert -> Helmert
inverseHelmert :: Helmert -> Helmert
inverseHelmert Helmert
h = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Helmert
Helmert (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
cX Helmert
h) (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
cY Helmert
h) (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
cZ Helmert
h)
                           (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
helmertScale Helmert
h)
                           (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
rX Helmert
h) (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
rY Helmert
h) (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Helmert -> Double
rZ Helmert
h)


-- | Earth-centred, Earth-fixed coordinates as a vector. The origin and axes are
-- not defined: use with caution.
type ECEF = Vec3 Double

-- | Apply a Helmert transformation to earth-centered coordinates.
applyHelmert:: Helmert -> ECEF -> ECEF
applyHelmert :: Helmert -> ECEF -> ECEF
applyHelmert Helmert
h (Double
x,Double
y,Double
z) = (
      Helmert -> Double
cX Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (                Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Helmert -> Double
rZ Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
rY Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z),
      Helmert -> Double
cY Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (        Helmert -> Double
rZ Helmert
h  Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+        Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Helmert -> Double
rX Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z),
      Helmert -> Double
cZ Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
negate (Helmert -> Double
rY Helmert
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
rX Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+        Double
z))
   where
      s :: Double
s = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Helmert -> Double
helmertScale Helmert
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
1e-6


-- | An Ellipsoid is defined by the major radius and the inverse flattening (which define its shape),
-- and its Helmert transform relative to WGS84 (which defines its position and orientation).
--
-- The inclusion of the Helmert parameters relative to WGS84 actually make this a Terrestrial
-- Reference Frame (TRF), but the term "Ellipsoid" will be used in this library for readability.
--
-- Minimum definition: @majorRadius@, @flatR@ & @helmert@.
--
-- Laws:
--
-- > helmertToWGS84 = applyHelmert . helmert
-- > helmertFromWGS84 e . helmertToWGS84 e = id
class (Show a, Eq a) => Ellipsoid a where
   majorRadius :: a -> Double
   flatR :: a -> Double
      -- ^ Inverse of the flattening.
   helmert :: a -> Helmert  -- ^ The Helmert parameters relative to WGS84,
   helmertToWGS84 :: a -> ECEF -> ECEF
      -- ^ The Helmert transform that will convert a position wrt
      -- this ellipsoid into a position wrt WGS84.
   helmertToWGS84 a
e = Helmert -> ECEF -> ECEF
applyHelmert (a -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert a
e)
   helmertFromWGS84 :: a -> ECEF -> ECEF
      -- ^ And its inverse.
   helmertFromWGS84 a
e = Helmert -> ECEF -> ECEF
applyHelmert (Helmert -> Helmert
inverseHelmert (Helmert -> Helmert) -> Helmert -> Helmert
forall a b. (a -> b) -> a -> b
$ a -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert a
e)


-- | The WGS84 geoid, major radius 6378137.0 meters, flattening = 1 / 298.257223563
-- as defined in \"Technical Manual DMA TM 8358.1 - Datums, Ellipsoids, Grids, and
-- Grid Reference Systems\" at the National Geospatial-Intelligence Agency (NGA).
--
-- The WGS84 has a special place in this library as the standard Ellipsoid against
-- which all others are defined.
data WGS84 = WGS84
  deriving (WGS84 -> WGS84 -> Bool
(WGS84 -> WGS84 -> Bool) -> (WGS84 -> WGS84 -> Bool) -> Eq WGS84
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: WGS84 -> WGS84 -> Bool
== :: WGS84 -> WGS84 -> Bool
$c/= :: WGS84 -> WGS84 -> Bool
/= :: WGS84 -> WGS84 -> Bool
Eq, Int -> WGS84 -> ShowS
[WGS84] -> ShowS
WGS84 -> String
(Int -> WGS84 -> ShowS)
-> (WGS84 -> String) -> ([WGS84] -> ShowS) -> Show WGS84
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> WGS84 -> ShowS
showsPrec :: Int -> WGS84 -> ShowS
$cshow :: WGS84 -> String
show :: WGS84 -> String
$cshowList :: [WGS84] -> ShowS
showList :: [WGS84] -> ShowS
Show)

instance Ellipsoid WGS84 where
   majorRadius :: WGS84 -> Double
majorRadius WGS84
_ = Double
6378137.0
   flatR :: WGS84 -> Double
flatR WGS84
_ = Double
298.257223563
   helmert :: WGS84 -> Helmert
helmert WGS84
_ = Helmert
forall a. Monoid a => a
mempty
   helmertToWGS84 :: WGS84 -> ECEF -> ECEF
helmertToWGS84 WGS84
_ = ECEF -> ECEF
forall a. a -> a
id
   helmertFromWGS84 :: WGS84 -> ECEF -> ECEF
helmertFromWGS84 WGS84
_ = ECEF -> ECEF
forall a. a -> a
id


-- | Ellipsoids other than WGS84, used within a defined geographical area where
-- they are a better fit to the local geoid. Can also be used for historical ellipsoids.
--
-- The @Show@ instance just returns the name.
-- Creating two different local ellipsoids with the same name is a Bad Thing.
data LocalEllipsoid = LocalEllipsoid {
   LocalEllipsoid -> String
nameLocal :: String,
   LocalEllipsoid -> Double
majorRadiusLocal :: Double,
   LocalEllipsoid -> Double
flatRLocal :: Double,
   LocalEllipsoid -> Helmert
helmertLocal :: Helmert
} deriving (LocalEllipsoid -> LocalEllipsoid -> Bool
(LocalEllipsoid -> LocalEllipsoid -> Bool)
-> (LocalEllipsoid -> LocalEllipsoid -> Bool) -> Eq LocalEllipsoid
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LocalEllipsoid -> LocalEllipsoid -> Bool
== :: LocalEllipsoid -> LocalEllipsoid -> Bool
$c/= :: LocalEllipsoid -> LocalEllipsoid -> Bool
/= :: LocalEllipsoid -> LocalEllipsoid -> Bool
Eq)

instance Show LocalEllipsoid where
    show :: LocalEllipsoid -> String
show = LocalEllipsoid -> String
nameLocal

instance Ellipsoid LocalEllipsoid where
   majorRadius :: LocalEllipsoid -> Double
majorRadius = LocalEllipsoid -> Double
majorRadiusLocal
   flatR :: LocalEllipsoid -> Double
flatR = LocalEllipsoid -> Double
flatRLocal
   helmert :: LocalEllipsoid -> Helmert
helmert = LocalEllipsoid -> Helmert
helmertLocal


-- | Flattening (f) of an ellipsoid.
flattening :: (Ellipsoid e) => e -> Double
flattening :: forall e. Ellipsoid e => e -> Double
flattening e
e = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ e -> Double
forall e. Ellipsoid e => e -> Double
flatR e
e

-- | The minor radius of an ellipsoid in meters.
minorRadius :: (Ellipsoid e) => e -> Double
minorRadius :: forall e. Ellipsoid e => e -> Double
minorRadius e
e = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
flattening e
e)


-- | The eccentricity squared of an ellipsoid.
eccentricity2 :: (Ellipsoid e) => e -> Double
eccentricity2 :: forall e. Ellipsoid e => e -> Double
eccentricity2 e
e = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Int
_2 where f :: Double
f = e -> Double
forall e. Ellipsoid e => e -> Double
flattening e
e

-- | The second eccentricity squared of an ellipsoid.
eccentricity'2 :: (Ellipsoid e) => e -> Double
eccentricity'2 :: forall e. Ellipsoid e => e -> Double
eccentricity'2 e
e = (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Int
_2) where f :: Double
f = e -> Double
forall e. Ellipsoid e => e -> Double
flattening e
e


-- | Distance in meters from the surface at the specified latitude to the
-- axis of the Earth straight down. Also known as the radius of
-- curvature in the prime vertical, and often denoted @N@.
normal :: (Ellipsoid e) => e -> Double -> Double
normal :: forall e. Ellipsoid e => e -> Double -> Double
normal e
e Double
lat = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)


-- | Radius of the circle of latitude: the distance from a point
-- at that latitude to the axis of the Earth, in meters.
latitudeRadius :: (Ellipsoid e) => e -> Double -> Double
latitudeRadius :: forall e. Ellipsoid e => e -> Double -> Double
latitudeRadius e
e Double
lat = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
normal e
e Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
lat


-- | Radius of curvature in the meridian at the specified latitude, in meters
-- Often denoted @M@.
meridianRadius :: (Ellipsoid e) => e -> Double -> Double
meridianRadius :: forall e. Ellipsoid e => e -> Double -> Double
meridianRadius e
e Double
lat =
   e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e)
   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3)


-- | Radius of curvature of the ellipsoid perpendicular to the meridian at the specified latitude, in meters.
primeVerticalRadius :: (Ellipsoid e) => e -> Double -> Double
primeVerticalRadius :: forall e. Ellipsoid e => e -> Double -> Double
primeVerticalRadius e
e Double
lat =
   e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)


-- | The isometric latitude. The isometric latitude is conventionally denoted by ψ
-- (not to be confused with the geocentric latitude): it is used in the development
-- of the ellipsoidal versions of the normal Mercator projection and the Transverse
-- Mercator projection. The name "isometric" arises from the fact that at any point
-- on the ellipsoid equal increments of ψ and longitude λ give rise to equal distance
-- displacements along the meridians and parallels respectively.
isometricLatitude :: (Ellipsoid e) => e -> Double -> Double
isometricLatitude :: forall e. Ellipsoid e => e -> Double -> Double
isometricLatitude e
ellipse Double
lat = Double -> Double
forall a. Floating a => a -> a
atanh Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
atanh (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat)
   where
      sinLat :: Double
sinLat = Double -> Double
forall a. Floating a => a -> a
sin Double
lat
      e :: Double
e = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
ellipse