```{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}

module Geodetics.Grid (
-- ** Grid types
GridClass (..),
GridPoint (..),
GridOffset (..),
-- ** Grid operations
polarOffset,
offsetScale,
offsetNegate,
applyOffset,
offsetDistance,
offsetDistanceSq,
offsetBearing,
gridOffset,
-- ** Unsafe conversion
unsafeGridCoerce,
-- ** Utility functions for grid references
fromGridDigits,
toGridDigits
) where

import Data.Char
import Data.Function
import Data.Monoid (Monoid)
import Data.Semigroup (Semigroup, (<>))
import Geodetics.Altitude
import Geodetics.Geodetic
import Numeric.Units.Dimensional.Prelude hiding ((.))
import qualified Prelude as P

-- | A Grid is a two-dimensional projection of the ellipsoid onto a plane. Any given type of grid can
-- usually be instantiated with parameters such as a tangential point or line, and these parameters
-- will include the terrestrial reference frame ("Ellipsoid" in this library) used as a foundation.
-- Hence conversion from a geodetic to a grid point requires the \"basis\" for the grid in question,
-- and grid points carry that basis with them because without it there is no defined relationship
-- between the grid points and terrestrial positions.
class GridClass r e | r->e where
fromGrid :: GridPoint r -> Geodetic e
toGrid :: r -> Geodetic e -> GridPoint r
gridEllipsoid :: r -> e

-- | A point on the specified grid.
data GridPoint r = GridPoint {
eastings, northings, altGP :: Length Double,
gridBasis :: r
} deriving (Show)

instance Eq (GridPoint r) where
p1 == p2  =
eastings p1 == eastings p2 &&
northings p1 == northings p2 &&
altGP p1 == altGP p2

instance HasAltitude (GridPoint g) where
altitude = altGP
setAltitude h gp = gp{altGP = h}

-- | A vector relative to a point on a grid.
-- Operations that use offsets will only give
-- meaningful results if all the points come from the same grid.
--
-- The monoid instance is the sum of offsets.
data GridOffset = GridOffset {
deltaEast, deltaNorth, deltaAltitude :: Length Double
} deriving (Eq, Show)

instance Semigroup GridOffset where
g1 <> g2 = GridOffset (deltaEast g1 + deltaEast g2)
(deltaNorth g1 + deltaNorth g2)
(deltaAltitude g1 + deltaAltitude g2)

instance Monoid GridOffset where
mempty = GridOffset _0 _0 _0
mappend = (<>)

-- | An offset defined by a distance and a bearing to the right of North.
--
-- There is no elevation parameter because we are using a plane to approximate an ellipsoid,
-- so elevation would not provide a useful result.  If you want to work with elevations
-- then "rayPath" will give meaningful results.
polarOffset :: Length Double -> Angle Double -> GridOffset
polarOffset r d = GridOffset (r * sin d) (r * cos d) _0

-- | Scale an offset by a scalar.
offsetScale :: Dimensionless Double -> GridOffset -> GridOffset
offsetScale s off = GridOffset (deltaEast off * s)
(deltaNorth off * s)
(deltaAltitude off * s)

-- | Invert an offset.
offsetNegate :: GridOffset -> GridOffset
offsetNegate off = GridOffset (negate \$ deltaEast off)
(negate \$ deltaNorth off)
(negate \$ deltaAltitude off)

-- Add an offset on to a point to get another point.
applyOffset :: GridOffset -> GridPoint g -> GridPoint g
applyOffset off p = GridPoint (eastings p + deltaEast off)
(northings p + deltaNorth off)
(altitude p + deltaAltitude off)
(gridBasis p)

-- | The distance represented by an offset.
offsetDistance :: GridOffset -> Length Double
offsetDistance = sqrt . offsetDistanceSq

-- | The square of the distance represented by an offset.
offsetDistanceSq :: GridOffset -> Area Double
offsetDistanceSq off =
deltaEast off ^ pos2 + deltaNorth off ^ pos2 + deltaAltitude off ^ pos2

-- | The direction represented by an offset, as bearing to the right of North.
offsetBearing :: GridOffset -> Angle Double
offsetBearing off = atan2 (deltaEast off) (deltaNorth off)

-- | The offset required to move from p1 to p2.
gridOffset :: GridPoint g -> GridPoint g -> GridOffset
gridOffset p1 p2 = GridOffset (eastings p2 - eastings p1)
(northings p2 - northings p1)
(altitude p2 - altitude p1)

-- | Coerce a grid point of one type into a grid point of a different type,
-- but with the same easting, northing and altitude. This is unsafe because it
-- will produce a different position unless the two grids are actually equal.
--
-- It should be used only to convert between distinguished grids (e.g. "UkNationalGrid") and
-- their equivalent numerical definitions.
unsafeGridCoerce :: b -> GridPoint a -> GridPoint b
unsafeGridCoerce base p = GridPoint (eastings p) (northings p) (altitude p) base

-- | Convert a list of digits to a distance. The first argument is the size of the
-- grid square within which these digits specify a position. The first digit is
-- in units of one tenth of the grid square, the second one hundredth, and so on.
-- The first result is the lower limit of the result, and the second is the size
-- of the specified offset.
-- So for instance @fromGridDigits (100 *~ kilo meter) "237"@ will return
--
-- > Just (23700 meters, 100 meters)
--
-- If there are any non-digits in the string then the function returns @Nothing@.
fromGridDigits :: Length Double -> String -> Maybe (Length Double, Length Double)
fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing
where
n = length ds
d = sum \$ zipWith (*)
(map ((*~ one) . fromIntegral . digitToInt) ds)
(tail \$ iterate (/ (10 *~ one)) sq)
p = sq / ((10 *~ one) ** (fromIntegral n *~ one))

-- | Convert a distance into a digit string suitable for printing as part
-- of a grid reference. The result is the nearest position to the specified
-- number of digits, expressed as an integer count of squares and a string of digits.
-- If any arguments are invalid then @Nothing@ is returned.
toGridDigits ::
Length Double    -- ^ Size of enclosing grid square. Must be at least 1 km.
-> Int           -- ^ Number of digits to return. Must be positive.
-> Length Double -- ^ Offset to convert into grid.
-> Maybe (Integer, String)
toGridDigits sq n d =
if sq < (1 *~ kilo meter) || n < 0 || d < _0
then Nothing
else