{-# LANGUAGE RankNTypes #-}
module Math.Geometry.Spherical
(
areaTriangle
, compactness
, relativeCompactness
, areaConvex
, perimeterPolygon
, areaPolygon
, distance
, albers
, littow
, craig
, winkel3
, mercator
, bonne
, washingtonDC
, mecca
, radians
, toRadians
, project
) where
import Control.Composition
radians :: Floating a => a -> a
radians = (*(pi/180))
toRadians :: Floating a => (a, a) -> (a, a)
toRadians = both radians
sinc :: (Eq a, Floating a) => a -> a
sinc 0 = 1
sinc x = sin x / x
washingtonDC :: (Floating a) => (a, a)
washingtonDC = toRadians (38.9072, -77.0369)
mecca :: (Floating a) => (a, a)
mecca = toRadians (21.3891, 39.8579)
littow :: Floating a => (a, a) -> (a, a)
littow (long, lat) = (sin (long - referenceLong) /cos lat, cos (long - referenceLong) * tan lat)
where referenceLong = radians 0
craig :: (Floating a, Eq a)
=> (a, a)
-> ((a, a) -> (a, a))
craig referencePoint (long, lat) = (long - referenceLong, y)
where (referenceLong, referenceLat) = referencePoint
expr = sin lat * cos (long - referenceLong) - tan referenceLat * cos lat
y | long - referenceLong == 0 = expr
| otherwise = (long - referenceLong) / sin (long - referenceLong) * expr
winkel3 :: (Eq a, Floating a) => (a, a) -> (a, a)
winkel3 (long, lat) = ((lambda * cos phi1 + (2 * cos lat * sin (lambda/2)/sinc alpha))/2, (lat + sin lat/sinc alpha)/2)
where lambda = long - lambda0
phi1 = acos $ 2 / pi
alpha = acos $ cos lat * cos (lambda/2)
lambda0 = radians (-77.0369)
mercator :: Floating a => (a, a) -> (a, a)
mercator (long, lat) = (long - meridian, asinh (tan lat))
where meridian = radians (-98.5795)
bonne :: Floating a
=> a
-> a
-> (a, a)
-> (a, a)
bonne phi1 meridian (long, lat) = (rho * sin e, cot phi1 - rho * cos e)
where rho = cot phi1 + phi1 - lat
e = (long - meridian) * cos lat / rho
cot = (1/) . tan
albers :: Floating a => (a, a)
-> ((a, a) -> (a, a))
albers referencePoint (long, lat) = (rho * sin theta, rho' - rho * cos theta)
where n = (sin phi1 + sin phi2) / 2
theta = n * (long - referenceLong)
c = cos phi1^(2 :: Int) + 2 * n * sin phi1
rho = sqrt (c - 2 * n * sin lat) / n
rho' = sqrt (c - 2 * n * sin referenceLat) / n
phi1 = radians 20
phi2 = radians 50
(referenceLong, referenceLat) = referencePoint
project :: (Floating a, Functor f) => ((a, a) -> (a, a))
-> f (a, a)
-> f (a, a)
project f = fmap (f . toRadians)
areaTriangle :: Floating a => a
-> (a, a)
-> (a, a)
-> (a, a)
-> a
areaTriangle r x1 x2 x3 = r^(2 :: Int) * e
where 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 :: (Floating a) => [(a, a)] -> a
relativeCompactness = (*scale) . compactness
where scale = 1/4*pi
compactness :: (Floating a)
=> [(a, a)]
-> a
compactness p = areaPolygon 1 p/(perimeterPolygon 1 p^(2 :: Int))
areaConvex :: (Floating a)
=> a
-> [(a, a)]
-> a
areaConvex r (base1:base2:pts) = fst $ foldr stepArea (0,base2) pts
where stepArea point (sum', base) = (sum' + areaTriangle r base1 base point, point)
areaConvex _ _ = error "attempted to take area of polygon with < 3 points"
areaPolygon :: Floating a
=> a
-> [(a, a)]
-> a
areaPolygon r = (*factor) . areaPolyRectangular . fmap (bonne (radians 25) (snd washingtonDC) . toRadians)
where factor = 1717856/4.219690791828533e-2 * ((r / 6371) ^ (2 :: Int))
perimeterPolygon :: Floating a
=> a
-> [(a, a)]
-> a
perimeterPolygon r [x1, x2] = distance r x1 x2
perimeterPolygon r (x1:x2:points) = perimeterPolygon r (x2:points) + distance r x1 x2
perimeterPolygon _ _ = error "Attempted to take perimeter of polygon with no points"
areaPolyRectangular :: (Floating a) => [(a, a)] -> a
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 :: Floating a
=> a
-> (a, a)
-> (a, a)
-> a
distance r = (*r) .* on centralAngle toRadians
centralAngle :: Floating a => (a, a) -> (a, a) -> a
centralAngle (long1, lat1) (long2, lat2) =
acos $ sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (long1 - long2)