module Geodetics.Geodetic (
   -- * Geodetic Coordinates
   Geodetic (..),
   readGroundPosition,
   toLocal,
   toWGS84,
   antipode,
   geometricalDistance,
   geometricalDistanceSq,
   groundDistance,
   properAngle,
   showAngle,
   -- * Earth Centred Earth Fixed Coordinates
   ECEF,
   geoToEarth,
   earthToGeo,
   -- * Re-exported for convenience
   WGS84 (..)
) where


import Data.Char (chr)
import Data.Maybe
import Geodetics.Altitude
import Geodetics.Ellipsoids
import Geodetics.LatLongParser
import Text.ParserCombinators.ReadP

-- | Defines a three-D position on or around the Earth using latitude,
-- longitude and altitude with respect to a specified ellipsoid, with
-- positive directions being North and East.  The default "show"
-- instance gives position in degrees, minutes and seconds to 5 decimal
-- places, which is a
-- resolution of about 1m on the Earth's surface. Internally latitude
-- and longitude are stored as double precision radians. Convert to
-- degrees using e.g.  @latitude g /~ degree@.
--
-- The functions here deal with altitude by assuming that the local
-- height datum is always co-incident with the ellipsoid in use,
-- even though the \"mean sea level\" (the usual height datum) can be tens
-- of meters above or below the ellipsoid, and two ellipsoids can
-- differ by similar amounts. This is because the altitude is
-- usually known with reference to a local datum regardless of the
-- ellipsoid in use, so it is simpler to preserve the altitude across
-- all operations. However if
-- you are working with ECEF coordinates from some other source then
-- this may give you the wrong results, depending on the altitude
-- correction your source has used.
--
-- There is no "Eq" instance because comparing two arbitrary
-- co-ordinates on the Earth is a non-trivial exercise. Clearly if all
-- the parameters are equal on the same ellipsoid then they are indeed
-- in the same place. However if different ellipsoids are used then
-- two co-ordinates with different numbers can still refer to the same
-- physical location.  If you want to find out if two co-ordinates are
-- the same to within a given tolerance then use "geometricDistance"
-- (or its squared variant to avoid an extra @sqrt@ operation).
data Geodetic e = Geodetic {
   forall e. Geodetic e -> Double
latitude, forall e. Geodetic e -> Double
longitude :: Double,  -- ^ In radians.
   forall e. Geodetic e -> Double
geoAlt :: Double,  -- ^ In meters.
   forall e. Geodetic e -> e
ellipsoid :: e
}

instance (Ellipsoid e) => Show (Geodetic e) where
   show :: Geodetic e -> String
show Geodetic e
g = [String] -> String
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [
      Double -> String
showAngle (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
g),  String
" ", String -> Double -> String
forall {a} {a}. (Ord a, Num a) => [a] -> a -> [a]
letter String
"SN" (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
g),  String
", ",
      Double -> String
showAngle (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
g), String
" ", String -> Double -> String
forall {a} {a}. (Ord a, Num a) => [a] -> a -> [a]
letter String
"WE" (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
g), String
", ",
      Double -> String
forall a. Show a => a -> String
show (Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
g), String
" ", e -> String
forall a. Show a => a -> String
show (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)]
      where letter :: [a] -> a -> [a]
letter [a]
s a
n = [[a]
s [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! (if a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then Int
0 else Int
1)]


-- | Read the latitude and longitude of a ground position and
-- return a Geodetic position on the specified ellipsoid.
--
-- The latitude and longitude may be in any of the following formats.
-- The comma between latitude and longitude is optional in all cases.
-- Latitude must always be first.
--
-- * Signed decimal degrees: 34.52327, -46.23234
--
-- * Decimal degrees NSEW: 34.52327N, 46.23234W
--
-- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92'
--
-- * Degrees, minutes and seconds (units optional): 34° 31' 23.52\" N, 46° 13' 56.43\" W
--
-- * DDDMMSS format with optional leading zeros: 343123.52N, 0461356.43W
readGroundPosition :: (Ellipsoid e) => e -> String -> Maybe (Geodetic e)
readGroundPosition :: forall e. Ellipsoid e => e -> String -> Maybe (Geodetic e)
readGroundPosition e
e String
str =
   case (((Double, Double), String) -> (Double, Double))
-> [((Double, Double), String)] -> [(Double, Double)]
forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double), String) -> (Double, Double)
forall a b. (a, b) -> a
fst ([((Double, Double), String)] -> [(Double, Double)])
-> [((Double, Double), String)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ (((Double, Double), String) -> Bool)
-> [((Double, Double), String)] -> [((Double, Double), String)]
forall a. (a -> Bool) -> [a] -> [a]
filter (String -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null (String -> Bool)
-> (((Double, Double), String) -> String)
-> ((Double, Double), String)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Double, Double), String) -> String
forall a b. (a, b) -> b
snd) ([((Double, Double), String)] -> [((Double, Double), String)])
-> [((Double, Double), String)] -> [((Double, Double), String)]
forall a b. (a -> b) -> a -> b
$ ReadP (Double, Double) -> ReadS (Double, Double)
forall a. ReadP a -> ReadS a
readP_to_S ReadP (Double, Double)
latLong String
str of
      [] -> Maybe (Geodetic e)
forall a. Maybe a
Nothing
      (Double
lat,Double
long) : [(Double, Double)]
_ -> Geodetic e -> Maybe (Geodetic e)
forall a. a -> Maybe a
Just (Geodetic e -> Maybe (Geodetic e))
-> Geodetic e -> Maybe (Geodetic e)
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Geodetic e
forall a. HasAltitude a => a -> a
groundPosition (Geodetic e -> Geodetic e) -> Geodetic e -> Geodetic e
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> e -> Geodetic e
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic (Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree) (Double
long Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree) Double
forall a. HasCallStack => a
undefined e
e


-- | Show an angle as degrees, minutes and seconds to two decimal places.
showAngle :: Double -> String
showAngle :: Double -> String
showAngle Double
a
   | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
a        = String
"NaN"  -- Not a Nangle
   | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
a   = String
sgn String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"Infinity"
   | Bool
otherwise      = [String] -> String
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [String
sgn, Integer -> String
forall a. Show a => a -> String
show Integer
d, [Int -> Char
chr Int
0xB0, Char
' '],
                              Integer -> String
forall a. Show a => a -> String
show Integer
m, String
"\8242 ",
                              Integer -> String
forall a. Show a => a -> String
show Integer
s, String
".", String
dstr, String
"\8243" ]
   where
      sgn :: String
sgn = if Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then String
"-" else String
""
      centisecs :: Integer
      centisecs :: Integer
centisecs = Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ (Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
arcsecond Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
100))
      (Integer
d, Integer
m1) = Integer
centisecs Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
360000
      (Integer
m, Integer
s1) = Integer
m1 Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
6000   -- hundredths of arcsec per arcmin
      (Integer
s, Integer
ds) = Integer
s1 Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
100
      dstr :: String
dstr = ShowS
forall a. [a] -> [a]
reverse ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$ Int -> ShowS
forall a. Int -> [a] -> [a]
take Int
2 ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$ ShowS
forall a. [a] -> [a]
reverse (Integer -> String
forall a. Show a => a -> String
show Integer
ds) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"00" -- Decimal fraction with zero padding.


instance (Ellipsoid e) => HasAltitude (Geodetic e) where
   altitude :: Geodetic e -> Double
altitude = Geodetic e -> Double
forall e. Geodetic e -> Double
geoAlt
   setAltitude :: Double -> Geodetic e -> Geodetic e
setAltitude Double
h Geodetic e
g = Geodetic e
g{geoAlt = h}



-- | The point on the Earth diametrically opposite the argument, with
-- the same altitude.
antipode :: (Ellipsoid e) => Geodetic e -> Geodetic e
antipode :: forall e. Ellipsoid e => Geodetic e -> Geodetic e
antipode Geodetic e
g = Double -> Double -> Double -> e -> Geodetic e
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
long (Geodetic e -> Double
forall e. Geodetic e -> Double
geoAlt Geodetic e
g) (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)
   where
      lat :: Double
lat = Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
g
      long' :: Double
long' = Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
g Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
180 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree
      long :: Double
long | Double
long' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0  = Double
long' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
360 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree
           | Bool
otherwise  = Double
long'



-- | Convert a geodetic coordinate into earth centered, relative to the
-- ellipsoid in use.
geoToEarth :: (Ellipsoid e) => Geodetic e -> ECEF
geoToEarth :: forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
geo = (
      (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslong,
      (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinlong,
      (Double
n 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. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinlat)
   where
      n :: Double
n = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
normal e
e (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      e :: e
e = Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
geo
      coslat :: Double
coslat = Double -> Double
forall a. Floating a => a -> a
cos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      coslong :: Double
coslong = Double -> Double
forall a. Floating a => a -> a
cos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo
      sinlat :: Double
sinlat = Double -> Double
forall a. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      sinlong :: Double
sinlong = Double -> Double
forall a. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo
      h :: Double
h = Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
geo


-- | Convert an earth centred coordinate into a geodetic coordinate on
-- the specified geoid.
--
-- Uses the closed form solution of H. Vermeille: Direct
-- transformation from geocentric coordinates to geodetic coordinates.
-- Journal of Geodesy Volume 76, Number 8 (2002), 451-454. Result is in the form
-- @(latitude, longitude, altitude)@.
earthToGeo :: (Ellipsoid e) => e -> ECEF -> (Double, Double, Double)
earthToGeo :: forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo e
e (Double
x,Double
y,Double
z) = (Double
phi, Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x, Double -> Double
forall a. Floating a => a -> a
sqrt (Double
l Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p2) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
norm)
   where
      -- Naming: numeric suffix inicates power. Hence x2 = x * x, x3 = x2 * x, etc.
      p2 :: Double
p2 = Double
x 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
* Double
y
      a :: Double
a = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e
      a2 :: Double
a2 = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a
      e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e
      e4 :: Double
e4 = Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e2
      zeta :: Double
zeta = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
e2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a2)
      rho :: Double
rho = (Double
p2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
zeta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e4) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6
      rho2 :: Double
rho2 = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho
      rho3 :: Double
rho3 = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho2
      s :: Double
s = Double
e4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zeta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a2)
      t :: Double
t = (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rho3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho3))) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3) -- Cube root
      u :: Double
u = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rho2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t
      v :: Double
v = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zeta)
      w :: Double
w = Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
zeta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
      kappa :: Double
kappa = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v)
      phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
atan (Double
kappa Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
p2)
      norm :: Double
norm = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
normal e
e Double
phi
      l :: Double
l = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
norm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi


-- | Convert a position from any geodetic to another one, assuming local altitude stays constant.
toLocal :: (Ellipsoid e1, Ellipsoid e2) => e2 -> Geodetic e1 -> Geodetic e2
toLocal :: forall e1 e2.
(Ellipsoid e1, Ellipsoid e2) =>
e2 -> Geodetic e1 -> Geodetic e2
toLocal e2
e2 Geodetic e1
g = Double -> Double -> Double -> e2 -> Geodetic e2
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
lon Double
alt e2
e2
   where
      alt :: Double
alt = Geodetic e1 -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e1
g
      (Double
lat, Double
lon, Double
_) = e2 -> ECEF -> ECEF
forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo e2
e2 (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Helmert -> ECEF -> ECEF
applyHelmert Helmert
h (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Geodetic e1 -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e1
g
      h :: Helmert
h = e1 -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert (Geodetic e1 -> e1
forall e. Geodetic e -> e
ellipsoid Geodetic e1
g) Helmert -> Helmert -> Helmert
forall a. Monoid a => a -> a -> a
`mappend` Helmert -> Helmert
inverseHelmert (e2 -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert e2
e2)

-- | Convert a position from any geodetic to WGS84, assuming local
-- altitude stays constant.
toWGS84 :: (Ellipsoid e) => Geodetic e -> Geodetic WGS84
toWGS84 :: forall e. Ellipsoid e => Geodetic e -> Geodetic WGS84
toWGS84 Geodetic e
g = Double -> Double -> Double -> WGS84 -> Geodetic WGS84
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
lon Double
alt WGS84
WGS84
   where
      alt :: Double
alt = Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
g
      (Double
lat, Double
lon, Double
_) = WGS84 -> ECEF -> ECEF
forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo WGS84
WGS84 (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Helmert -> ECEF -> ECEF
applyHelmert Helmert
h (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g
      h :: Helmert
h = e -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)


-- | The absolute distance in a straight line between two geodetic
-- points. They must be on the same ellipsoid.
-- Note that this is not the geodetic distance taken by following
-- the curvature of the earth.
geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
geometricalDistance :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistance Geodetic e
g1 Geodetic e
g2 = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Geodetic e -> Double
forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq Geodetic e
g1 Geodetic e
g2

-- | The square of the absolute distance.
geometricalDistanceSq :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq Geodetic e
g1 Geodetic e
g2 = (Double
x1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
y1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
z1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
z2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
   where
      (Double
x1,Double
y1,Double
z1) = Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g1
      (Double
x2,Double
y2,Double
z2) = Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g2


-- | The shortest ellipsoidal distance between two points on the
-- ground with reference to the same ellipsoid. Altitude is ignored.
--
-- The results are the distance between the points, the bearing of
-- the second point from the first, and (180 degrees - the bearing
-- of the first point from the second).
--
-- The algorithm can fail to converge where the arguments are near to
-- antipodal. In this case it returns @Nothing@.
--
-- Uses Vincenty's formula. \"Direct and inverse solutions of
-- geodesics on the ellipsoid with application of nested
-- equations\". T. Vincenty. Survey Review XXII 176, April
-- 1975. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>
groundDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e ->
                  Maybe (Double, Double, Double)
groundDistance :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Maybe ECEF
groundDistance Geodetic e
p1 Geodetic e
p2 = do
     ((Double, (Double, Double, Double, Double, Double))
_, (Double
lambda, (Double
cos2Alpha, Double
delta, Double
sinDelta, Double
cosDelta, Double
cos2DeltaM))) <-
       [((Double, (Double, Double, Double, Double, Double)),
  (Double, (Double, Double, Double, Double, Double)))]
-> Maybe
     ((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))
forall a. [a] -> Maybe a
listToMaybe ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> Maybe
      ((Double, (Double, Double, Double, Double, Double)),
       (Double, (Double, Double, Double, Double, Double))))
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> Maybe
     ((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))
forall a b. (a -> b) -> a -> b
$ (((Double, (Double, Double, Double, Double, Double)),
  (Double, (Double, Double, Double, Double, Double)))
 -> Bool)
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((Double, (Double, Double, Double, Double, Double)),
 (Double, (Double, Double, Double, Double, Double)))
-> Bool
forall {a} {b} {b}.
(Ord a, Fractional a) =>
((a, b), (a, b)) -> Bool
converging ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ Int
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a. Int -> [a] -> [a]
take Int
100 ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ [(Double, (Double, Double, Double, Double, Double))]
-> [(Double, (Double, Double, Double, Double, Double))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Double, (Double, Double, Double, Double, Double))]
lambdas ([(Double, (Double, Double, Double, Double, Double))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [(Double, (Double, Double, Double, Double, Double))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ Int
-> [(Double, (Double, Double, Double, Double, Double))]
-> [(Double, (Double, Double, Double, Double, Double))]
forall a. Int -> [a] -> [a]
drop Int
1 [(Double, (Double, Double, Double, Double, Double))]
lambdas
     let
       uSq :: Double
uSq = Double
cos2Alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
aDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
bDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
       bigA :: Double
bigA = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSqDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
16384 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
768) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
320 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
uSq))))
       bigB :: Double
bigB =     Double
uSqDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
1024  Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256  Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
128) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
74 Double -> Double -> Double
forall a. Num a => a -> a -> a
-  Double
47Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq))))
       deltaDelta :: Double
deltaDelta =
         Double
bigB Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
+
                             Double
bigBDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaMDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
                                       Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bigBDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDeltaDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3)
                                          Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3)))
       s :: Double
s = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bigA Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
deltaDelta)
       alpha1 :: Double
alpha1 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2(Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lambda) (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
lambda)
       alpha2 :: Double
alpha2 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2(Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lambda) (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2)
     ECEF -> Maybe ECEF
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double
s, Double
alpha1, Double
alpha2)
  where
    f :: Double
f = e -> Double
forall e. Ellipsoid e => e -> Double
flattening (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    a :: Double
a = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    b :: Double
b = e -> Double
forall e. Ellipsoid e => e -> Double
minorRadius (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    l :: Double
l = Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
p1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
p2
    u1 :: Double
u1 = Double -> Double
forall a. Floating a => a -> a
atan ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
f) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
p1))
    u2 :: Double
u2 = Double -> Double
forall a. Floating a => a -> a
atan ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
f) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
p2))
    sinU1 :: Double
sinU1 = Double -> Double
forall a. Floating a => a -> a
sin Double
u1
    cosU1 :: Double
cosU1 = Double -> Double
forall a. Floating a => a -> a
cos Double
u1
    sinU2 :: Double
sinU2 = Double -> Double
forall a. Floating a => a -> a
sin Double
u2
    cosU2 :: Double
cosU2 = Double -> Double
forall a. Floating a => a -> a
cos Double
u2

    nextLambda :: Double -> (Double, (Double, Double, Double, Double, Double))
nextLambda Double
lambda = (Double
lambda1, (Double
cos2Alpha, Double
delta, Double
sinDelta, Double
cosDelta, Double
cos2DeltaM))
      where
        sinLambda :: Double
sinLambda = Double -> Double
forall a. Floating a => a -> a
sin Double
lambda
        cosLambda :: Double
cosLambda = Double -> Double
forall a. Floating a => a -> a
cos Double
lambda
        sinDelta :: Double
sinDelta = Double -> Double
forall a. Floating a => a -> a
sqrt((Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLambda) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+
                        (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLambda) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)
        cosDelta :: Double
cosDelta = Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLambda
        delta :: Double
delta = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinDelta Double
cosDelta
        sinAlpha :: Double
sinAlpha = if Double
sinDelta Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLambda Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sinDelta
        cos2Alpha :: Double
cos2Alpha = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
        cos2DeltaM :: Double
cos2DeltaM = if Double
cos2Alpha Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                     then Double
0
                     else Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
cos2Alpha
        c :: Double
c = (Double
fDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
16) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Alpha))
        lambda1 :: Double
lambda1 = Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
c) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha
                  Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDelta
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)))
    lambdas :: [(Double, (Double, Double, Double, Double, Double))]
lambdas = ((Double, (Double, Double, Double, Double, Double))
 -> (Double, (Double, Double, Double, Double, Double)))
-> (Double, (Double, Double, Double, Double, Double))
-> [(Double, (Double, Double, Double, Double, Double))]
forall a. (a -> a) -> a -> [a]
iterate (Double -> (Double, (Double, Double, Double, Double, Double))
nextLambda (Double -> (Double, (Double, Double, Double, Double, Double)))
-> ((Double, (Double, Double, Double, Double, Double)) -> Double)
-> (Double, (Double, Double, Double, Double, Double))
-> (Double, (Double, Double, Double, Double, Double))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, (Double, Double, Double, Double, Double)) -> Double
forall a b. (a, b) -> a
fst) (Double
l, (Double, Double, Double, Double, Double)
forall a. HasCallStack => a
undefined)
    converging :: ((a, b), (a, b)) -> Bool
converging ((a
l1,b
_),(a
l2,b
_)) = a -> a
forall a. Num a => a -> a
abs (a
l1 a -> a -> a
forall a. Num a => a -> a -> a
- a
l2) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
1e-14


-- | Add or subtract multiples of 2*pi so that for all @t@, @-pi < properAngle t < pi@.
properAngle :: Double -> Double
properAngle :: Double -> Double
properAngle Double
t
   | Double
r1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double -> Double
forall a. Num a => a -> a
negate Double
forall a. Floating a => a
pi    = Double
r1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
pi2
   | Double
r1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi            = Double
r1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
pi2
   | Bool
otherwise          = Double
r1
   where
      pf :: Double -> (Int, Double)
      pf :: Double -> (Int, Double)
pf = Double -> (Int, Double)
forall b. Integral b => Double -> (b, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction  -- Shut up GHC warning about defaulting to Integer.
      (Int
_,Double
r) = Double -> (Int, Double)
pf (Double
tDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
pi2)
      r1 :: Double
r1 = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
pi2
      pi2 :: Double
pi2 = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2