-- |

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