{-# LANGUAGE CPP #-}
{-|
  To represent an Ordnance Survey of Great Britain (OSGB) grid reference.

  British National Grid
  Projection: Transverse Mercator
  Reference ellipsoid: Airy 1830
  Units: metres
  Origin: 49°N, 2°W
  False co-ordinates of origin: 400000m east, -100000m north
  
  A full reference includes a two-character code identifying a particular
  100,000m grid square. The table below shows how the two-character 100,000m
  grid squares are identified. The bottom left corner is at the false origin of
  the grid. Squares without values fall outside the boundaries of the British
  National Grid.

   km
  0
  100
  200
  300
  400
  500
  600
  700
  
  
  1200
  HL
  HM
  HN
  HO
  HP
  JL
  JM
   
  
  
  1100
  HQ
  HR
  HS
  HT
  HU
  JQ
  JR
   
  
  
  1000
  HV
  HW
  HX
  HY
  HZ
  JV
  JW
   
  
  
   900
  NA
  NB
  NC
  ND
  NE
  OA
  OB
   
  
  
   800
  NF
  NG
  NH
  NJ
  NK
  OF
  OG
  OH
  
  
   700
  NL
  NM
  NN
  NO
  NP
  OL
  OM
  ON
  
  
   600
  NQ
  NR
  NS
  NT
  NU
  OQ
  OR
  OS
  
  
   500
   
  NW
  NX
  NY
  NZ
  OV
  OW
  OX
  
  
   400
   
  SB
  SC
  SD
  SE
  TA
  TB
  TC
  
  
   300
   
  SG
  SH
  SJ
  SK
  TF
  TG
  TH
  
  
   200
   
  SM
  SN
  SO
  SP
  TL
  TM
  TN
  
  
   100
  SQ
  SR
  SS
  ST
  SU
  TQ
  TR
  TS
  
  
   0
  SV
  SW
  SX
  SY
  SZ
  TV
  TW
   
  
  Within each 100,000m square, the grid is further subdivided into 1000m
  squares. These 1km squares are shown on Ordnance Survey 1:25000 and 1:50000
  mapping as the main grid. To reference a 1km square, give the easting and
  then the northing, e.g. TR2266. In this example, TR represents the 100,000m
  square, 22 represents the easting (in km) and 66 represents the northing (in
  km). This is commonly called a four-figure grid reference.
  
  When providing local references, the 2 characters representing the 100,000m
  square are often omitted.
-}
module OSRef 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 OSRef = OSRef { 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 (Eq, Show)

data Precision = SixDigits -- ^ A six-figure representation this OSGB grid reference i.e XY123456
                 | EightDigits -- ^ A eight-figure representation this OSGB grid reference i.e XY12345678


scaleFactor, falseOriginLatitude, falseOriginLongitude, falseOriginEasting, falseOriginNorthing :: Double
scaleFactor = 0.9996012717 -- OSGB_F0
falseOriginLatitude = 49.0
falseOriginLongitude = -2.0
falseOriginEasting = 400000.0
falseOriginNorthing = -100000.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.
-}
mkOSRef :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 800000.0.
           -> Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than 1400000.0.
           -> Except String OSRef -- ^ Throws an exception if either the easting or the northing are invalid.
mkOSRef e n = do
  est <- withExcept (const "Invalid easting") (evalEasting e)
  nrt <- withExcept (const "Invalid northing") (evalNorthing n)
  pure OSRef { easting = est, northing = nrt, datum = osgb36Datum }


{-|
  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.
-}
mkOSRef' :: String -- ^ a String representing a six-figure Ordnance Survey grid reference in the form XY123456
            -> Except String OSRef -- ^ Throws an exception if ref is not of the form XY123456.
mkOSRef' ref = do
  -- TODO 2006-02-05 : check format
  let coords = findCoords (fromLabel ref 2) (fromLabel ref 5) (ref !! 0)
               where
                   fromLabel :: String -> Int -> Int
                   fromLabel r p = (read (take 3 $ drop p r) :: Int) * 100

                   findCoords :: Int -> Int -> Char -> (Int, Int)
                   findCoords e n c
                     | c == 'H' = (e, n + 1000000)
                     | c == 'N' = (e, n + 500000)
                     | c == 'O' = (e + 500000, n + 500000)
                     | c == 'T' = (e + 500000, n)

  let char2Ord = ord $ ref !! 1
  let c2 = if (char2Ord > 73) then char2Ord - 66 else char2Ord - 65

  let nx = c2 `mod` 5 * 100000
  let ny = (4 - floor ((fromIntegral c2 :: Double) / 5)) * 100000

  mkOSRef (fromIntegral (fst coords + nx)) (fromIntegral (snd coords + ny))


{-|
  Convert latitude and longitude into an OSGB (Ordnance Survey of Great Britain) grid reference.
-}
toOSRef :: L.LatLng -> Except String OSRef
toOSRef (L.LatLng latitude longitude _ _) = do
    let osgb_f0 = 0.9996012717 :: Double
        n0 = -100000.0 :: Double
        e0 = 400000.0 :: Double
        phi0 = toRadians 49.0
        lambda0 = toRadians (-2.0)

        a = semiMajorAxis airy1830Ellipsoid
        b = semiMinorAxis airy1830Ellipsoid
        eSquared = eccentricitySquared airy1830Ellipsoid

        phi = toRadians latitude
        lambda = toRadians longitude

        n = (a - b) / (a + b)
        vc = a * osgb_f0 * (1.0 - eSquared * sinSquared phi) ** (-0.5)
        rho = a * osgb_f0 * (1.0 - eSquared) * (1.0 - eSquared * sinSquared phi) ** (-1.5)
        etaSquared = vc / rho - 1.0

        cp = cos phi
        sp = sin phi
        ts = tan phi ** 2

        m = (b * osgb_f0) * (((1.0 + n + (1.25 * n * n) + (1.25 * n * n * n)) * (phi - phi0))
            - (((3 * n) + (3 * n * n) + (2.625 * n * n * n)) * sin(phi - phi0) * cos(phi + phi0))
            + (((1.875 * n * n) + (1.875 * n * n * n)) * sin(2.0 * (phi - phi0)) * cos(2.0 * (phi + phi0)))
            - (((35.0 / 24.0) * n * n * n) * sin(3.0 * (phi - phi0)) * cos(3.0 * (phi + phi0))))
        i = m + n0
        ii = (vc / 2.0) * sp * cp
        iii = (vc / 24.0) * sp * (cp ** 3) * (5.0 - ts + 9.0 * etaSquared)
        iiia = (vc / 720.0) * sp * (cp ** 5) * (61.0 - 58.0 * ts + ts ** 2)
        iv = vc * cp
        v = (vc / 6.0) * (cp ** 3) * (vc / rho - ts)
        vi = (vc / 120.0) * (cp ** 5.0) * (5.0 - 18.0 * ts + ts ** 2 + 14 * etaSquared - 58 * ts * etaSquared)

    mkOSRef
      (e0 + iv * (lambda - lambda0) + v * (lambda - lambda0) ** 3 + vi * (lambda - lambda0) ** 5)
      (i + ii * (lambda - lambda0) ** 2 + iii * (lambda - lambda0) ** 4 + iiia * (lambda - lambda0) ** 6)


getOsRefWithPrecisionOf :: Precision -> OSRef -> String
getOsRefWithPrecisionOf SixDigits (OSRef e n _) = evalOsRef 100 e n
getOsRefWithPrecisionOf EightDigits (OSRef e n _) = evalOsRef 10 e n

evalOsRef :: Double -> Double -> Double -> String
evalOsRef precision easting northing = do
  let
      hundredkmE = floor (easting / 100000)
      hundredkmN = floor (northing / 100000)

      firstLetter
        | hundredkmN < 5 = if (hundredkmE < 5) then 'S' else 'T'
        | hundredkmN < 10 = if (hundredkmE < 5) then 'N' else 'O'
        | otherwise = 'H'

      i = 85 - 5 * (hundredkmN `mod` 5) + (hundredkmE `mod` 5)
      secondLetter = chr $ if (i >= 73) then i + 1 else i

      e = floor ((easting - 100000 * fromIntegral hundredkmE) / precision)
      n = floor ((northing - 100000 * fromIntegral hundredkmN) / precision)

  firstLetter : secondLetter : 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 OSGB grid reference to a latitude/longitude pair using the
  OSGB36 datum. Note that, the LatLng object may need to be converted to the
  WGS84 datum depending on the application.
-}
toLatLng :: OSRef
            -> Except String L.LatLng -- ^ To represent OSGB grid reference using the OSGB36 datum.
toLatLng (OSRef 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 800000.0.
               -> Except String Double -- ^ Throws an Exception if the easting is invalid.
evalEasting e | e < 0.0 || e >= 800000.0 = throwError ("Easting (" ++ show e ++ ") is invalid. Must be greater than or equal to 0.0 and less than 800000.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 1400000.0.
                -> Except String Double -- ^ Throws an Exception if the northing is invalid.
evalNorthing n | n < 0.0 || n >= 1400000.0 = throwError ("Northing (" ++ show n ++ ") is invalid. Must be greather than or equal to 0.0 and less than or equal to 1400000.0.")
               | otherwise = pure (n)