module IrishRef where
#if __GLASGOW_HASKELL__ < 710
import Control.Applicative
#endif
import Control.Monad.Except
import Data.Char
import Datum
import Ellipsoid
import qualified LatLng as L
import MathExtensions
data IrishRef = IrishRef { easting :: Double
, northing :: Double
, datum :: Datum
} deriving (Show)
scaleFactor, falseOriginLatitude, falseOriginLongitude, falseOriginEasting, falseOriginNorthing :: Double
scaleFactor = 1.000035
falseOriginLatitude = 53.5
falseOriginLongitude = 8.0
falseOriginEasting = 200000.0
falseOriginNorthing = 250000.0
mkIrishRef :: Double
-> Double
-> Except String IrishRef
mkIrishRef e n = do
est <- withExcept (const "Invalid easting") (evalEasting e)
nrt <- withExcept (const "Invalid northing") (evalNorthing n)
pure IrishRef { easting = est, northing = nrt, datum = ireland1965Datum }
mkIrishRef' :: String
-> Except String IrishRef
mkIrishRef' ref = do
let
east = (read (take 3 $ drop 1 ref) :: Int) * 100
north = (read (take 3 $ drop 4 ref) :: Int) * 100
firstCh = ord (ref !! 0)
ch = if firstCh > 73 then firstCh 1 else firstCh
nx = ((ch 65) `mod` 5) * 100000
ny = (4 floor ((fromIntegral ch 65) / 5)) * 100000
mkIrishRef (fromIntegral (east + nx)) (fromIntegral (north + ny))
mkIrishRef'' :: L.LatLng -> Except String IrishRef
mkIrishRef'' (L.LatLng latitude longitude height _) = do
let
n0 = falseOriginNorthing
e0 = falseOriginEasting
phi0 = toRadians falseOriginLatitude
lambda0 = toRadians falseOriginLongitude
el = ellipsoid ireland1965Datum
a = scaleFactor * semiMajorAxis el
b = scaleFactor * semiMinorAxis el
eSquared = eccentricitySquared el
phi = toRadians latitude
lambda = toRadians longitude
n = (a b) / (a + b)
va = a * (1 eSquared * sinSquared phi) ** (0.5)
rho = a * (1 eSquared) * (1 eSquared * sinSquared phi) ** (1.5)
etaSquared = va / rho 1
m = b * ((1 + n + 1.25 * n ** 2 + 1.25 * n ** 3) * (phi phi0)
(3 * n + 3 * n ** 2 + 21.0 / 8.0 * n ** 3) * sin (phi phi0) * cos (phi + phi0)
+ (15.0 / 8.0 * n ** 2 + 15.0 / 8.0 * n ** 3) * sin(2.0 * (phi phi0)) * cos(2.0 * (phi + phi0))
(35.0 / 24.0 * n ** 3) * sin(3.0 * (phi phi0)) * cos(3.0 * (phi + phi0)))
i = m + n0
ii = va / 2.0 * sin phi * cos phi
iii = va / 24.0 * sin phi * cos phi ** 3 * (5.0 tanSquared phi + 9.0 * etaSquared)
iiia = va / 720.0 * sin phi * cos phi ** 5 * (61.0 58.0 * tanSquared phi + tan phi ** 4)
iv = va * cos(phi)
v = va / 6.0 * cos phi ** 3 * (va / rho tanSquared phi)
vi = va / 120.0 * cos phi ** 5.0 * (5.0 18.0 * tanSquared phi + tan phi ** 4.0
+ 14 * etaSquared 58 * tanSquared phi * etaSquared)
north = i + ii * (lambda lambda0) ** 2
+ iii * (lambda lambda0) ** 4
+ iiia * (lambda lambda0) ** 6
east = e0 + iv * (lambda lambda0)
+ v * (lambda lambda0) ** 3
+ vi * (lambda lambda0) ** 5
mkIrishRef east north
toSixFigureString :: IrishRef -> String
toSixFigureString (IrishRef easting northing datum) = do
let
hundredkmE = floor (easting / 100000)
hundredkmN = floor (northing / 100000)
charOffset = 4 hundredkmN
i = 65 + 5 * charOffset + hundredkmE
index = if (i >= 73) then i + 1 else i
e = floor ((easting 100000 * fromIntegral hundredkmE) / 100)
n = floor ((northing 100000 * fromIntegral hundredkmN) / 100)
chr index : compose e ++ compose n
where compose :: Int -> String
compose x = (if (x < 100) then "0" else "")
++ (if (x < 10) then "0" else "")
++ show x
toLatLng :: IrishRef
-> Except String L.LatLng
toLatLng (IrishRef easting northing datum) = do
let
n0 = falseOriginNorthing
e0 = falseOriginEasting
phi0 = toRadians falseOriginLatitude
lambda0 = toRadians falseOriginLongitude
el = ellipsoid datum
a = semiMajorAxis el
b = semiMinorAxis el
eSquared = eccentricitySquared el
n = (a b) / (a + b)
phiPrime = calcPhiPrime ((northing n0) / (a * scaleFactor) + phi0) a b n phi0 n0 northing
va = a * scaleFactor * (1 eSquared * sinSquared phiPrime) ** (0.5)
rho = a * scaleFactor * (1 eSquared) * (1 eSquared * sinSquared phiPrime) ** (1.5)
etaSquared = va / rho 1
vii = tan phiPrime / (2 * rho * va)
viii = tan phiPrime / (24 * rho * va ** 3)
* (5 + 3 * tanSquared phiPrime + etaSquared 9 * tanSquared phiPrime * etaSquared)
ix = tan phiPrime / (720 * rho * va ** 5)
* (61 + 90 * tanSquared phiPrime + 45 * tanSquared phiPrime ** 2)
x = sec phiPrime / va
xi = sec phiPrime / (6 * va ** 3) * (va / rho + 2 * tanSquared phiPrime)
xii = sec phiPrime / (120 * va ** 5)
* (5 + 28 * tanSquared phiPrime + 24 * tanSquared phiPrime ** 2)
xiia = sec phiPrime / (5040 * va ** 7)
* (61 + 662 * tanSquared phiPrime + 1320 * tanSquared phiPrime ** 2
+ 720 * tanSquared phiPrime ** 3)
phi = phiPrime vii * (easting e0) ** 2 + viii * (easting e0) ** 4.0 ix * (easting e0) ** 6.0
lambda = lambda0 + x * (easting e0) xi * (easting e0) ** 3.0 + xii * (easting e0) ** 5.0 xiia * (easting e0) ** 7.0
L.mkLatLng (toDegrees phi) (toDegrees lambda) 0 wgs84Datum
where calcPhiPrime :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
calcPhiPrime phi a b n phi0 n0 northing = do
let m = b * scaleFactor
* ((1 + n + 1.25 * n ** 2 + 1.25 * n ** 3) * (phi phi0)
(3 * n + 3 * n ** 2 + 21.0 / 8.0 * n ** 3) * sin (phi phi0) * cos (phi + phi0)
+ (15.0 / 8.0 * n ** 2 + 15.0 / 8.0 * n ** 3) * sin(2.0 * (phi phi0)) * cos(2.0 * (phi + phi0))
(35.0 / 24.0 * n ** 3) * sin(3.0 * (phi phi0)) * cos(3.0 * (phi + phi0)))
delta = northing n0 m
phiPrime = phi + delta / (a * scaleFactor)
if (delta >= 0.001) then calcPhiPrime phiPrime a b n phi0 n0 northing
else phiPrime
evalEasting :: Double
-> Except String Double
evalEasting e | e < 0.0 || e >= 400000.0 = throwError ("Easting (" ++ show e ++ ") is invalid. Must be greater than or equal to 0.0 and less than 400000.0.")
| otherwise = pure (e)
evalNorthing :: Double
-> Except String Double
evalNorthing n | n < 0.0 || n > 500000.0 = throwError ("Northing (" ++ show n ++ ") is invalid. Must be greather than or equal to 0.0 and less than or equal to 500000.0.")
| otherwise = pure (n)