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