module LatLng where
#if __GLASGOW_HASKELL__ < 710
import Control.Applicative
#endif
import Control.Monad.Except
import Datum
import Ellipsoid
import MathExtensions
data LatLng = LatLng { latitude :: Double
, longitude :: Double
, height :: Double
, datum :: Datum
} deriving (Show)
mkLatLng :: Double -> Double -> Double -> Datum -> Except String LatLng
mkLatLng lat lng h dtm = do
lt <- withExcept (const $ "Latitude (" ++ show lat ++ ") is invalid. Must be between -90.0 and 90.0 inclusive.") (validateLatitude lat)
ln <- withExcept (const $ "Longitude (" ++ show lng ++ ") is invalid. Must be between -180.0 and 180.0 inclusive.") (validateLongitude lng)
pure LatLng { latitude = lt, longitude = ln, height = h, datum = dtm }
where
validateLatitude :: Double -> Except String Double
validateLatitude l
| l < 90.0 || l > 90.0 = throwError "Invalid latitude"
| otherwise = pure l
validateLongitude :: Double -> Except String Double
validateLongitude l
| l < 180.0 || l > 180.0 = throwError "Invalid longitude"
| otherwise = pure l
calcPhiN' :: Double -> Double -> Double -> Double -> Double -> Double
calcPhiN' zB eSquared2 a2 p phiN = do
let v = a2 / sqrt (1 eSquared2 * sinSquared phiN)
atan $ (zB + eSquared2 * v * sin phiN) / p
toWGS84 :: LatLng -> LatLng
toWGS84 (LatLng latitude longitude height datum) = do
let
a = semiMajorAxis airy1830Ellipsoid
eSquared = eccentricitySquared airy1830Ellipsoid
phi = toRadians latitude
lambda = toRadians longitude
v = a / sqrt(1 eSquared * sinSquared phi)
x = v * cos phi * cos lambda
y = v * cos phi * sin lambda
z = (1 eSquared) * v * sin phi
tx = 446.448 :: Double
ty = 125.157 :: Double
tz = 542.060 :: Double
s = 0.0000204894 :: Double
rx = toRadians (0.00004172222) :: Double
ry = toRadians (0.00006861111) :: Double
rz = toRadians (0.00023391666) :: Double
xB = tx + x * (1 + s) + (rx) * y + ry * z
yB = ty + rz * x + y * (1 + s) + (rx) * z
zB = tz + (ry) * x + rx * y + z * (1 + s)
a2 = semiMajorAxis wgs84Ellipsoid
eSquared2 = eccentricitySquared wgs84Ellipsoid
p = sqrt (xB ** 2 + yB ** 2)
phiN = atan(zB / (p * (1 eSquared2)))
calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN
LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum
toDatum :: LatLng -> Datum -> LatLng
toDatum ll@(LatLng latitude longitude height datum) newDatum
| datum /= wgs84Datum && newDatum /= wgs84Datum = mkToDatum ll newDatum 1
| datum == wgs84Datum && newDatum /= wgs84Datum = mkToDatum ll newDatum (1)
| datum /= wgs84Datum && newDatum == wgs84Datum = ll
| otherwise = mkToDatum ll newDatum 1
where mkToDatum :: LatLng -> Datum -> Double -> LatLng
mkToDatum ll@(LatLng latitude longitude height datum) newDatum invert = do
let
a = semiMajorAxis $ ellipsoid $ datum
eSquared = eccentricitySquared $ ellipsoid $ datum
phi = toRadians latitude
lambda = toRadians longitude
v = a / sqrt(1 eSquared * sinSquared phi)
x = (v + height) * cos phi * cos lambda
y = (v + height) * cos phi * sin lambda
z = ((1 eSquared) * v + height) * sin phi
dx2 = invert * dx newDatum
dy2 = invert * dy newDatum
dz2 = invert * dz newDatum
ds2 = invert * ds newDatum / 1000000.0
rx2 = invert * toRadians (rx newDatum / 3600.0)
ry2 = invert * toRadians (ry newDatum / 3600.0)
rz2 = invert * toRadians (rz newDatum / 3600.0)
sc = 1 + ds2
xB = dx2 + sc * (x + (rx2) * y + ry2 * z)
yB = dy2 + sc * (rz2 * x + y + (rx2) * z)
zB = dz2 + sc * ((ry2) * x + rx2 * y + z)
a2 = semiMajorAxis $ ellipsoid $ newDatum
eSquared2 = eccentricitySquared $ ellipsoid $ newDatum
p = sqrt(xB ** 2 + yB ** 2)
phiN = atan (zB / (p * (1 eSquared2)))
calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN
LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum
toOSGB36 :: LatLng -> LatLng
toOSGB36 (LatLng latitude longitude height datum) = do
let
a = semiMajorAxis wgs84Ellipsoid
eSquared = eccentricitySquared wgs84Ellipsoid
phi = toRadians latitude
lambda = toRadians longitude
v = a / sqrt(1 eSquared * sinSquared phi)
x = v * cos phi * cos lambda
y = v * cos phi * sin lambda
z = (1 eSquared) * v * sin phi
tx = 446.448 :: Double
ty = 125.157 :: Double
tz = 542.060 :: Double
s = 0.0000204894 :: Double
rx = toRadians (0.00004172222) :: Double
ry = toRadians (0.00006861111) :: Double
rz = toRadians (0.00023391666) :: Double
xB = tx + x * (1 + s) + (rx) * y + ry * z
yB = ty + rz * x + y * (1 + s) + (rx) * z
zB = tz + (ry) * x + rx * y + z * (1 + s)
a2 = semiMajorAxis airy1830Ellipsoid
eSquared2 = eccentricitySquared airy1830Ellipsoid
p = sqrt (xB ** 2 + yB ** 2)
phiN = atan(zB / (p * (1 eSquared2)))
calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN
LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum
distance :: LatLng -> LatLng -> Double
distance (LatLng latitudeFrom longitudeFrom _ _) (LatLng latitudeTo longitudeTo _ _) = do
let
er = 6366.707 :: Double
latFrom = toRadians latitudeFrom
latTo = toRadians latitudeTo
lngFrom = toRadians longitudeFrom
lngTo = toRadians longitudeTo
acos (sin latFrom * sin latTo + cos latFrom * cos latTo * cos (lngTo lngFrom)) * er
distanceMiles :: LatLng -> LatLng -> Double
distanceMiles from to = (distance from to) / 1.609344
latitudeDegrees, longitudeDegrees :: LatLng -> Int
latitudeDegrees (LatLng l _ _ _) = calcDegrees l
longitudeDegrees (LatLng _ l _ _) = calcDegrees l
calcDegrees :: Double -> Int
calcDegrees l = do
let deg = floor l
minx = (l fromIntegral deg) :: Double
case l of _ | l < 0 && minx /= 0 -> (deg :: Int) + 1
| otherwise -> deg :: Int
latitudeMinutes, longitudeMinutes :: LatLng -> Int
latitudeMinutes (LatLng l _ _ _) = calcMinutes l
longitudeMinutes (LatLng _ l _ _) = calcMinutes l
calcMinutes :: Double -> Int
calcMinutes l = do
let deg = floor l
minx = (l fromIntegral deg) :: Double
min | l < 0 && minx /= 0 = 1 minx
| otherwise = minx
floor (60 * min) :: Int
latitudeSeconds, longitudeSeconds :: LatLng -> Double
latitudeSeconds (LatLng l _ _ _) = calcSeconds l
longitudeSeconds (LatLng _ l _ _) = calcSeconds l
calcSeconds :: Double -> Double
calcSeconds l = do
let deg = floor l
minx = (l fromIntegral deg) :: Double
min | l < 0 && minx /= 0 = 1 minx
| otherwise = minx
60 * (60 * min fromIntegral (floor (60 * min)))