GreatCircle.hs

-- |
-- Module:      Data.Geo.Jord.GreatCircle
-- Copyright:   (c) 2018 Cedric Liegeois
-- License:     BSD3
-- Maintainer:  Cedric Liegeois <ofmooseandmen@yahoo.fr>
-- Stability:   experimental
-- Portability: portable
--
-- Types and functions for working with <https://en.wikipedia.org/wiki/Great_circle Great Circle>.
--
-- All functions are implemented using the vector-based approached described in
-- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation>
--
-- This module assumes a spherical earth.
--
module Data.Geo.Jord.GreatCircle
    (
    -- * The 'GreatCircle' type
    GreatCircle
    -- * Smart constructors
    , greatCircle
    , greatCircleE
    , greatCircleF
    , greatCircleBearing
    -- * Geodesic calculations
    , crossTrackDistance
    , crossTrackDistance'
    , intersections
    , isInside
    ) where

import Control.Monad.Fail
import Data.Geo.Jord.Angle
import Data.Geo.Jord.LatLong
import Data.Geo.Jord.Length
import Data.Geo.Jord.NVector
import Data.Geo.Jord.Position
import Data.Geo.Jord.Quantity
import Data.Maybe (fromMaybe)
import Prelude hiding (fail)

-- | A circle on the surface of the Earth which lies in a plane passing through
-- the Earth's centre. Every two distinct and non-antipodal points on the surface
-- of the Earth define a Great Circle.
--
-- It is internally represented as its normal vector - i.e. the normal vector
-- to the plane containing the great circle.
--
-- See 'greatCircle', 'greatCircleE', 'greatCircleF' or 'greatCircleBearing' constructors.
--
data GreatCircle = GreatCircle
    { normal :: NVector
    , dscr :: String
    } deriving (Eq)

instance Show GreatCircle where
    show = dscr

-- | 'GreateCircle' passing by both given 'Position's. 'error's if given positions are
-- equal or antipodal.
greatCircle :: (Eq a, Position a, Show a) => a -> a -> GreatCircle
greatCircle p1 p2 =
    fromMaybe
        (error (show p1 ++ " and " ++ show p2 ++ " do not define a unique Great Circle"))
        (greatCircleF p1 p2)

-- | 'GreateCircle' passing by both given 'Position's. A 'Left' indicates that given positions are
-- equal or antipodal.
greatCircleE :: (Eq a, Position a) => a -> a -> Either String GreatCircle
greatCircleE p1 p2
    | p1 == p2 = Left "Invalid Great Circle: positions are equal"
    | p1 == antipode p2 = Left "Invalid Great Circle: positions are antipodal"
    | otherwise =
        Right
            (GreatCircle
                 (cross v1 v2)
                 ("passing by " ++
                  show (fromNVector v1 :: LatLong) ++ " & " ++ show (fromNVector v2 :: LatLong)))
  where
    v1 = toNVector p1
    v2 = toNVector p2

-- | 'GreateCircle' passing by both given 'Position's. 'fail's if given positions are
-- equal or antipodal.
greatCircleF :: (Eq a, MonadFail m, Position a) => a -> a -> m GreatCircle
greatCircleF p1 p2 =
    case e of
        Left err -> fail err
        Right gc -> return gc
  where
    e = greatCircleE p1 p2

-- | 'GreatCircle' passing by the given 'Position' and heading on given bearing.
greatCircleBearing :: (Position a) => a -> Angle -> GreatCircle
greatCircleBearing p b =
    GreatCircle
        (sub n' e')
        ("passing by " ++ show (fromNVector v :: LatLong) ++ " heading on " ++ show b)
  where
    v = toNVector p
    e = cross northPole v -- easting
    n = cross v e -- northing
    e' = scale e (cos' b / norm e)
    n' = scale n (sin' b / norm n)

-- | 'crossTrackDistance'' assuming a radius of 'meanEarthRadius'.
crossTrackDistance :: (Position a) => a -> GreatCircle -> Length
crossTrackDistance p gc = crossTrackDistance' p gc meanEarthRadius

-- | Signed distance from given 'Position' to given 'GreatCircle'.
-- Returns a negative 'Length' if position if left of great circle,
-- positive 'Length' if position if right of great circle; the orientation of the
-- great circle is therefore important:
--
-- @
--     let gc1 = greatCircle (latLongDecimal 51 0) (latLongDecimal 52 1)
--     let gc2 = greatCircle (latLongDecimal 52 1) (latLongDecimal 51 0)
--     crossTrackDistance p gc1 == (- crossTrackDistance p gc2)
-- @
crossTrackDistance' :: (Position a) => a -> GreatCircle -> Length -> Length
crossTrackDistance' p gc =
    arcLength (sub (angularDistance (normal gc) (toNVector p) Nothing) (decimalDegrees 90))

-- | Computes the intersections between the two given 'GreatCircle's.
-- Two 'GreatCircle's intersect exactly twice unless there are equal (regardless of orientation),
-- in which case 'Nothing' is returned.
intersections :: (Position a) => GreatCircle -> GreatCircle -> Maybe (a, a)
intersections gc1 gc2
    | norm i == 0.0 = Nothing
    | otherwise
    , let ni = unit i = Just (fromNVector ni, fromNVector (antipode ni))
  where
    i = cross (normal gc1) (normal gc2)