{-# LANGUAGE CPP #-}
{-|
  To represent an Irish National Grid reference.

  Projection: Transverse Mercator
  Reference ellipsoid: Modified Airy
  Units: metres
  Origin: 53°30'N, 8°W
  False co-ordinates of origin: 200000m east, 250000m north
-}
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 -- ^ The easting in metres relative to the origin of the British National Grid.
                         , northing :: Double -- ^ The northing in metres relative to the origin of the British National Grid.
                         , 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

{-|
  Create a new Ordnance Survey grid reference using a given easting and
  northing. The easting and northing must be in metres and must be relative
  to the origin of the British National Grid.
-}
mkIrishRef :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 400000.0.
              -> Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than or equal to 500000.0.
              -> Except String IrishRef -- ^ Throws an exception if either the easting or the northing are invalid.
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 }


{-|
  Take a string formatted as a six-figure OS grid reference (e.g. "TG514131")
  and create a new OSRef object that represents that grid reference. The
  first character must be H, N, S, O or T. The second character can be any
  uppercase character from A through Z excluding I.
-}
mkIrishRef' :: String -- ^ A String representing a six-figure Ordnance Survey grid reference in the form XY123456.
               -> Except String IrishRef -- ^ Throws an exception if ref is not of the form XY123456.
mkIrishRef' ref = do
  -- TODO 2006-02-05 : check format
  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 -- Adjust for no I
      nx = ((ch - 65) `mod` 5) * 100000
      ny = (4 - floor ((fromIntegral ch - 65) / 5)) * 100000

  mkIrishRef (fromIntegral (east + nx)) (fromIntegral (north + ny))


-- | Create an IrishRef object from the given latitude and longitude.
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


{-|
  Return a String representation of this Irish grid reference using the
  six-figure notation in the form X123456
-}
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


{-|
  Convert this Irish grid reference to a latitude/longitude pair using the
  Ireland 1965 datum. Note that, the LatLng object may need to be converted to the
  WGS84 datum depending on the application.
-}
toLatLng :: IrishRef
            -> Except String L.LatLng -- ^ To represent Irish grid reference using the Ireland 1965 datum.
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


-- | Validate the easting.
evalEasting :: Double                  -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 400000.0.
               -> Except String Double -- ^ Throws an Exception if the easting is invalid.
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)


-- | Validate the northing.
evalNorthing :: Double                  -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than or equal to 500000.0.
                -> Except String Double -- ^ Throws an Exception if the northing is invalid.
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)