-- |

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

--

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

    (

    -- * The 'Position' type

      Position(..)

    -- * Geodetic calculations

    , angularDistance

    , antipode

    , destination

    , destination'

    , distance

    , distance'

    , finalBearing

    , initialBearing

    , interpolate

    , isInside

    , mean

    -- * Misc.

    , meanEarthRadius

    , northPole

    , southPole

    ) where



import Data.Geo.Jord.Angle

import Data.Geo.Jord.LatLong

import Data.Geo.Jord.Length

import Data.Geo.Jord.NVector

import Data.Geo.Jord.Quantity

import Data.List (subsequences)

import Prelude hiding (fail)



-- | The 'Position' class defines 2 functions to convert a position to and from a 'NVector'.

-- All functions in this module first convert 'Position' to 'NVector' and any resulting 'NVector' back

-- to a 'Position'. This allows the call site to pass either 'NVector' or 'LatLong' and to get back

-- the same class instance.

class Position a where

    -- | Converts a 'NVector' into 'Position' instance.

    fromNVector :: NVector -> a

    -- | Converts the 'Position' instance into a 'NVector'.

    toNVector :: a -> NVector



-- | 'LatLong' to/from 'NVector'.

instance Position LatLong where

    fromNVector v = latLong lat lon

      where

        lat = atan2' (z v) (sqrt (x v * x v + y v * y v))

        lon = atan2' (y v) (x v)

    toNVector g = nvector x' y' z'

      where

        lat = latitude g

        lon = longitude g

        cl = cos' lat

        x' = cl * cos' lon

        y' = cl * sin' lon

        z' = sin' lat



-- | Identity.

instance Position NVector where

    fromNVector v = v

    toNVector v = v



-- | Angle between the two given 'NVector's.

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

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

angularDistance :: NVector -> NVector -> Maybe NVector -> Angle

angularDistance v1 v2 n = atan2' sinO cosO

  where

    sign = maybe 1 (signum . dot (cross v1 v2)) n

    sinO = sign * norm (cross v1 v2)

    cosO = dot v1 v2



-- | Returns the antipodal 'Position' of the given 'Position' - i.e. the position on the surface

-- of the Earth which is diametrically opposite to the given position.

antipode :: (Position a) => a -> a

antipode p = fromNVector (scale (toNVector p) (-1.0))



-- | 'destination'' assuming a radius of 'meanEarthRadius'.

destination :: (Position a) => a -> Angle -> Length -> a

destination p b d = destination' p b d meanEarthRadius



-- | Computes the destination 'Position' from the given 'Position' having travelled the given distance on the

-- given initial bearing (bearing will normally vary before destination is reached) and using the given earth radius.

--

-- This is known as the direct geodetic problem.

destination' :: (Position a) => a -> Angle -> Length -> Length -> a

destination' p b d r

    | isZero d = p

    | otherwise = fromNVector (add (scale v (cos' ta)) (scale de (sin' ta)))

  where

    v = toNVector p

    ed = unit (cross northPole v) -- east direction vector at v

    nd = cross v ed -- north direction vector at v

    ta = central d r -- central angle

    de = add (scale nd (cos' b)) (scale ed (sin' b)) -- unit vector in the direction of the azimuth



-- | 'distance'' assuming a radius of 'meanEarthRadius'.

distance :: (Position a) => a -> a -> Length

distance p1 p2 = distance' p1 p2 meanEarthRadius



-- | Computes the surface distance (length of geodesic) in 'Meters' assuming a

-- spherical Earth between the two given 'Position's and using the given earth radius.

distance' :: (Position a) => a -> a -> Length -> Length

distance' p1 p2 = arcLength (angularDistance v1 v2 Nothing)

  where

    v1 = toNVector p1

    v2 = toNVector p2



-- | Computes the final bearing arriving at given destination  @p2@ 'Position' from given 'Position' @p1@.

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

-- Returns 180 if both position are equals.

finalBearing :: (Position a) => a -> a -> Angle

finalBearing p1 p2 = normalise (initialBearing p2 p1) (decimalDegrees 180)



-- | Computes the initial bearing from given @p1@ 'Position' to given @p2@ 'Position', in compass degrees.

-- Returns 0 if both position are equals.

initialBearing :: (Position a) => a -> a -> Angle

initialBearing p1 p2 = normalise (angularDistance gc1 gc2 (Just v1)) (decimalDegrees 360)

  where

    v1 = toNVector p1

    v2 = toNVector p2

    gc1 = cross v1 v2 -- great circle through p1 & p2

    gc2 = cross v1 northPole -- great circle through p1 & north pole



-- | Computes the 'Position' at given fraction @f@ between the two given 'Position's @p0@ and @p1@.

--

-- Special conditions:

--

-- @

--     interpolate p0 p1 0.0 => p0

--     interpolate p0 p1 1.0 => p1

-- @

--

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

--

interpolate :: (Position a) => a -> a -> Double -> a

interpolate p0 p1 f

    | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f)

    | f == 0 = p0

    | f == 1 = p1

    | otherwise = fromNVector (unit (add v0 (scale (sub v1 v0) f)))

  where

    v0 = toNVector p0

    v1 = toNVector p1



-- | Determines whether the given 'Position' is inside the polygon defined by the given list of 'Position's.

-- The polygon is closed if needed (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 positions does not at least defines a triangle.

--

isInside :: (Eq a, Position a) => a -> [a] -> Bool

isInside p ps

    | null ps = False

    | head ps == last ps = isInside p (init ps)

    | length ps < 3 = False

    | otherwise =

        let aSum = foldl (\a v' -> add a (uncurry angularDistance v' (Just v))) (decimalDegrees 0) es

         in abs (toDecimalDegrees aSum) > 180.0

  where

    v = toNVector p

    es = egdes (map (sub v . toNVector) ps)



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

egdes :: [NVector] -> [(NVector, NVector)]

egdes ps = zip ps ps'

  where

    ps' = tail ps ++ [head ps]



-- | Computes the geographic mean 'Position' of the given 'Position's if it is defined.

--

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

-- cancel each other).

--

-- Special conditions:

--

-- @

--     mean [] == Nothing

--     mean [p] == Just p

--     mean [p1, p2, p3] == Just circumcentre

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

-- @

--

mean :: (Position a) => [a] -> Maybe a

mean [] = Nothing

mean [p] = Just p

mean ps =

    if null antipodals

        then Just (fromNVector (unit (foldl add zero vs)))

        else Nothing

  where

    vs = map toNVector ps

    ts = filter (\l -> length l == 2) (subsequences vs)

    antipodals =

        filter

            (\t -> (fromNVector (antipode (head t)) :: LatLong) == (fromNVector (last t) :: LatLong))

            ts



-- | Mean Earth radius: 6,371,008.8 metres.

meanEarthRadius :: Length

meanEarthRadius = metres 6371008.8



-- | 'Position' of the North Pole.

northPole :: (Position a) => a

northPole = fromNVector (nvector 0.0 0.0 1.0)



-- | 'Position' of the South Pole.

southPole :: (Position a) => a

southPole = fromNVector (nvector 0.0 0.0 (-1.0))