{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls #-}
-- |Author: Thomas DuBuisson
-- Copyright: Thomas DuBuisson
-- License: BSD3
--
-- A basic GPS library with calculations for distance and speed along
-- with helper functions for filtering/smoothing trails.  All distances
-- are in meters and time is in seconds.  Speed is thus meters/second
--
-- The current intent of this library is to
--   1) Provide a standard type class interface for coordinates
--   2) Include fleshed out support for relevant libraries, a task
--      too often neglected in modern Hackage packages.
--
--  Integration includes KML support via the xml package, pretty printing,
--  anb binary.

module Data.GPS
	( -- * Types and Classes
	  Distance
	, Heading
	, Speed
	, Vector
	, Trail
	-- * Constants
	, north
	, south
	, east
	, west
	, radiusOfEarth
	-- * Helper Functions
	, restLocations
	, closestDistance
	, getRadianPair
	, addVector
	, divideArea
	, convexHull
	) where

import Data.Ord (comparing)
import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy)
import Data.Geo.GPX.PtType
import Data.Geo.GPX.LongitudeType
import Data.Geo.GPX.LatitudeType
import Data.Geo.GPX.Accessor.Time
import Data.Geo.GPX.Accessor.Lon
import Data.Geo.GPX.Accessor.Lat

import Text.XML.XSD.DateTime(DateTime, toUTCTime)

import Data.Time
import Data.Maybe (listToMaybe)
import Data.Fixed (mod')

-- |Distances are expressed in meters
type Distance = Double

-- |Angles are expressed in radians from North.
-- 	0	== North
-- 	pi/2	== West
-- 	pi 	== South
-- 	(3/2)pi	== East    == - (pi / 2)
type Heading = Double

-- |Speed is hard coded as meters per second
type Speed = Double

type Vector = (Distance, Heading)
type Trail a = [a]

getUTCTime :: (Lat a, Lon a, Time a) => a -> Maybe UTCTime
getUTCTime = fmap toUTCTime . time

distance :: (Lat a, Lon a) => a -> a -> Distance
distance a b =
	radiusOfEarth * acos( sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2 - lon1) )
 where
  (lat1, lon1) = getRadianPairD a
  (lat2, lon2) = getRadianPairD b

heading         :: (Lat a, Lon a) => a -> a -> Heading	-- ^ 0 = North, pi/2 = West...
heading a b =
	atan2	(sin (diffLon) * cos (lat2)) 
		(cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))
 where
  (lat1, lon1) = getRadianPairD a
  (lat2, lon2) = getRadianPairD b
  diffLon = lon1 - lon2

getVector :: (Lat a, Lon a) => a -> a -> Vector
getVector a b = (distance a b, heading a b)

-- | Speed in meters per second
speed :: (Lat loc, Lon loc, Time loc) => loc -> loc -> Maybe Speed
speed a b = 
	case (getUTCTime b, getUTCTime a) of
		(Just x, Just y) -> Just $ realToFrac (diffUTCTime x y) / (distance a b)
		_ -> Nothing

data TempTrail a = T (Trail a) a

-- |radius of the earth in meters
radiusOfEarth :: Double
radiusOfEarth = 6378700

-- |North is 0 radians
north :: Heading
north = 0

-- |South, being 180 degrees from North, is pi.
south :: Heading
south = pi

-- |East is 270 degrees from North
east :: Heading
east = (3 / 2) * pi

-- |West is 90 degrees (pi/2)
west :: Heading
west = pi / 2

toDecimal = (*) (180 / pi)


-- |Given a vector and coordinate, computes a new coordinate.
-- Within some epsilon it should hold that if
--
-- 	@dest = addVector (dist,heading) start@
--
-- then
--
-- 	@heading == dmsHeading start dest@
-- 	
-- 	@dist    == distance start dest@
addVector :: (Lat c, Lon c) => Vector -> c -> c
addVector (d,h) p = setLon (longitudeType lon2) . setLat (latitudeType lat2) $ p
  where
	(lat,lon) = getRadianPairD p
	lat2 = lat + (cos h) * (d / radiusOfEarth)
	lon2 = lon + acos ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2))

getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double)
getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair

-- |Provides a lat/lon pair of doubles in radians
getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)
getRadianPair p = (t, g)
  where
    t = toRadians (lat p)
    g = toRadians (lon p)

toRadians :: Floating f => f -> f
toRadians = (*) (pi / 180)

-- |Filters out any points that go backward in time (thus must not be valid if this is a trail)
linearTime :: (Lat a, Lon a, Time a) => Trail a -> Trail a
linearTime [] = []
linearTime (p:ps) = go (getUTCTime p) ps
  where
  go _ [] = []
  go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps

-- |Creates a list of trails all of which are within the given distance of each
-- other spanning atleast the given amount of time.
--
-- For example @restLocations 50 600@
-- would return lists of all points that are within 50 meters of each other and
-- span at least 10 minutes (600 seconds).
--
-- Note this gives points within fifty meters of the earliest point - wandering
-- in a rest area with a 50 meter radius could result in several rest points
-- ([a,b..]) or even none if the distance between individual points exceeds 50m.
restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a]
restLocations d s xs = go xs
  where
  go [] = []
  go (a:as) =
	let (lst, close, far) = takeWhileLast ((<= d) . distance a) as
	in case lst of
		Just x -> case (getUTCTime a, getUTCTime x) of
				(Just a', Just x') ->
					let d = diffUTCTime a' x'
			  		in if d >= s then close : go far else go as
				_ -> go as
		Nothing -> go as

takeWhileLast :: (a -> Bool) -> [a] -> (Maybe a, [a], [a])
takeWhileLast p [] = (Nothing, [], [])
takeWhileLast p (x:xs)
	| not (p x) = (Nothing, [], x:xs)
	| otherwise = go x xs
  where
  go a [] = (Just a, [a], [])
  go a (b:bs)
	| p b = let (c,d,f) = go b bs in (c, a:d, f)
	| otherwise = (Just a, [a], b:bs)

-- |Returns the closest distance between two trails (or Nothing if a trail is empty)
-- O( (n * m) * log (n * m) )
closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance
closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs]


-- |@divideArea vDist hDist nw se@ divides an area into a grid of equally
-- spaced coordinates within the box drawn by the northwest point (nw) and
-- southeast point (se).  Because this uses floating point there might be a
-- different number of points in some rows (the last might be too far east based
-- on a heading from the se point).
divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]
divideArea vDist hDist nw se =
	let (top,left)  = (lat nw, lon nw)
	    (btm,right) = (lat se, lon se)
	    columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw
	    buildRow  = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))
	in map buildRow columnOne

-- | Uses Grahams scan to compute the convex hull of the given points.
-- This operation requires sorting of the points, so don't try it unless
-- you have notably more memory than the list of points will consume.
convexHull :: (Eq c, Lat c, Lon c) => [c] -> [c]
convexHull xs =
	let first = southMost xs
	in case first of
		Nothing -> []
		Just f  ->
	    	     let sorted = sortBy (comparing ((`mod'` (2*pi)). (+ pi/2). heading f)) (filter (/= f) xs)
		     in case sorted of
			(a:b:cs) -> grahamScan (b:a:f:[]) cs
			cs       -> f : cs
  where
  grahamScan [] _ = []
  grahamScan ps [] = ps
  grahamScan (x:[]) _ = [x]
  grahamScan (p2:p1:ps) (x:xs) =
	case turn p1 p2 x of
		LeftTurn  -> grahamScan (x:p2:p1:ps) xs
		Straight  -> grahamScan (x:p2:p1:ps) xs
		_	  -> grahamScan (p1:ps) (x:xs)

data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)

turn :: (Lat c, Lon c) => c -> c -> c -> Turn
turn a b c =
	let h1 = heading a b
	    h2 = heading b c
	in case compare (h2 - h1) 0 of
			LT -> RightTurn
			GT -> LeftTurn
			EQ -> Straight

-- | Find the southmost point
southMost :: (Lat c) => [c] -> Maybe c
southMost []  = Nothing
southMost cs = Just . minimumBy (comparing lat) $ cs