-- |

-- Module:      Data.Geo.Jord.GreatCircle

-- Copyright:   (c) 2020 Cedric Liegeois

-- License:     BSD3

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

-- Stability:   experimental

-- Portability: portable

--

-- Geographical Position calculations on great circles, i.e. using a __sphere__ to represent

-- the celestial body that positions refer to.

--

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

--

-- @

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

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

-- @

--

-- All functions are implemented using the vector-based approached described in

-- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Point_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation>

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

      GreatCircle
    , through
    , headingOn
    -- * The 'MinorArc' type

    , MinorArc
    , minorArc
    , minorArcStart
    , minorArcEnd
    -- * Calculations

    , Side(..)
    , alongTrackDistance
    , alongTrackDistance'
    , angularDistance
    , crossTrackDistance
    , crossTrackDistance'
    , destination
    , distance
    , enclosedBy
    , finalBearing
    , initialBearing
    , interpolated
    , intersection
    , intersections
    , mean
    , projection
    , side
    , turn
    ) where

import Data.Geo.Jord.Angle (Angle)
import qualified Data.Geo.Jord.Angle as Angle
import Data.Geo.Jord.Ellipsoid (equatorialRadius)
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.Math3d (V3)
import qualified Data.Geo.Jord.Math3d as Math3d
import Data.Geo.Jord.Model (Spherical, surface)

-- | A circle on the surface of a __sphere__ which lies in a plane

-- passing through the sphere centre. Every two distinct and non-antipodal points

-- define a unique Great Circle.

--

-- It is internally represented as its normal vector - i.e. the normal vector

-- to the plane containing the great circle.

data GreatCircle a =
    GreatCircle !V3 !a String
    deriving (GreatCircle a -> GreatCircle a -> Bool
(GreatCircle a -> GreatCircle a -> Bool)
-> (GreatCircle a -> GreatCircle a -> Bool) -> Eq (GreatCircle a)
forall a. Eq a => GreatCircle a -> GreatCircle a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: GreatCircle a -> GreatCircle a -> Bool
$c/= :: forall a. Eq a => GreatCircle a -> GreatCircle a -> Bool
== :: GreatCircle a -> GreatCircle a -> Bool
$c== :: forall a. Eq a => GreatCircle a -> GreatCircle a -> Bool
Eq)

instance (Spherical a) => Show (GreatCircle a) where
    show :: GreatCircle a -> String
show (GreatCircle V3
_ a
_ String
s) = String
s

-- | @through p1 p2@ returns the 'GreatCircle' passing by both positions @p1@ and @p2@ (in this direction). For example:

--

-- >>> let p1 = Geodetic.s84Pos 45.0 (-143.5)

-- >>> let p2 = Geodetic.s84Pos 46.0 14.5

-- >>> GreatCircle.through p1 p2

-- Just Great Circle { through 45°0'0.000"N,143°30'0.000"W (S84) & 46°0'0.000"N,14°30'0.000"E (S84) }

--

-- Returns 'Nothing' if given positions are equal or @p1@ is antipode of @p2@.

through :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Maybe (GreatCircle a)
through :: HorizontalPosition a
-> HorizontalPosition a -> Maybe (GreatCircle a)
through HorizontalPosition a
p1 HorizontalPosition a
p2 = (V3 -> GreatCircle a) -> Maybe V3 -> Maybe (GreatCircle a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\V3
n -> V3 -> a -> String -> GreatCircle a
forall a. V3 -> a -> String -> GreatCircle a
GreatCircle V3
n (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p1) String
dscr) (HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
p1 HorizontalPosition a
p2)
  where
    dscr :: String
dscr = String
"Great Circle { through " String -> ShowS
forall a. [a] -> [a] -> [a]
++ HorizontalPosition a -> String
forall a. Show a => a -> String
show HorizontalPosition a
p1 String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" & " String -> ShowS
forall a. [a] -> [a] -> [a]
++ HorizontalPosition a -> String
forall a. Show a => a -> String
show HorizontalPosition a
p2 String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" }"

-- | @headingOn p b@ returns the 'GreatCircle' passing by position @p@ and heading on bearing @b@. For example:

--

-- >>> let p = Geodetic.s84Pos 45.0 (-143.5)

-- >>> let b = Angle.decimalDegrees 33.0

-- >>> GreatCircle.headingOn p b

-- Great Circle { by 45°0'0.000"N,143°30'0.000"W (S84) & heading on 33°0'0.000" }

headingOn :: (Spherical a) => HorizontalPosition a -> Angle -> GreatCircle a
headingOn :: HorizontalPosition a -> Angle -> GreatCircle a
headingOn HorizontalPosition a
p Angle
b = V3 -> a -> String -> GreatCircle a
forall a. V3 -> a -> String -> GreatCircle a
GreatCircle (V3 -> V3 -> V3
Math3d.subtract V3
n' V3
e') a
m String
dscr
  where
    m :: a
m = HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p
    v :: V3
v = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p
    e :: V3
e = V3 -> V3 -> V3
Math3d.cross (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector (HorizontalPosition a -> V3)
-> (a -> HorizontalPosition a) -> a -> V3
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> HorizontalPosition a
forall a. Model a => a -> HorizontalPosition a
Geodetic.northPole (a -> V3) -> a -> V3
forall a b. (a -> b) -> a -> b
$ a
m) V3
v -- easting

    n :: V3
n = V3 -> V3 -> V3
Math3d.cross V3
v V3
e -- northing

    e' :: V3
e' = V3 -> Double -> V3
Math3d.scale V3
e (Angle -> Double
Angle.cos Angle
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ V3 -> Double
Math3d.norm V3
e)
    n' :: V3
n' = V3 -> Double -> V3
Math3d.scale V3
n (Angle -> Double
Angle.sin Angle
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ V3 -> Double
Math3d.norm V3
n)
    dscr :: String
dscr = String
"Great Circle { by " String -> ShowS
forall a. [a] -> [a] -> [a]
++ HorizontalPosition a -> String
forall a. Show a => a -> String
show HorizontalPosition a
p String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" & heading on " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Angle -> String
forall a. Show a => a -> String
show Angle
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" }"

-- | Oriented minor arc of a great circle between two positions: shortest path between positions on a great circle.

data MinorArc a =
    MinorArc !V3 (HorizontalPosition a) (HorizontalPosition a)
    deriving (MinorArc a -> MinorArc a -> Bool
(MinorArc a -> MinorArc a -> Bool)
-> (MinorArc a -> MinorArc a -> Bool) -> Eq (MinorArc a)
forall a. Model a => MinorArc a -> MinorArc a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MinorArc a -> MinorArc a -> Bool
$c/= :: forall a. Model a => MinorArc a -> MinorArc a -> Bool
== :: MinorArc a -> MinorArc a -> Bool
$c== :: forall a. Model a => MinorArc a -> MinorArc a -> Bool
Eq)

instance (Spherical a) => Show (MinorArc a) where
    show :: MinorArc a -> String
show (MinorArc V3
_ HorizontalPosition a
s HorizontalPosition a
e) = String
"Minor Arc { from: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ HorizontalPosition a -> String
forall a. Show a => a -> String
show HorizontalPosition a
s String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
", to: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ HorizontalPosition a -> String
forall a. Show a => a -> String
show HorizontalPosition a
e String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" }"

-- | @minorArc p1 p2@ returns the 'MinorArc' from @p1@ to @p2@.  For example:

--

-- >>> let p1 = Geodetic.s84Pos 45.0 (-143.5)

-- >>> let p2 = Geodetic.s84Pos 46.0 14.5

-- >>> GreatCircle.minorArc p1 p2

-- Just Minor Arc { from: 45°0'0.000"N,143°30'0.000"W (S84), to: 46°0'0.000"N,14°30'0.000"E (S84) }

--

-- Returns 'Nothing' if given positions are equal.

minorArc :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
minorArc :: HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
minorArc HorizontalPosition a
p1 HorizontalPosition a
p2 = (V3 -> MinorArc a) -> Maybe V3 -> Maybe (MinorArc a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\V3
n -> V3 -> HorizontalPosition a -> HorizontalPosition a -> MinorArc a
forall a.
V3 -> HorizontalPosition a -> HorizontalPosition a -> MinorArc a
MinorArc V3
n HorizontalPosition a
p1 HorizontalPosition a
p2) (HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
p1 HorizontalPosition a
p2)

-- | @minorArcStart ma@ returns the start position of minor arc @ma@.

minorArcStart :: (Spherical a) => MinorArc a -> HorizontalPosition a
minorArcStart :: MinorArc a -> HorizontalPosition a
minorArcStart (MinorArc V3
_ HorizontalPosition a
s HorizontalPosition a
_) = HorizontalPosition a
s

-- | @minorArcEnd ma@ returns the end position of minor arc @ma@.

minorArcEnd :: (Spherical a) => MinorArc a -> HorizontalPosition a
minorArcEnd :: MinorArc a -> HorizontalPosition a
minorArcEnd (MinorArc V3
_ HorizontalPosition a
_ HorizontalPosition a
e) = HorizontalPosition a
e

-- | Side of a position w.r.t. to a great circle.

data Side
    = LeftOf -- ^ position is left of the great circle.

    | RightOf -- ^ position is right of the great circle.

    | None -- ^ position is on the great circle, or the great circle is undefined.

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

-- | @alongTrackDistance p a@ computes how far position @p@ is along a path described by the minor arc @a@: if a

-- perpendicular is drawn from @p@  to the path, the along-track distance is the signed distance from the start point to

-- where the perpendicular crosses the path. For example:

--

-- >>> let p = Geodetic.s84Pos 53.2611 (-0.7972)

-- >>> let mas = Geodetic.s84Pos 53.3206 (-1.7297)

-- >>> let mae = Geodetic.s84Pos 53.1887 0.1334

-- >>> fmap (GreatCircle.alongTrackDistance p) (GreatCircle.minorArc mas mae)

-- Just 62.3315791km

alongTrackDistance :: (Spherical a) => HorizontalPosition a -> MinorArc a -> Length
alongTrackDistance :: HorizontalPosition a -> MinorArc a -> Length
alongTrackDistance HorizontalPosition a
p (MinorArc V3
n HorizontalPosition a
s HorizontalPosition a
_) = HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
alongTrackDistance'' HorizontalPosition a
p HorizontalPosition a
s V3
n

-- | @alongTrackDistance' p s b@ computes how far Position @p@ is along a path starting at @s@ and heading on

-- bearing @b@: if a perpendicular is drawn from @p@  to the path, the along-track distance is the signed distance from

-- the start point to where the perpendicular crosses the path. For example:

--

-- >>> let p = Geodetic.s84Pos 53.2611 (-0.7972)

-- >>> let s = Geodetic.s84Pos 53.3206 (-1.7297)

-- >>> let b = Angle.decimalDegrees 96.0017325

-- >>> GreatCircle.alongTrackDistance' p s b

-- 62.3315791km

alongTrackDistance' ::
       (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Angle -> Length
alongTrackDistance' :: HorizontalPosition a -> HorizontalPosition a -> Angle -> Length
alongTrackDistance' HorizontalPosition a
p HorizontalPosition a
s Angle
b = HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
alongTrackDistance'' HorizontalPosition a
p HorizontalPosition a
s V3
n
  where
    (GreatCircle V3
n a
_ String
_) = HorizontalPosition a -> Angle -> GreatCircle a
forall a.
Spherical a =>
HorizontalPosition a -> Angle -> GreatCircle a
headingOn HorizontalPosition a
s Angle
b

-- | @angularDistance p1 p2 n@ computes the angle between the horizontal Points @p1@ and @p2@.

-- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], signed + if @p1@ is clockwis

-- looking along @n@, - in opposite direction.

angularDistance ::
       (Spherical a)
    => HorizontalPosition a
    -> HorizontalPosition a
    -> Maybe (HorizontalPosition a)
    -> Angle
angularDistance :: HorizontalPosition a
-> HorizontalPosition a -> Maybe (HorizontalPosition a) -> Angle
angularDistance HorizontalPosition a
p1 HorizontalPosition a
p2 Maybe (HorizontalPosition a)
n = V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween V3
v1 V3
v2 Maybe V3
vn
  where
    v1 :: V3
v1 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p1
    v2 :: V3
v2 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p2
    vn :: Maybe V3
vn = (HorizontalPosition a -> V3)
-> Maybe (HorizontalPosition a) -> Maybe V3
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector Maybe (HorizontalPosition a)
n

-- | @crossTrackDistance p gc@ computes the signed distance from horizontal Position @p@ to great circle @gc@.

-- Returns a negative 'Length' if Position is left of great circle, positive 'Length' if Position is right

-- of great circle; the orientation of the great circle is therefore important. For example:

--

-- >>> let p = Geodetic.s84Pos 53.2611 (-0.7972)

-- >>> let gc1 = GreatCircle.through (Geodetic.s84Pos 51 0) (Geodetic.s84Pos 52 1)

-- >>> fmap (GreatCircle.crossTrackDistance p) gc1

-- Just -176.756870526km

-- >>> let gc2 = GreatCircle.through (Geodetic.s84Pos 52 1) (Geodetic.s84Pos 51 0)

-- >>> fmap (GreatCircle.crossTrackDistance p) gc2

-- Just 176.7568725km

-- >>> let gc3 = GreatCircle.headingOn (Geodetic.s84Pos 53.3206 (-1.7297)) (Angle.decimalDegrees 96.0)

-- >>> GreatCircle.crossTrackDistance p gc3

-- -305.665267m metres

crossTrackDistance :: (Spherical a) => HorizontalPosition a -> GreatCircle a -> Length
crossTrackDistance :: HorizontalPosition a -> GreatCircle a -> Length
crossTrackDistance HorizontalPosition a
p (GreatCircle V3
n a
_ String
_) =
    Angle -> Length -> Length
Angle.arcLength (Angle -> Angle -> Angle
Angle.subtract Angle
a (Double -> Angle
Angle.decimalDegrees Double
90)) (HorizontalPosition a -> Length
forall a. Spherical a => HorizontalPosition a -> Length
radius HorizontalPosition a
p)
  where
    a :: Angle
a = V3 -> V3 -> Angle
angleBetween V3
n (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p)

-- | @crossTrackDistance' p s b@ computes the signed distance from horizontal Position @p@ to the great circle passing

-- by @s@ and heading on bearing @b@.

--

-- This is equivalent to:

--

-- > GreatCircle.crossTrackDistance p (GreatCircle.headingOn s b)

crossTrackDistance' ::
       (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Angle -> Length
crossTrackDistance' :: HorizontalPosition a -> HorizontalPosition a -> Angle -> Length
crossTrackDistance' HorizontalPosition a
p HorizontalPosition a
s Angle
b = HorizontalPosition a -> GreatCircle a -> Length
forall a.
Spherical a =>
HorizontalPosition a -> GreatCircle a -> Length
crossTrackDistance HorizontalPosition a
p (HorizontalPosition a -> Angle -> GreatCircle a
forall a.
Spherical a =>
HorizontalPosition a -> Angle -> GreatCircle a
headingOn HorizontalPosition a
s Angle
b)

-- | @destination p b d@ computes the position along the great circle, reached from position @p@ having travelled the

-- distance @d@ on the initial bearing (compass angle) @b@. For example:

--

-- >>> let p = Geodetic.s84Pos 54 154

-- >>> GreatCircle.destination p (Angle.decimalDegrees 33) (Length.kilometres 1000)

-- 61°10'44.188"N,164°10'19.254"E (S84)

--

-- Note that the bearing will normally vary before destination is reached.

destination :: (Spherical a) => HorizontalPosition a -> Angle -> Length -> HorizontalPosition a
destination :: HorizontalPosition a -> Angle -> Length -> HorizontalPosition a
destination HorizontalPosition a
p Angle
b Length
d
    | Length
d Length -> Length -> Bool
forall a. Eq a => a -> a -> Bool
== Length
Length.zero = HorizontalPosition a
p
    | Bool
otherwise = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' V3
nvd a
m
  where
    m :: a
m = HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p
    nv :: V3
nv = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p
    ed :: V3
ed = V3 -> V3
Math3d.unit (V3 -> V3 -> V3
Math3d.cross (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector (HorizontalPosition a -> V3)
-> (a -> HorizontalPosition a) -> a -> V3
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> HorizontalPosition a
forall a. Model a => a -> HorizontalPosition a
Geodetic.northPole (a -> V3) -> a -> V3
forall a b. (a -> b) -> a -> b
$ a
m) V3
nv) -- east direction vector at v

    nd :: V3
nd = V3 -> V3 -> V3
Math3d.cross V3
nv V3
ed -- north direction vector at v

    r :: Length
r = HorizontalPosition a -> Length
forall a. Spherical a => HorizontalPosition a -> Length
radius HorizontalPosition a
p
    ta :: Angle
ta = Length -> Length -> Angle
Angle.central Length
d Length
r -- central angle

    de :: V3
de = V3 -> V3 -> V3
Math3d.add (V3 -> Double -> V3
Math3d.scale V3
nd (Angle -> Double
Angle.cos Angle
b)) (V3 -> Double -> V3
Math3d.scale V3
ed (Angle -> Double
Angle.sin Angle
b)) -- unit vector in the direction of the azimuth

    nvd :: V3
nvd = V3 -> V3 -> V3
Math3d.add (V3 -> Double -> V3
Math3d.scale V3
nv (Angle -> Double
Angle.cos Angle
ta)) (V3 -> Double -> V3
Math3d.scale V3
de (Angle -> Double
Angle.sin Angle
ta))

-- | @distance p1 p2@ computes the surface distance on the great circle between the positions @p1@ and @p2@. For example:

--

-- >>> GreatCircle.distance (Geodetic.northPole S84) (Geodetic.southPole S84)

-- 20015.114352233km

-- >>> GreatCircle.distance (Geodetic.northPole S84) (Geodetic.northPole S84)

-- 0.0m

distance :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Length
distance :: HorizontalPosition a -> HorizontalPosition a -> Length
distance HorizontalPosition a
p1 HorizontalPosition a
p2 = Angle -> Length -> Length
Angle.arcLength Angle
a (HorizontalPosition a -> Length
forall a. Spherical a => HorizontalPosition a -> Length
radius HorizontalPosition a
p1)
  where
    a :: Angle
a = V3 -> V3 -> Angle
angleBetween (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p1) (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p2)

-- | @enclosedBy p ps@ determines whether position @p@ is enclosed by the polygon defined by horizontal positions @ps@.

-- The polygon can be opened or closed (i.e. if @head ps /= last ps@).

--

-- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles

-- will sum to less than 360°, and exterior point angles will be small but non-zero.

--

-- Always returns 'False' if @ps@ does not at least defines a triangle or if @p@ is any of the @ps@.

--

-- ==== __Examples__

--

-- >>> let malmo = Geodetic.s84Pos 55.6050 13.0038

-- >>> let ystad = Geodetic.s84Pos 55.4295 13.82

-- >>> let lund = Geodetic.s84Pos 55.7047 13.1910

-- >>> let helsingborg = Geodetic.s84Pos 56.0465 12.6945

-- >>> let kristianstad = Geodetic.s84Pos 56.0294 14.1567

-- >>> let polygon = [malmo, ystad, kristianstad, helsingborg, lund]

-- >>> let hoor = Geodetic.s84Pos 55.9295 13.5297

-- >>> let hassleholm = Geodetic.s84Pos 56.1589 13.7668

-- >>> GreatCircle.enclosedBy hoor polygon

-- True

-- >>> GreatCircle.enclosedBy hassleholm polygon

-- False

enclosedBy :: (Spherical a) => HorizontalPosition a -> [HorizontalPosition a] -> Bool
enclosedBy :: HorizontalPosition a -> [HorizontalPosition a] -> Bool
enclosedBy HorizontalPosition a
p [HorizontalPosition a]
ps
    | [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [HorizontalPosition a]
ps = Bool
False
    | HorizontalPosition a
p HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [HorizontalPosition a]
ps = Bool
False
    | [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head [HorizontalPosition a]
ps HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
last [HorizontalPosition a]
ps = HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall a.
Spherical a =>
HorizontalPosition a -> [HorizontalPosition a] -> Bool
enclosedBy HorizontalPosition a
p ([HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
init [HorizontalPosition a]
ps)
    | [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = Bool
False
    | Bool
otherwise =
        let aSum :: Angle
aSum =
                (Angle -> (V3, V3) -> Angle) -> Angle -> [(V3, V3)] -> Angle
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                    (\Angle
a (V3, V3)
v' -> Angle -> Angle -> Angle
Angle.add Angle
a ((V3 -> V3 -> Maybe V3 -> Angle) -> (V3, V3) -> Maybe V3 -> Angle
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween (V3, V3)
v' (V3 -> Maybe V3
forall a. a -> Maybe a
Just V3
v)))
                    Angle
Angle.zero
                    ([V3] -> [(V3, V3)]
egdes ((V3 -> V3) -> [V3] -> [V3]
forall a b. (a -> b) -> [a] -> [b]
map (V3 -> V3 -> V3
Math3d.subtract V3
v) [V3]
vs))
         in Double -> Double
forall a. Num a => a -> a
abs (Angle -> Double
Angle.toDecimalDegrees Angle
aSum) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
180.0
  where
    v :: V3
v = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p
    vs :: [V3]
vs = (HorizontalPosition a -> V3) -> [HorizontalPosition a] -> [V3]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector [HorizontalPosition a]
ps

-- | @finalBearing p1 p2@ computes the final bearing arriving at @p2@ from @p1@ in compass angle.

-- 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.

-- For example:

--

-- >>> let p1 = Geodetic.s84Pos 0 1

-- >>> let p2 = Geodetic.s84Pos 0 0

-- >>> GreatCircle.finalBearing p1 p2

-- Just 270°0'0.000"

--

-- Returns 'Nothing' if both positions are equals.

finalBearing :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Maybe Angle
finalBearing :: HorizontalPosition a -> HorizontalPosition a -> Maybe Angle
finalBearing HorizontalPosition a
p1 HorizontalPosition a
p2
    | HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Maybe Angle
forall a. Maybe a
Nothing
    | Bool
otherwise = Angle -> Maybe Angle
forall a. a -> Maybe a
Just (Angle -> Angle -> Angle
Angle.normalise (HorizontalPosition a -> HorizontalPosition a -> Angle
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Angle
initialBearing' HorizontalPosition a
p2 HorizontalPosition a
p1) (Double -> Angle
Angle.decimalDegrees Double
180))

-- | @initialBearing p1 p2@ computes the initial bearing from @p1@ to @p2@ in compass angle.

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

-- For example:

--

-- >>> let p1 = Geodetic.s84Pos 58.643889 (-5.714722)

-- >>> let p2 = Geodetic.s84Pos 50.066389 (-5.714722)

-- >>> GreatCircle.initialBearing p1 p2

-- Just 180°0'0.000"

--

-- Returns 'Nothing' if both positions are equals.

initialBearing :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Maybe Angle
initialBearing :: HorizontalPosition a -> HorizontalPosition a -> Maybe Angle
initialBearing HorizontalPosition a
p1 HorizontalPosition a
p2
    | HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Maybe Angle
forall a. Maybe a
Nothing
    | Bool
otherwise = Angle -> Maybe Angle
forall a. a -> Maybe a
Just (HorizontalPosition a -> HorizontalPosition a -> Angle
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Angle
initialBearing' HorizontalPosition a
p1 HorizontalPosition a
p2)

-- | @interpolated p0 p1 f# computes the position at fraction @f@ between the @p0@ and @p1@. For example:

--

-- >>> let p1 = Geodetic.s84Pos 53.479444 (-2.245278)

-- >>> let p2 = Geodetic.s84Pos 55.605833 13.035833

-- >>> GreatCircle.interpolated p1 p2 0.5

-- 54°47'0.805"N,5°11'41.947"E (S84)

--

-- Special conditions:

--

-- @

-- interpolated p0 p1 0.0 = p0

-- interpolated p0 p1 1.0 = p1

-- 'error's if @f < 0 || f > 1@

-- @

interpolated ::
       (Spherical a)
    => HorizontalPosition a
    -> HorizontalPosition a
    -> Double
    -> HorizontalPosition a
interpolated :: HorizontalPosition a
-> HorizontalPosition a -> Double -> HorizontalPosition a
interpolated HorizontalPosition a
p0 HorizontalPosition a
p1 Double
f
    | Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 = String -> HorizontalPosition a
forall a. HasCallStack => String -> a
error (String
"fraction must be in range [0..1], was " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
f)
    | Double
f Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = HorizontalPosition a
p0
    | Double
f Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1 = HorizontalPosition a
p1
    | Bool
otherwise = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' V3
iv (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p0)
  where
    nv0 :: V3
nv0 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p0
    nv1 :: V3
nv1 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p1
    iv :: V3
iv = V3 -> V3
Math3d.unit (V3 -> V3 -> V3
Math3d.add V3
nv0 (V3 -> Double -> V3
Math3d.scale (V3 -> V3 -> V3
Math3d.subtract V3
nv1 V3
nv0) Double
f))

-- | Computes the intersection between the two given minor arcs of great circle. For example:

--

-- >>> let a1s = Geodetic.s84Pos 51.885 0.235

-- >>> let a1e = Geodetic.s84Pos 48.269 13.093

-- >>> let a2s = Geodetic.s84Pos 49.008 2.549

-- >>> let a2e = Geodetic.s84Pos 56.283 11.304

-- >>> GreatCircle.intersection <$> (GreatCircle.minorArc a1s a1e) <*> (GreatCircle.minorArc a2s a2e)

-- Just (Just 50°54'6.260"N,4°29'39.052"E (S84))

--

-- see also 'intersections'

intersection :: (Spherical a) => MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
intersection :: MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
intersection a1 :: MinorArc a
a1@(MinorArc V3
n1 HorizontalPosition a
s1 HorizontalPosition a
e1) a2 :: MinorArc a
a2@(MinorArc V3
n2 HorizontalPosition a
s2 HorizontalPosition a
e2) =
    case V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
forall a.
Spherical a =>
V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections' V3
n1 V3
n2 (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
s1) of
        Maybe (HorizontalPosition a, HorizontalPosition a)
Nothing -> Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
        (Just (HorizontalPosition a
i1, HorizontalPosition a
i2))
            | HorizontalPosition a -> MinorArc a -> Bool
forall a. Spherical a => HorizontalPosition a -> MinorArc a -> Bool
onMinorArc HorizontalPosition a
pot MinorArc a
a1 Bool -> Bool -> Bool
&& HorizontalPosition a -> MinorArc a -> Bool
forall a. Spherical a => HorizontalPosition a -> MinorArc a -> Bool
onMinorArc HorizontalPosition a
pot MinorArc a
a2 -> HorizontalPosition a -> Maybe (HorizontalPosition a)
forall a. a -> Maybe a
Just HorizontalPosition a
pot
            | Bool
otherwise -> Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
            where mid :: V3
mid =
                      [V3] -> V3
meanV
                          [ HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
s1
                          , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
e1
                          , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
s2
                          , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
e2
                          ]
                  pot :: HorizontalPosition a
pot =
                      if V3 -> V3 -> Double
Math3d.dot (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
i1) V3
mid Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                          then HorizontalPosition a
i1
                          else HorizontalPosition a
i2

-- | Computes the intersections between the two given 'GreatCircle's.

-- Two great circles intersect exactly twice unless there are equal (regardless of orientation), in which case 'Nothing'

-- is returned. For example:

--

-- >>> let gc1 = GreatCircle.headingOn (Geodetic.s84Pos 51.885 0.235) (Angle.decimalDegrees 108.63)

-- >>> let gc2 = GreatCircle.headingOn (Geodetic.s84Pos 49.008 2.549) (Angle.decimalDegrees 32.72)

-- >>> GreatCircle.intersections gc1 gc2

-- Just (50°54'6.201"N,4°29'39.401"E (S84),50°54'6.201"S,175°30'20.598"W (S84))

-- >>> let is = GreatCircle.intersections gc1 gc2

-- >>> fmap fst is == fmap (Geodetic.antipode . snd) is

-- True

intersections ::
       (Spherical a)
    => GreatCircle a
    -> GreatCircle a
    -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections :: GreatCircle a
-> GreatCircle a
-> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections (GreatCircle V3
n1 a
m String
_) (GreatCircle V3
n2 a
_ String
_) = V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
forall a.
Spherical a =>
V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections' V3
n1 V3
n2 a
m

-- | @mean ps@ computes the geographic mean horizontal position of @ps@, if it is defined. For example:

--

-- >>> let ps =

--             [ Geodetic.s84Pos 90 0

--             , Geodetic.s84Pos 60 10

--             , Geodetic.s84Pos 50 (-20)

--             ]

-- >>> GreatCircle.mean ps

-- Just 67°14'10.150"N,6°55'3.040"W (S84)

--

-- The geographic mean is not defined for antipodals positions (since they

-- cancel each other).

--

-- Special conditions:

--

-- @

-- mean [] = Nothing

-- mean [p] = Just p

-- mean [p1, .., antipode p1] = Nothing

-- @

mean :: (Spherical a) => [HorizontalPosition a] -> Maybe (HorizontalPosition a)
mean :: [HorizontalPosition a] -> Maybe (HorizontalPosition a)
mean [] = Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
mean [HorizontalPosition a
p] = HorizontalPosition a -> Maybe (HorizontalPosition a)
forall a. a -> Maybe a
Just HorizontalPosition a
p
mean [HorizontalPosition a]
ps =
    if (HorizontalPosition a -> Bool) -> [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [HorizontalPosition a]
ps) [HorizontalPosition a]
antipodes
        then Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
        else HorizontalPosition a -> Maybe (HorizontalPosition a)
forall a. a -> Maybe a
Just
                 (V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos'
                      ([V3] -> V3
meanV ((HorizontalPosition a -> V3) -> [HorizontalPosition a] -> [V3]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector [HorizontalPosition a]
ps))
                      (HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model (HorizontalPosition a -> a)
-> ([HorizontalPosition a] -> HorizontalPosition a)
-> [HorizontalPosition a]
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head ([HorizontalPosition a] -> a) -> [HorizontalPosition a] -> a
forall a b. (a -> b) -> a -> b
$ [HorizontalPosition a]
ps))
  where
    antipodes :: [HorizontalPosition a]
antipodes = (HorizontalPosition a -> HorizontalPosition a)
-> [HorizontalPosition a] -> [HorizontalPosition a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap HorizontalPosition a -> HorizontalPosition a
forall a. Model a => HorizontalPosition a -> HorizontalPosition a
Geodetic.antipode [HorizontalPosition a]
ps

-- | @projection p ma@ computes the projection of the position @p@ on the minor arc @ma@. Returns 'Nothing' if the

-- position @p@ is the normal of the minor arc or if the projection is not within the minor arc @ma@ (including start

-- & end). For example:

--

-- >>> let p = Geodetic.s84Pos 53.2611 (-0.7972)

-- >>> let ma = fromJust (GreatCircle.minorArc (Geodetic.s84Pos 53.3206 (-1.7297)) (Geodetic.s84Pos 53.1887 0.1334))

-- >>> GreatCircle.projection p ma

-- Just Geodetic.s84Pos 53.25835330666666 (-0.7977433863888889)

projection :: (Spherical a) => HorizontalPosition a -> MinorArc a -> Maybe (HorizontalPosition a)
projection :: HorizontalPosition a -> MinorArc a -> Maybe (HorizontalPosition a)
projection HorizontalPosition a
p ma :: MinorArc a
ma@(MinorArc V3
na HorizontalPosition a
mas HorizontalPosition a
mae) =
    case Maybe V3
mnb of
        Maybe V3
Nothing -> Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
        (Just V3
nb) ->
            case Maybe (HorizontalPosition a, HorizontalPosition a)
is of
                Maybe (HorizontalPosition a, HorizontalPosition a)
Nothing -> Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
                (Just (HorizontalPosition a
i1, HorizontalPosition a
i2)) ->
                    if HorizontalPosition a -> MinorArc a -> Bool
forall a. Spherical a => HorizontalPosition a -> MinorArc a -> Bool
onMinorArc HorizontalPosition a
pot MinorArc a
ma
                        then HorizontalPosition a -> Maybe (HorizontalPosition a)
forall a. a -> Maybe a
Just HorizontalPosition a
pot
                        else Maybe (HorizontalPosition a)
forall a. Maybe a
Nothing
                    where mid :: V3
mid =
                              [V3] -> V3
meanV
                                  [ HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
mas
                                  , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
mae
                                  , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
nap
                                  , HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
nbp
                                  ]
                          pot :: HorizontalPosition a
pot =
                              if V3 -> V3 -> Double
Math3d.dot V3
mid (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
i1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                                  then HorizontalPosition a
i1
                                  else HorizontalPosition a
i2
            where nbp :: HorizontalPosition a
nbp = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' V3
nb a
m -- ensure resolution of lat, lon

                  is :: Maybe (HorizontalPosition a, HorizontalPosition a)
is = V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
forall a.
Spherical a =>
V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections' (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
nap) (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
nbp) a
m
  where
    m :: a
m = HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p
    nap :: HorizontalPosition a
nap = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' V3
na a
m -- ensure resolution of lat, lon

    mnb :: Maybe V3
mnb = HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
nap HorizontalPosition a
p -- normal to great circle (na, p) - if na is p or antipode of p, then projection is not possible.


-- | @side p0 p1 p2@ determines whether @p0@ is left of, right of or on the great circle passing through @p1@ and @p2@ (in this

-- direction). For example:

--

-- >>> GreatCircle.side (Geodetic.s84Pos 10 10) (Geodetic.s84Pos 0 0) (Geodetic.northPole S84)

-- RightOf

-- >>> GreatCircle.side (Geodetic.s84Pos 10 (-10)) (Geodetic.s84Pos 0 0) (Geodetic.northPole S84)

-- LeftOf

-- >>> GreatCircle.side (Geodetic.s84Pos 10 0) (Geodetic.s84Pos 0 0) (Geodetic.northPole S84)

-- None

-- Returns 'None' if @p1@ & @p2@ do not define a great circle (see 'through') or if any of the three position are equal.

side ::
       (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> HorizontalPosition a -> Side
side :: HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Side
side HorizontalPosition a
p0 HorizontalPosition a
p1 HorizontalPosition a
p2
    | HorizontalPosition a
p0 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p1 Bool -> Bool -> Bool
|| HorizontalPosition a
p0 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Side
None
    | Bool
otherwise = Side -> (Double -> Side) -> Maybe Double -> Side
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Side
None Double -> Side
toSide Maybe Double
ms
  where
    ms :: Maybe Double
ms = (V3 -> Double) -> Maybe V3 -> Maybe Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (V3 -> V3 -> Double
Math3d.dot (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p0)) (HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
p1 HorizontalPosition a
p2)

-- | @turn a b c@ computes the angle turned from AB to BC; the angle is positive for left turn, negative for right turn

-- and 0 if all 3 positions are aligned or if any 2 positions are equal. For example:

--

-- >>> GreatCircle.turn (Geodetic.s84Pos 0 0) (Geodetic.s84Pos 45 0) (Geodetic.s84Pos 60 (-10))

-- 18°11'33.741"

turn ::
       (Spherical a)
    => HorizontalPosition a
    -> HorizontalPosition a
    -> HorizontalPosition a
    -> Angle
turn :: HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Angle
turn HorizontalPosition a
a HorizontalPosition a
b HorizontalPosition a
c =
    case (Maybe V3, Maybe V3)
ns of
        (Just V3
n1, Just V3
n2) -> V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween V3
n1 V3
n2 (V3 -> Maybe V3
forall a. a -> Maybe a
Just (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
b))
        (Maybe V3
_, Maybe V3
_) -> Angle
Angle.zero
  where
    ns :: (Maybe V3, Maybe V3)
ns = ((V3 -> V3) -> Maybe V3 -> Maybe V3
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap V3 -> V3
Math3d.unit (HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
a HorizontalPosition a
b), (V3 -> V3) -> Maybe V3 -> Maybe V3
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap V3 -> V3
Math3d.unit (HorizontalPosition a -> HorizontalPosition a -> Maybe V3
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
b HorizontalPosition a
c))

-- private

alongTrackDistance'' ::
       (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
alongTrackDistance'' :: HorizontalPosition a -> HorizontalPosition a -> V3 -> Length
alongTrackDistance'' HorizontalPosition a
p HorizontalPosition a
s V3
n = Angle -> Length -> Length
Angle.arcLength Angle
a (HorizontalPosition a -> Length
forall a. Spherical a => HorizontalPosition a -> Length
radius HorizontalPosition a
s)
  where
    a :: Angle
a =
        V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween
            (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
s)
            (V3 -> V3 -> V3
Math3d.cross (V3 -> V3 -> V3
Math3d.cross V3
n (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p)) V3
n)
            (V3 -> Maybe V3
forall a. a -> Maybe a
Just V3
n)

-- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)]

egdes :: [V3] -> [(V3, V3)]
egdes :: [V3] -> [(V3, V3)]
egdes [V3]
ps = [V3] -> [V3] -> [(V3, V3)]
forall a b. [a] -> [b] -> [(a, b)]
zip [V3]
ps ([V3] -> [V3]
forall a. [a] -> [a]
tail [V3]
ps [V3] -> [V3] -> [V3]
forall a. [a] -> [a] -> [a]
++ [[V3] -> V3
forall a. [a] -> a
head [V3]
ps])

intersections' ::
       (Spherical a) => V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections' :: V3 -> V3 -> a -> Maybe (HorizontalPosition a, HorizontalPosition a)
intersections' V3
n1 V3
n2 a
s
    | (V3 -> Double
Math3d.norm V3
i :: Double) Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 = Maybe (HorizontalPosition a, HorizontalPosition a)
forall a. Maybe a
Nothing
    | Bool
otherwise
    , let ni :: HorizontalPosition a
ni = V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
Geodetic.nvectorPos' (V3 -> V3
Math3d.unit V3
i) a
s = (HorizontalPosition a, HorizontalPosition a)
-> Maybe (HorizontalPosition a, HorizontalPosition a)
forall a. a -> Maybe a
Just (HorizontalPosition a
ni, HorizontalPosition a -> HorizontalPosition a
forall a. Model a => HorizontalPosition a -> HorizontalPosition a
Geodetic.antipode HorizontalPosition a
ni)
  where
    i :: V3
i = V3 -> V3 -> V3
Math3d.cross V3
n1 V3
n2

initialBearing' :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Angle
initialBearing' :: HorizontalPosition a -> HorizontalPosition a -> Angle
initialBearing' HorizontalPosition a
p1 HorizontalPosition a
p2 = Angle -> Angle -> Angle
Angle.normalise Angle
a (Double -> Angle
Angle.decimalDegrees Double
360)
  where
    m :: a
m = HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p1
    v1 :: V3
v1 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p1
    v2 :: V3
v2 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p2
    gc1 :: V3
gc1 = V3 -> V3 -> V3
Math3d.cross V3
v1 V3
v2 -- great circle through p1 & p2

    gc2 :: V3
gc2 = V3 -> V3 -> V3
Math3d.cross V3
v1 (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector (HorizontalPosition a -> V3)
-> (a -> HorizontalPosition a) -> a -> V3
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> HorizontalPosition a
forall a. Model a => a -> HorizontalPosition a
Geodetic.northPole (a -> V3) -> a -> V3
forall a b. (a -> b) -> a -> b
$ a
m) -- great circle through p1 & north pole

    a :: Angle
a = V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween V3
gc1 V3
gc2 (V3 -> Maybe V3
forall a. a -> Maybe a
Just V3
v1)

arcNormal :: (Spherical a) => HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal :: HorizontalPosition a -> HorizontalPosition a -> Maybe V3
arcNormal HorizontalPosition a
p1 HorizontalPosition a
p2
    | HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Maybe V3
forall a. Maybe a
Nothing
    | HorizontalPosition a -> HorizontalPosition a
forall a. Model a => HorizontalPosition a -> HorizontalPosition a
Geodetic.antipode HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Maybe V3
forall a. Maybe a
Nothing
    | Bool
otherwise = V3 -> Maybe V3
forall a. a -> Maybe a
Just (V3 -> V3 -> V3
Math3d.cross (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p1) (HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p2))

-- | @onMinorArc p a@ determines whether position @p@ is on the minor arc @a@.

-- Warning: this function assumes that @p@ is on great circle of the minor arc.

-- return true if chord lengths between (p, start) & (p, end) are both less than

-- chord length between (start, end)

onMinorArc :: (Spherical a) => HorizontalPosition a -> MinorArc a -> Bool
onMinorArc :: HorizontalPosition a -> MinorArc a -> Bool
onMinorArc HorizontalPosition a
p (MinorArc V3
_ HorizontalPosition a
s HorizontalPosition a
e) =
    V3 -> V3 -> Double
Math3d.squaredDistance V3
v0 V3
v1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
l Bool -> Bool -> Bool
&& V3 -> V3 -> Double
Math3d.squaredDistance V3
v0 V3
v2 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
l
  where
    v0 :: V3
v0 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
p
    v1 :: V3
v1 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
s
    v2 :: V3
v2 = HorizontalPosition a -> V3
forall a. HasCoordinates a => a -> V3
Geodetic.nvector HorizontalPosition a
e
    l :: Double
l = V3 -> V3 -> Double
Math3d.squaredDistance V3
v1 V3
v2

-- | reference sphere radius.

radius :: (Spherical a) => HorizontalPosition a -> Length
radius :: HorizontalPosition a -> Length
radius = Ellipsoid -> Length
equatorialRadius (Ellipsoid -> Length)
-> (HorizontalPosition a -> Ellipsoid)
-> HorizontalPosition a
-> Length
forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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

-- | angle between 2 vectors.

angleBetween :: V3 -> V3 -> Angle
angleBetween :: V3 -> V3 -> Angle
angleBetween V3
v1 V3
v2 = V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween V3
v1 V3
v2 Maybe V3
forall a. Maybe a
Nothing

-- | Signed angle between 2 vectors.

-- If @n@ is 'Nothing', the angle is always in [0..pi], otherwise it is in [-pi, +pi],

-- signed + if @v1@ is clockwise looking along @n@, - in opposite direction.

signedAngleBetween :: V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween :: V3 -> V3 -> Maybe V3 -> Angle
signedAngleBetween V3
v1 V3
v2 Maybe V3
n = Double -> Double -> Angle
Angle.atan2 Double
sinO Double
cosO
  where
    sign :: Double
sign = Double -> (V3 -> Double) -> Maybe V3 -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Double
1 (Double -> Double
forall a. Num a => a -> a
signum (Double -> Double) -> (V3 -> Double) -> V3 -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. V3 -> V3 -> Double
Math3d.dot (V3 -> V3 -> V3
Math3d.cross V3
v1 V3
v2)) Maybe V3
n
    sinO :: Double
sinO = Double
sign Double -> Double -> Double
forall a. Num a => a -> a -> a
* V3 -> Double
Math3d.norm (V3 -> V3 -> V3
Math3d.cross V3
v1 V3
v2)
    cosO :: Double
cosO = V3 -> V3 -> Double
Math3d.dot V3
v1 V3
v2

meanV :: [Math3d.V3] -> V3
meanV :: [V3] -> V3
meanV [V3]
vs = V3 -> V3
Math3d.unit (V3 -> V3) -> V3 -> V3
forall a b. (a -> b) -> a -> b
$ (V3 -> V3 -> V3) -> V3 -> [V3] -> V3
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl V3 -> V3 -> V3
Math3d.add V3
Math3d.zero [V3]
vs

toSide :: Double -> Side
toSide :: Double -> Side
toSide Double
s
    | Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = Side
RightOf
    | Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 = Side
LeftOf
    | Bool
otherwise = Side
None