module GIS.Math.Spherical where
import Control.Lens
import Control.Lens.Tuple
import Data.Function
import Data.List
import GIS.Types
import Data.Composition
import GIS.Math.Projections
import GIS.Math.Utils
shittyCentroid :: Polygon -> Point
shittyCentroid poly = (avg $ map fst poly, avg $ map snd poly)
avg :: (RealFrac a, Foldable t) => t a -> a
avg list = sum list / (fromIntegral . length $ list)
areaTriangle :: Point -> Point -> Point -> Double
areaTriangle x1 x2 x3 = r^2 * 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
relativeCompactness :: Polygon -> Double
relativeCompactness = (*scale) . compactness1
where scale = (1/4*pi)
compactness1 :: Polygon -> Double
compactness1 p = (areaPolygon p)/(perimeterPolygon p)^2
areaConvex :: Polygon -> Double
areaConvex (base1:base2:pts) = (view _1) $ foldr stepArea (0,base2) pts
where stepArea point (sum, base) = (sum + (areaTriangle base1 base point), point)
areaPolygon :: Polygon -> Double
areaPolygon = (*factor) . areaPolyRectangular . fmap (bonne . toRadians)
where factor = 1717856/4.219690791828533e-2
totalPerimeter :: [Polygon] -> Double
totalPerimeter lines = sum $ fmap perimeterPolygon lines
perimeterPolygon [x1, x2] = distance x1 x2
perimeterPolygon (x1:x2:points) = perimeterPolygon (x2:points) + distance x1 x2
areaPolyRectangular :: Polygon -> Double
areaPolyRectangular (pt:pts) = abs . (*0.5) . fst $ (foldl' areaPolyCalc (0,pt) pts)
where areaPolyCalc (sum,(x1,y1)) (x2, y2) = (sum + (x1 * y2 x2 * y1),(x2,y2))
((x1, y1),(xn, yn)) = (pt, last pts)
distance :: (Double, Double) -> (Double, Double) -> Double
distance = (*6371) .* (on centralAngle toRadians)
centralAngle :: (Double, Double) -> (Double, Double) -> Double
centralAngle (long1, lat1) (long2, lat2) = centralAngle
where centralAngle = acos $ (sin lat1) * (sin lat2) + (cos lat1) * (cos lat2) * (cos (long1 long2))