```{-# LANGUAGE TupleSections #-}
module Geo.Computations.Basic
( -- * Types
Distance
, Speed
, Vector
, Trail
, Circle
, Arc
, Point(..)
, pt
-- * Constants
, north
, south
, east
, west
, circumferenceOfEarth
-- * Coordinate Functions
, distance
, speed
, getVector
, interpolate
, circleIntersectionPoints
, intersectionArcsOf
, maximumDistanceOfArc
-- * IO helpers
) where

import Data.Time
import Data.Maybe
import Geo.GPX.Conduit

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

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

-- | Genearlly a circle indicates a known area in which we are searching
-- (so a center point and maximum possible distance from that point)
type Circle a = (a, Distance)

-- | An arc is represented as a circle, starting heading and ending heading

type Trail a = [a]

distance :: Point -> Point -> Distance
distance x y =
deltaLat    = lat2 - lat1
deltaLon    = lon2 - lon1
a = (sin (deltaLat / 2))^(2::Int) + cos lat1 * cos lat2 * (sin (deltaLon / 2))^(2::Int)
c = 2 * atan2 (a**0.5) ((1-a)**0.5)

-- | Direction two points aim toward (0 = North, pi/2 = West, pi = South, 3pi/2 = East)
atan2	(sin (diffLon) * cos (lat2))
(cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))
where
diffLon = lon2 - lon1

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

-- |Given a vector and coordinate, computes a new coordinate.
-- Within some epsilon it should hold that if
--
--
-- then
--
--
-- 	@dist    == distance start dest@
addVector :: Vector -> Point -> Point
p { pntLon = toDegrees lon2
, pntLat = toDegrees lat2
}
where
lat2 = asin (sin lat * cos (d / radiusOfEarth) + cos lat
lon2 = lon - atan2 (sin h * sin (d / radiusOfEarth) * cos lat)
(cos (d/radiusOfEarth) - sin lat * sin lat2)

-- | Speed in meters per second, only if a 'Time' was recorded for each waypoint.
speed :: Point -> Point -> Maybe Speed
speed a b =
case (pntTime b, pntTime a) of
(Just x, Just y) ->
let timeDiff = realToFrac (diffUTCTime x y)
in if timeDiff == 0 then Nothing else Just \$ (distance a b) / timeDiff
_ -> Nothing

-- |radius of the earth in meters

-- |Circumference of earth (meters)
circumferenceOfEarth :: Double
circumferenceOfEarth = radiusOfEarth * 2 * pi

north = 0

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

-- |East is 270 degrees (3 pi / 2)
east = (3 / 2) * pi

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

toDegrees :: Double -> Double
toDegrees = (*) (180 / pi)

-- get latitude and longituide in Radians as Double's

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

-- | @interpolate c1 c2 w@ where @0 <= w <= 1@ Gives a point on the line
-- between c1 and c2 equal to c1 when @w == 0@ (weighted linearly
-- toward c2).
interpolate :: Point -> Point -> Double -> Point
interpolate c1 c2 w
| w < 0 = c1
| w > 1 = c2
| otherwise =
let (h,d) = (heading c1 c2, distance c1 c2)
v = (d * w, h)

-- | Compute the points at which two circles intersect (assumes a flat plain).  If
-- the circles do not intersect or are identical then the result is @Nothing@.
circleIntersectionPoints :: (Point, Distance) -> (Point, Distance) -> Maybe (Point,Point)
circleIntersectionPoints (a,r1) (b,r2)
| pntLat a == pntLat b && pntLon a == pntLon b && r1 == r2 = Nothing -- FIXME need approx eq
| r1 + r2 < ab = Nothing
| any isNaN (map pntLat pts) || any isNaN (map pntLon pts) = Nothing
| otherwise = Just (p1, p2)
where
ab = distance a b
angABX = acos ( (r1^(2::Int) + ab^(2::Int) - r2^(2::Int)) / (2 * r1 * ab) )
ang1 = heading a b + angABX
ang2 = heading a b - angABX
p1 = addVector (r1, ang1) a
p2 = addVector (r1, ang2) a
pts = [p1,p2]

-- | Find the area in which all given circles intersect.  The resulting
-- area is described in terms of the bounding arcs.   All cirlces must
-- intersect at two points.
intersectionArcsOf :: [Circle Point] -> [Arc Point]
intersectionArcsOf cs =
let isArcWithinCircle circ arc = maximumDistanceOfArc (fst circ) arc <= (snd circ)
isArcWithinAllCircles arc = all (\$ arc) (map isArcWithinCircle cs)
-- getArcs :: Circle a -> Circle a -> [Arc a]
getArcs c1 c2 = concatMap (buildArcsFromPoints c1 c2) . maybeToList \$ circleIntersectionPoints c1 c2
-- buildArcsFromPoints :: (a, a) -> [Arc a]
buildArcsFromPoints c1 c2 (p1,p2) =
let c1h1 = heading (fst c1) p1
c1h2 = heading (fst c1) p2
c2h1 = heading (fst c2) p1
c2h2 = heading (fst c2) p2
in [(c1,c1h1,c1h2), (c1,c1h2, c1h1), (c2,c2h1,c2h2), (c2,c2h2,c2h1)]
in filter isArcWithinAllCircles . concatMap (uncurry getArcs) . choose2 \$ cs

maximumDistanceOfArc :: Point -> Arc Point -> Distance
maximumDistanceOfArc pnt ((c,r), h1, h2) =