{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}

module Geodetics.TransverseMercator(
   GridTM (trueOrigin, falseOrigin, gridScale),
   mkGridTM
) where

import Data.Function
import Data.Monoid
import Geodetics.Ellipsoids
import Geodetics.Geodetic
import Geodetics.Grid
import Numeric.Units.Dimensional.Prelude hiding ((.))
import Prelude ()

-- | A Transverse Mercator projection gives an approximate mapping of the ellipsoid on to a 2-D grid. It models
-- a sheet curved around the ellipsoid so that it touches it at one north-south line (hence making it part of
-- a slightly elliptical cylinder).
data GridTM e = GridTM {
   trueOrigin :: Geodetic e,
      -- ^ A point on the line where the projection touches the ellipsoid (altitude is ignored).
   falseOrigin :: GridOffset,
      -- ^ The grid position of the true origin. Used to avoid negative coordinates over 
      -- the area of interest. The altitude gives a vertical offset from the ellipsoid.
   gridScale :: Dimensionless Double,
      -- ^ A scaling factor that balances the distortion between the east & west edges and the middle 
      -- of the projection.
      
   -- Remaining elements are memoised parameters computed from the ellipsoid underlying the true origin.
   gridN1, gridN2, gridN3, gridN4 :: Dimensionless Double
} deriving (Show)


-- | Create a Transverse Mercator grid.
mkGridTM :: (Ellipsoid e) => 
   Geodetic e               -- ^ True origin.
   -> GridOffset            -- ^ Vector from true origin to false origin.
   -> Dimensionless Double  -- ^ Scale factor.
   -> GridTM e
mkGridTM origin offset sf =
   GridTM {trueOrigin = origin,
           falseOrigin = offset,
           gridScale = sf,
           gridN1 = _1 + n + (_5/_4) * n^pos2 + (_5/_4) * n^pos3,
           gridN2 = _3 * n + _3 * n^pos2 + ((21*~one)/_8) * n^pos3,
           gridN3 = ((15*~one)/_8) * (n^pos2 + n^pos3),
           gridN4 = ((35*~one)/(24*~one)) * n^pos3
        }
    where 
       f = flattening $ ellipsoid origin
       n = f / (_2-f)  -- Equivalent to (a-b)/(a+b) where b = (1-f)*a




-- | Equation C3 from reference [1].
m :: (Ellipsoid e) => GridTM e -> Dimensionless Double -> Length Double
m grid lat = bF0 * (gridN1 grid * dLat 
                    - gridN2 grid * sin dLat * cos sLat
                    + gridN3 grid * sin (_2 * dLat) * cos (_2 * sLat) 
                    - gridN4 grid * sin (_3 * dLat) * cos (_3 * sLat))
   where
      dLat = lat - latitude (trueOrigin grid)
      sLat = lat + latitude (trueOrigin grid)
      bF0 = minorRadius (gridEllipsoid grid) * gridScale grid


instance (Ellipsoid e) => GridClass (GridTM e) e where
   fromGrid p = Geodetic
      (lat' - east' ^ pos2 * tanLat / (_2 * rho * v)  -- Term VII
            + east' ^ pos4 * (tanLat / ((24 *~ one) * rho * v ^ pos3)) 
                           * (_5 + _3 * tanLat ^ pos2 + eta2 - _9 * tanLat ^ pos2 * eta2)  -- Term VIII
            - east' * east' ^ pos5 * (tanLat / ((720 *~ one) * rho * v ^ pos5))
                           * (61 *~ one + (90 *~ one) * tanLat ^ pos2 + (45 *~ one) * tanLat ^ pos4)) -- Term IX
      (longitude (trueOrigin grid) 
            + east' / (cosLat * v)  -- Term X
            - (east' ^ pos3 / (_6 * cosLat * v ^ pos3)) * (v / rho + _2 * tanLat ^ pos2)  -- Term XI
            + (east' ^ pos5 / ((120 *~ one) * cosLat * v ^ pos5)) 
                 * (_5 + (28 *~ one) * tanLat ^ pos2  + (24 *~ one) * tanLat ^ pos4)  -- Term XII
            - (east' ^ pos5 * east' ^ pos2 / ((5040 *~ one) * cosLat * v * v * v ^ pos5))
                 * ((61 *~ one) + (662 *~ one) * tanLat ^ pos2 + (1320 *~ one) * tanLat ^ pos4 + (720 *~ one) * tanLat * tanLat ^ pos5)) -- Term XIIa
     (0 *~ meter) (gridEllipsoid grid)
            
            
      where
         GridPoint east' north' _ _ = falseOrigin grid `applyOffset` p
         lat' = fst $ head $ dropWhile ((> 0.01 *~ milli meter) . snd) 
               $ tail $ iterate next (latitude $ trueOrigin grid, 1 *~ meter) 
            where
               next (phi, _) = let delta = north' - m grid phi in (phi + delta / aF0, delta) 
               -- head and tail are safe because iterate returns an infinite list.
          
         sinLat = sin lat'
         cosLat = cos lat'
         tanLat = tan lat'
         sinLat2 = sinLat ^ pos2
         v = aF0 / sqrt (_1 - e2 * sinLat2)
         rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one)
         eta2 = v / rho - _1
               
               
         aF0 = majorRadius (gridEllipsoid grid) * gridScale grid
         e2 = eccentricity2 $ gridEllipsoid grid
         grid = gridBasis p
         
   toGrid grid geo = applyOffset (off  `mappend` (offsetNegate $ falseOrigin grid)) $ 
                     GridPoint _0 _0 _0 grid
      where
         v = aF0 / sqrt (_1 - e2 * sinLat2)
         rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one)
         eta2 = v / rho - _1
         off = GridOffset
                  (dLong * term_IV
                   + dLong ^ pos3 * term_V
                   + dLong ^ pos5 * term_VI)
                  (m grid lat + dLong ^ pos2 * term_II
                     + dLong ^ pos4 * term_III 
                     + dLong * dLong ^ pos5 * term_IIIa)
                  (0 *~ meter)
         -- Terms defined in [1].
         term_II   = (v/_2) * sinLat * cosLat
         term_III  = (v/(24*~one)) * sinLat * cosLat ^ pos3 
                     * (_5 - tanLat ^ pos2 + _9 * eta2)
         term_IIIa = (v/(720*~one)) * sinLat * cosLat ^ pos5 
                     * ((61 *~ one) - (58 *~ one) * tanLat ^ pos2 + tanLat ^ pos4)
         term_IV   = v * cosLat
         term_V    = (v/_6) * cosLat ^ pos3 * (v/rho - tanLat ^ pos2)
         term_VI   = (v/(120*~one)) * cosLat ^ pos5 
                     * (_5 - (18*~one) * tanLat ^ pos2 
                              + tanLat ^ pos4 + (14*~one) * eta2
                              - (58*~one) * tanLat ^ pos2 * eta2)
         {- 
         -- Trace message for debugging. Uncomment this code for easy access to intermediate values.
         traceMsg = concat [
            "v    = ", show v, "\n",
            "rho  = ", show rho, "\n",
            "eta2 = ", show eta2, "\n",
            "M    = ", show $ m grid lat, "\n",
            "I    = ", show $ m grid lat + deltaNorth (falseOrigin grid), "\n",
            "II   = ", show term_II, "\n",
            "III  = ", show term_III, "\n",
            "IIIa = ", show term_IIIa, "\n",
            "IV   = ", show term_IV, "\n",
            "V    = ", show term_V, "\n",
            "VI   = ", show term_VI, "\n"]
         -}
         -- Common subexpressions
         lat = latitude geo
         long = longitude geo
         dLong = long - longitude (trueOrigin grid)
         sinLat = sin lat
         cosLat = cos lat
         tanLat = tan lat
         sinLat2 = sinLat ^ pos2
         aF0 = (majorRadius $ gridEllipsoid grid) * gridScale grid
         e2 = eccentricity2 $ gridEllipsoid grid
         
   gridEllipsoid = ellipsoid . trueOrigin