-- |

-- Module:      Data.Geo.Jord.Geodesic

-- Copyright:   (c) 2020 Cedric Liegeois

-- License:     BSD3

-- Maintainer:  Cedric Liegeois <ofmooseandmen@yahoo.fr>

-- Stability:   experimental

-- Portability: portable

--

-- Solutions to the direct and inverse geodesic problems on ellipsoidal models using Vincenty formulaes.

-- A geodesic is the shortest path between two points on a curved surface - here an ellispoid. Using these

-- functions improves on the accuracy available using "Data.Geo.Jord.GreatCircle" at the expense of higher

-- CPU usage.

--

-- In order to use this module you should start with the following imports:

--

-- @

-- import qualified Data.Geo.Jord.Angle as Angle

-- import Data.Geo.Jord.Geodesic (Geodesic)

-- import qualified Data.Geo.Jord.Geodesic as Geodesic

-- import qualified Data.Geo.Jord.Geodetic as Geodetic

-- import qualified Data.Geo.Jord.Length as Length

-- @

--

-- <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf T Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations", Survey Review, vol XXIII no 176, 1975.>

module Data.Geo.Jord.Geodesic
    (
    -- * The 'Geodesic' type

      Geodesic
    , startPosition
    , endPosition
    , initialBearing
    , finalBearing
    , length
    -- * Calculations

    , direct
    , inverse
    ) where

import Prelude hiding (length)

import Data.Geo.Jord.Angle (Angle)
import qualified Data.Geo.Jord.Angle as Angle
import Data.Geo.Jord.Ellipsoid
import Data.Geo.Jord.Geodetic (HorizontalPosition)
import qualified Data.Geo.Jord.Geodetic as Geodetic
import Data.Geo.Jord.Length (Length)
import qualified Data.Geo.Jord.Length as Length
import Data.Geo.Jord.Model (Ellipsoidal, surface)

-- | Geodesic line: shortest route between two positions on the surface of a model.

--

-- Bearing are in compass angles.

-- Compass angles are clockwise angles from true north: 0° = north, 90° = east, 180° = south, 270° = west.

-- The final bearing will differ from the initial bearing by varying degrees according to distance and latitude.

data Geodesic a =
    Geodesic
        { Geodesic a -> HorizontalPosition a
startPosition :: HorizontalPosition a -- ^ geodesic start position.

        , Geodesic a -> HorizontalPosition a
endPosition :: HorizontalPosition a -- ^ geodesic end position.

        , Geodesic a -> Maybe Angle
initialBearing :: Maybe Angle -- ^ initial bearing from @startPosition@ to @endPosition@, if both are different.

        , Geodesic a -> Maybe Angle
finalBearing :: Maybe Angle -- ^ final bearing from @startPosition@ to @endPosition@, if both are different.

        , Geodesic a -> Length
length :: Length -- ^ length of the geodesic: the surface distance between @startPosition@ to @endPosition@.

        }
    deriving (Geodesic a -> Geodesic a -> Bool
(Geodesic a -> Geodesic a -> Bool)
-> (Geodesic a -> Geodesic a -> Bool) -> Eq (Geodesic a)
forall a. Model a => Geodesic a -> Geodesic a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Geodesic a -> Geodesic a -> Bool
$c/= :: forall a. Model a => Geodesic a -> Geodesic a -> Bool
== :: Geodesic a -> Geodesic a -> Bool
$c== :: forall a. Model a => Geodesic a -> Geodesic a -> Bool
Eq, Int -> Geodesic a -> ShowS
[Geodesic a] -> ShowS
Geodesic a -> String
(Int -> Geodesic a -> ShowS)
-> (Geodesic a -> String)
-> ([Geodesic a] -> ShowS)
-> Show (Geodesic a)
forall a. Model a => Int -> Geodesic a -> ShowS
forall a. Model a => [Geodesic a] -> ShowS
forall a. Model a => Geodesic a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Geodesic a] -> ShowS
$cshowList :: forall a. Model a => [Geodesic a] -> ShowS
show :: Geodesic a -> String
$cshow :: forall a. Model a => Geodesic a -> String
showsPrec :: Int -> Geodesic a -> ShowS
$cshowsPrec :: forall a. Model a => Int -> Geodesic a -> ShowS
Show)

-- | @direct p1 b1 d@ solves the direct geodesic problem using Vicenty formula: position

-- along the geodesic, reached from position @p1@ having travelled the __surface__ distance @d@ on

-- the initial bearing (compass angle) @b1@ at __constant__ height; it also returns the final bearing

-- at the reached position. For example:

--

-- >>> Geodesic.direct (Geodetic.northPole WGS84) Angle.zero (Length.kilometres 20003.931458623)

-- Just (Geodesic { startPosition = 90°0'0.000"N,0°0'0.000"E (WGS84)

--                , endPosition = 90°0'0.000"S,180°0'0.000"E (WGS84)

--                , initialBearing = Just 0°0'0.000"

--                , finalBearing = Just 180°0'0.000"

--                , length = 20003.931458623km})

--

-- The Vincenty formula for the direct problem should always converge, however this function returns

-- 'Nothing' if it would ever fail to do so (probably thus indicating a bug in the implementation).

direct :: (Ellipsoidal a) => HorizontalPosition a -> Angle -> Length -> Maybe (Geodesic a)
direct :: HorizontalPosition a -> Angle -> Length -> Maybe (Geodesic a)
direct HorizontalPosition a
p1 Angle
b1 Length
d
    | Length
d Length -> Length -> Bool
forall a. Eq a => a -> a -> Bool
== Length
Length.zero = Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p1 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) Length
Length.zero)
    | Bool
otherwise =
        case Maybe (Double, Double, Double, Double)
rec of
            Maybe (Double, Double, Double, Double)
Nothing -> Maybe (Geodesic a)
forall a. Maybe a
Nothing
            (Just (Double
s, Double
cosS, Double
sinS, Double
cos2S')) -> Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b2) Length
d)
                where x :: Double
x = Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1
                      lat2 :: Double
lat2 =
                          Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2
                              (Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1)
                              ((Double
1.0 Double -> 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
sqrt (Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
                      lambda :: Double
lambda = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha1) (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1)
                      _C :: Double
_C = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha))
                      _L :: Double
_L =
                          Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
-
                          (Double
1.0 Double -> 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
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S')))
                      lon2 :: Double
lon2 = Double
lon1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_L
                      b2 :: Angle
b2 =
                          Angle -> Angle -> Angle
Angle.normalise
                              (Double -> Angle
Angle.radians (Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinAlpha (-Double
x)))
                              (Double -> Angle
Angle.decimalDegrees Double
360.0)
                      p2 :: HorizontalPosition a
p2 =
                          Angle -> Angle -> a -> HorizontalPosition a
forall a. Model a => Angle -> Angle -> a -> HorizontalPosition a
Geodetic.latLongPos'
                              (Double -> Angle
Angle.radians Double
lat2)
                              (Double -> Angle
Angle.radians Double
lon2)
                              (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p1)
  where
    lat1 :: Double
lat1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    lon1 :: Double
lon1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    ell :: Ellipsoid
ell = a -> Ellipsoid
forall a. Model a => a -> Ellipsoid
surface (a -> Ellipsoid)
-> (HorizontalPosition a -> a) -> HorizontalPosition a -> Ellipsoid
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model (HorizontalPosition a -> Ellipsoid)
-> HorizontalPosition a -> Ellipsoid
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    (Double
a, Double
b, Double
f) = Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
ell
    br1 :: Double
br1 = Angle -> Double
Angle.toRadians Angle
b1
    cosAlpha1 :: Double
cosAlpha1 = Double -> Double
forall a. Floating a => a -> a
cos Double
br1
    sinAlpha1 :: Double
sinAlpha1 = Double -> Double
forall a. Floating a => a -> a
sin Double
br1
    (Double
tanU1, Double
cosU1, Double
sinU1) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat1 Double
f
    sigma1 :: Double
sigma1 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
tanU1 Double
cosAlpha1 -- angular distance on the sphere from the equator to p1

    sinAlpha :: Double
sinAlpha = Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha1 -- alpha = azimuth of the geodesic at the equator

    cosSqAlpha :: Double
cosSqAlpha = Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha
    uSq :: Double
uSq = Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
    _A :: Double
_A = Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16384.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
768.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
320.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
    _B :: Double
_B = Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1024.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
128.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
74.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
47.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
    dm :: Double
dm = Length -> Double
Length.toMetres Length
d
    sigma :: Double
sigma = Double
dm Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A)
    rec :: Maybe (Double, Double, Double, Double)
rec = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dm Double
_A Double
_B Double
b Double
sigma Int
0

-- | @inverse p1 p2@ solves the inverse geodesic problem using Vicenty formula: __surface__ distance,

-- and initial/final bearing between the geodesic line between positions @p1@ and @p2@. For example:

--

-- >>> Geodesic.inverse (Geodetic.latLongPos 0 0 WGS84) (Geodetic.latLongPos 0.5 179.5 WGS84)

-- Just (Geodesic { startPosition = 0°0'0.000"N,0°0'0.000"E 0.0m (WGS84)

--                , endPosition = 0°30'0.000"N,179°30'0.000"E 0.0m (WGS84)

--                , initialBearing = Just 25°40'18.742"

--                , finalBearing = Just 154°19'37.507"

--                , length = 19936.288578981km})

--

-- The Vincenty formula for the inverse problem can fail to converge for nearly antipodal points in which

-- case this function returns 'Nothing'. For example:

--

-- >>> Geodesic.inverse (Geodetic.latLongPos 0 0 WGS84) (Geodetic.latLongPos 0.5 179.7 WGS84)

-- Nothing

inverse :: (Ellipsoidal a) => HorizontalPosition a -> HorizontalPosition a -> Maybe (Geodesic a)
inverse :: HorizontalPosition a -> HorizontalPosition a -> Maybe (Geodesic a)
inverse HorizontalPosition a
p1 HorizontalPosition a
p2
    | HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 Maybe Angle
forall a. Maybe a
Nothing Maybe Angle
forall a. Maybe a
Nothing Length
Length.zero)
    | Bool
otherwise =
        case Maybe
  (Double, Double, Double, Double, Double, Double, Double, Double)
rec of
            Maybe
  (Double, Double, Double, Double, Double, Double, Double, Double)
Nothing -> Maybe (Geodesic a)
forall a. Maybe a
Nothing
            (Just (Double
cosL, Double
sinL, Double
s, Double
cosS, Double
sinS, Double
sinSqS, Double
cos2S', Double
cosSqA)) ->
                Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b2) Length
d)
                where uSq :: Double
uSq = Double
cosSqA Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
                      _A :: Double
_A =
                          Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+
                          Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16384.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
768.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
320.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
                      _B :: Double
_B = Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1024.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
128.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
74.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
47.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
                      deltaSigma :: Double
deltaSigma =
                          Double
_B Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                          (Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
+
                           Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                           (Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S') Double -> Double -> Double
forall a. Num a => a -> a -> a
-
                            Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS) Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                            (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S')))
                      d :: Length
d = Double -> Length
Length.metres (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
deltaSigma))
                      a1R :: Double
a1R =
                          if Double -> Double
forall a. Num a => a -> a
abs Double
sinSqS Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
                              then Double
0.0
                              else Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) (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
cosL)
                      a2R :: Double
a2R =
                          if Double -> Double
forall a. Num a => a -> a
abs Double
sinSqS Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
                              then Double
forall a. Floating a => a
pi
                              else Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) (-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
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL)
                      b1 :: Angle
b1 = Angle -> Angle -> Angle
Angle.normalise (Double -> Angle
Angle.radians Double
a1R) (Double -> Angle
Angle.decimalDegrees Double
360.0)
                      b2 :: Angle
b2 = Angle -> Angle -> Angle
Angle.normalise (Double -> Angle
Angle.radians Double
a2R) (Double -> Angle
Angle.decimalDegrees Double
360.0)
  where
    lat1 :: Double
lat1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    lon1 :: Double
lon1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    lat2 :: Double
lat2 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p2
    lon2 :: Double
lon2 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p2
    ell :: Ellipsoid
ell = a -> Ellipsoid
forall a. Model a => a -> Ellipsoid
surface (a -> Ellipsoid)
-> (HorizontalPosition a -> a) -> HorizontalPosition a -> Ellipsoid
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model (HorizontalPosition a -> Ellipsoid)
-> HorizontalPosition a -> Ellipsoid
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
    (Double
a, Double
b, Double
f) = Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
ell
    _L :: Double
_L = Double
lon2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lon1 -- difference in longitude

    (Double
_, Double
cosU1, Double
sinU1) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat1 Double
f
    (Double
_, Double
cosU2, Double
sinU2) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat2 Double
f
    antipodal :: Bool
antipodal = Double -> Double
forall a. Num a => a -> a
abs Double
_L Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs (Double
lat2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lat1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0
    rec :: Maybe
  (Double, Double, Double, Double, Double, Double, Double, Double)
rec = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
     (Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
_L Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal Int
0

directRec ::
       Double
    -> Double
    -> Double
    -> Double
    -> Double
    -> Double
    -> Int
    -> Maybe (Double, Double, Double, Double)
directRec :: Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dist Double
_A Double
_B Double
b Double
sigma Int
i
    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
100 = Maybe (Double, Double, Double, Double)
forall a. Maybe a
Nothing
    | Double -> Double
forall a. Num a => a -> a
abs (Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
newSigma) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-12 = (Double, Double, Double, Double)
-> Maybe (Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
newSigma, Double
cosSigma, Double
sinSigma, Double
cos2Sigma')
    | Bool
otherwise = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dist Double
_A Double
_B Double
b Double
newSigma (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
  where
    cos2Sigma' :: Double
cos2Sigma' = Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sigma1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sigma)
    sinSigma :: Double
sinSigma = Double -> Double
forall a. Floating a => a -> a
sin Double
sigma
    cosSigma :: Double
cosSigma = Double -> Double
forall a. Floating a => a -> a
cos Double
sigma
    deltaSigma :: Double
deltaSigma =
        Double
_B Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
*
        (Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
+
         Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
*
         (Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma') Double -> Double -> Double
forall a. Num a => a -> a -> a
-
          Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma) Double -> Double -> Double
forall a. Num a => a -> a -> a
*
          (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma')))
    newSigma :: Double
newSigma = Double
dist Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
deltaSigma

inverseRec ::
       Double
    -> Double
    -> Double
    -> Double
    -> Double
    -> Double
    -> Double
    -> Bool
    -> Int
    -> Maybe (Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec :: Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
     (Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
lambda Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal Int
i
    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1000 = Maybe
  (Double, Double, Double, Double, Double, Double, Double, Double)
forall a. Maybe a
Nothing
    -- co-incident/antipodal points (falls back on λ/σ = L)

    | Double
sinSqSigma Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon = (Double, Double, Double, Double, Double, Double, Double, Double)
-> Maybe
     (Double, Double, Double, Double, Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
-> Double
-> Double
-> Bool
-> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback Double
cosL Double
sinL Double
sinSqSigma Bool
antipodal)
    | Double
iterationCheck Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi = Maybe
  (Double, Double, Double, Double, Double, Double, Double, Double)
forall a. Maybe a
Nothing
    | Double -> Double
forall a. Num a => a -> a
abs (Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
newLambda) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-12 =
        (Double, Double, Double, Double, Double, Double, Double, Double)
-> Maybe
     (Double, Double, Double, Double, Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
cosL, Double
sinL, Double
sigma, Double
cosSigma, Double
sinSigma, Double
sinSqSigma, Double
cos2Sigma', Double
cosSqAlpha)
    | Bool
otherwise = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
     (Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
newLambda Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
  where
    sinL :: Double
sinL = Double -> Double
forall a. Floating a => a -> a
sin Double
lambda
    cosL :: Double
cosL = Double -> Double
forall a. Floating a => a -> a
cos Double
lambda
    sinSqSigma :: Double
sinSqSigma =
        (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) 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
cosL) 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
cosL)
    sinSigma :: Double
sinSigma = Double -> Double
forall a. Floating a => a -> a
sqrt Double
sinSqSigma
    cosSigma :: Double
cosSigma = 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
cosL
    sigma :: Double
sigma = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinSigma Double
cosSigma
    sinAlpha :: Double
sinAlpha = 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
sinL Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sinSigma
    cosSqAlpha :: Double
cosSqAlpha = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha
    cos2Sigma' :: Double
cos2Sigma' =
        if Double
cosSqAlpha Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0
            then Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2.0 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
cosSqAlpha
            else Double
0
    _C :: Double
_C = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha))
    newLambda :: Double
newLambda =
        Double
_L Double -> Double -> Double
forall a. Num a => a -> a -> a
+
        (Double
1.0 Double -> 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
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
+
         Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma')))
    iterationCheck :: Double
iterationCheck =
        if Bool
antipodal
            then Double -> Double
forall a. Num a => a -> a
abs Double
newLambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
pi
            else Double -> Double
forall a. Num a => a -> a
abs Double
newLambda

inverseFallback ::
       Double
    -> Double
    -> Double
    -> Bool
    -> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback :: Double
-> Double
-> Double
-> Bool
-> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback Double
cosL Double
sinL Double
sinSqSigma Bool
antipodal =
    (Double
cosL, Double
sinL, Double
sigma, Double
cosSigma, Double
sinSigma, Double
sinSqSigma, Double
cos2Sigma', Double
cosSqAlpha)
  where
    sigma :: Double
sigma =
        if Bool
antipodal
            then Double
forall a. Floating a => a
pi
            else Double
0
    cosSigma :: Double
cosSigma =
        if Bool
antipodal
            then (-Double
1)
            else Double
1
    sinSigma :: Double
sinSigma = Double
0
    cos2Sigma' :: Double
cos2Sigma' = Double
1
    cosSqAlpha :: Double
cosSqAlpha = Double
1

-- | see Numeric.Limits

epsilon :: Double
epsilon :: Double
epsilon = Double
r
  where
    r :: Double
r = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat (Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) Int
e
    (Integer
m, Int
e) = Double -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat (Double
1 :: Double)

reducedLat :: Double -> Double -> (Double, Double, Double)
reducedLat :: Double -> Double -> (Double, Double, Double)
reducedLat Double
lat Double
f = (Double
tanU, Double
cosU, Double
sinU)
  where
    tanU :: Double
tanU = (Double
1.0 Double -> 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 Double
lat
    cosU :: Double
cosU = Double
1.0 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
+ Double
tanU Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanU)
    sinU :: Double
sinU = Double
tanU Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU

abf :: Ellipsoid -> (Double, Double, Double)
abf :: Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
e = (Length -> Double
Length.toMetres (Length -> Double) -> (Ellipsoid -> Length) -> Ellipsoid -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ellipsoid -> Length
equatorialRadius (Ellipsoid -> Double) -> Ellipsoid -> Double
forall a b. (a -> b) -> a -> b
$ Ellipsoid
e, Length -> Double
Length.toMetres (Length -> Double) -> (Ellipsoid -> Length) -> Ellipsoid -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ellipsoid -> Length
polarRadius (Ellipsoid -> Double) -> Ellipsoid -> Double
forall a b. (a -> b) -> a -> b
$ Ellipsoid
e, Ellipsoid -> Double
flattening Ellipsoid
e)