-- | Utilities to compute area, perimeterPolygon, etc. on the surface of a sphere.
module GIS.Math.Spherical ( areaPolygon
                          , perimeterPolygon
                          , relativeCompactness
                          , totalPerimeter
                          -- * Spherical geometry
                          , areaTriangle
                          , areaConvex
                          -- * Internal functions
                          , distance
                          ) where

import           Control.Composition
import           GIS.Math.Projections
import           GIS.Math.Utils
import           GIS.Types

-- | Compute the area of a triangle using L'Huillier's formula
areaTriangle :: Point -> Point -> Point -> Double
areaTriangle x1 x2 x3 = r^(2 :: Int) * e
    where r = 6371
          e = 4 * atan(sqrt(tan(s/2) * tan((s - a)/2) * tan((s - b)/2) * tan((s - c)/2)))
          s = (a + b + c) / 2
          a = distanceRad x1 x2
          b = distanceRad x1 x3
          c = distanceRad x2 x3
          distanceRad = on centralAngle toRadians

-- TODO mandelbrot/fractal dimension?
-- consider "area of largest circumscribable circle" as well.

-- | Relative compactness. Dimensionless.
relativeCompactness :: Polygon -> Double
relativeCompactness = (*scale) . compactness1
    where scale = 1/4*pi

-- | Take the area of the polygon and divide by the perimeter squared. Dimensionless.
compactness1 :: Polygon -> Double
compactness1 p = areaPolygon p/perimeterPolygon p^(2 :: Int)

-- | Compute the area of a convex polygon on the surface of a sphere.
areaConvex :: Polygon -> Double
areaConvex (base1:base2:pts) = fst $ foldr stepArea (0,base2) pts
    where stepArea point (sum', base) = (sum' + areaTriangle base1 base point, point)
areaConvex _ = error "attempted to take area of polygon with < 3 points"

-- | Uses areal projection; then finds area of the polygon.
-- Result is in km^2
areaPolygon :: Polygon -> Double
areaPolygon = (*factor) . areaPolyRectangular . fmap (bonne . toRadians)
    where factor = 1717856/4.219690791828533e-2

-- | Given a list of polygons, return the total perimeter.
totalPerimeter :: [Polygon] -> Double
totalPerimeter = sum . fmap perimeterPolygon

perimeterPolygon :: Polygon -> Double
perimeterPolygon [x1, x2]       = distance x1 x2
perimeterPolygon (x1:x2:points) = perimeterPolygon (x2:points) + distance x1 x2
perimeterPolygon _              = error "Attempted to take area of polygon with no points"

-- | Find the area of a polygon with rectangular coördinates given. This
-- function will error if given a polygon with no points.
areaPolyRectangular :: Polygon -> Double
areaPolyRectangular (pt:pts) = abs . (*0.5) . fst $ foldr areaPolyCalc (0,pt) pts
    where areaPolyCalc (x2, y2) (sum',(x1,y1)) = (sum' + (x1 * y2 - x2 * y1),(x2,y2))
areaPolyRectangular _ = error "Attempted to take area of polygon with no points"

-- | Distance in kilometers between two points given in degrees.
distance :: Point -> Point -> Double
distance = (*6371) .* on centralAngle toRadians

-- | Compute central angle from points given in radians
centralAngle :: Point -> Point -> Double
centralAngle (long1, lat1) (long2, lat2) =
    acos $ sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (long1 - long2)