module Chorale.Geo.Coordinates (
radiusEarth,
degreeToRad,
radToDegree,
Coordinates(..),
CoordinatesCartesian(..),
CoordinatesWGS84(..)) where
import Chorale.Common
import qualified Data.Binary as B
import GHC.Generics (Generic)
radiusEarth :: Double
radiusEarth = 6371.009
degreeToRad :: Double -> Double
degreeToRad d = d * pi / 180
radToDegree :: Double -> Double
radToDegree r = r * 180 / pi
norm :: Floating a => [a] -> a
norm = sqrt . sum . uncurry (zipWith (*)) . make2
class (Ord c, Eq c, Show c) => Coordinates c where
toTuple :: c -> (Double, Double)
distance :: c -> c -> Double
azimuth :: c -> c -> Double
newtype CoordinatesCartesian = CoordinatesCartesian (Double, Double) deriving (Ord, Eq, Show, Generic)
newtype CoordinatesWGS84 = CoordinatesWGS84 (Double, Double) deriving (Ord, Eq, Show, Generic)
instance B.Binary CoordinatesCartesian
instance Coordinates CoordinatesCartesian where
toTuple (CoordinatesCartesian x) = x
distance = curry $ norm . uncurry (zipWith ()) . map12 (tupleToList2 . toTuple)
azimuth (CoordinatesCartesian (x1, y1)) (CoordinatesCartesian (x2, y2)) = alpha (x2 x1, y2 y1) where
alpha (x, y)
| x == 0 && y == 0 = error "azimuth cannot be computed between two identical vectors"
| az < 0 = az + 2 * pi
| otherwise = az where
az = atan2 y x + pi / 2
instance Coordinates (Double, Double) where
toTuple = id
distance = curry $ uncurry distance . map12 CoordinatesCartesian
azimuth = curry $ uncurry azimuth . map12 CoordinatesCartesian
instance B.Binary CoordinatesWGS84
instance Coordinates CoordinatesWGS84 where
toTuple (CoordinatesWGS84 x) = x
distance cs1 cs2 = radiusEarth * sqrt (diffAngle lat2 lat1 ** 2 + (cos ((lat2 + lat1) / 2) * diffAngle lon2 lon1)**2) where
[lat1, lon1] = map degreeToRad . tupleToList2 . toTuple $ cs1
[lat2, lon2] = map degreeToRad . tupleToList2 . toTuple $ cs2
diffAngle y x
| abs (y x) < pi / 2 = y x
| otherwise = diffAngle y (x + signum (y x) * pi)
azimuth (CoordinatesWGS84 (lat1', lon1')) (CoordinatesWGS84 (lat2', lon2'))
| lat1' lat2' == 0 && lon1' lon2' == 0 = error "azimuth cannot be computed between two identical vectors"
| az < 0 = az + 2 * pi
| otherwise = az where
az = atan2 (cos u2 * sin lambda) (cos u1 * sin u2 sin u1 * cos u2 * cos lambda)
lat1 = degreeToRad lat1'
lat2 = degreeToRad lat2'
lon1 = degreeToRad lon1'
lon2 = degreeToRad lon2'
f = 1 / 298.257223563
u1 = atan ((1 f) * tan lat1)
u2 = atan ((1 f) * tan lat2)
l = lon2 lon1
lambda = recursiveFunction l
recursiveFunction x
| abs (x x') < 10e-12 = x'
| otherwise = recursiveFunction x' where
sins = sqrt $ (cos u2 * sin x)**2 + (cos u1 * sin u2 sin u1 * cos u2 * cos x)**2
coss = sin u1 * sin u2 + cos u1 * cos u2 * cos x
s = atan (sins / coss)
sina = (cos u1 * cos u2 * sin x) / sins
cos2a = 1 sina**2
cos2sm = sqrt cos2a 2 * sin u1 * sin u2 / cos2a
c = f / 16 * cos2a * (4 + f * (4 3 * cos2a))
x' = l + (1 c) * f * sina * (s + c * sins * (cos2sm + c * coss * (1 + 2 * cos2sm**2)))