{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}

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

import Geodetics.Ellipsoids
import Geodetics.Geodetic
import Geodetics.Grid

import qualified Data.Stream as Stream

-- | 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).
--
-- The calculations here are based on \"Transverse Mercator Projection: Constants, Formulae and Methods\"
-- by the Ordnance Survey, March 1983.
-- Retrieved from http://www.threelittlemaids.co.uk/magdec/transverse_mercator_projection.pdf
data GridTM e = GridTM {
   forall e. GridTM e -> Geodetic e
trueOrigin :: Geodetic e,
      -- ^ A point on the line where the projection touches the ellipsoid (altitude is ignored).
   forall e. GridTM e -> GridOffset
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.
   forall e. GridTM e -> Double
gridScale :: 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.
   forall e. GridTM e -> Double
gridN1, forall e. GridTM e -> Double
gridN2, forall e. GridTM e -> Double
gridN3, forall e. GridTM e -> Double
gridN4 :: Double
} deriving (Int -> GridTM e -> ShowS
[GridTM e] -> ShowS
GridTM e -> String
(Int -> GridTM e -> ShowS)
-> (GridTM e -> String) -> ([GridTM e] -> ShowS) -> Show (GridTM e)
forall e. Ellipsoid e => Int -> GridTM e -> ShowS
forall e. Ellipsoid e => [GridTM e] -> ShowS
forall e. Ellipsoid e => GridTM e -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall e. Ellipsoid e => Int -> GridTM e -> ShowS
showsPrec :: Int -> GridTM e -> ShowS
$cshow :: forall e. Ellipsoid e => GridTM e -> String
show :: GridTM e -> String
$cshowList :: forall e. Ellipsoid e => [GridTM e] -> ShowS
showList :: [GridTM e] -> ShowS
Show)


-- | Create a Transverse Mercator grid.
mkGridTM :: (Ellipsoid e) =>
   Geodetic e               -- ^ True origin.
   -> GridOffset            -- ^ Vector from true origin to false origin.
   -> Double                -- ^ Scale factor.
   -> GridTM e
mkGridTM :: forall e.
Ellipsoid e =>
Geodetic e -> GridOffset -> Double -> GridTM e
mkGridTM Geodetic e
origin GridOffset
offset Double
sf =
   GridTM {trueOrigin :: Geodetic e
trueOrigin = Geodetic e
origin,
           falseOrigin :: GridOffset
falseOrigin = GridOffset
offset,
           gridScale :: Double
gridScale = Double
sf,
           gridN1 :: Double
gridN1 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
5Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
5Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3,
           gridN2 :: Double
gridN2 = Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
21Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
8) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3,
           gridN3 :: Double
gridN3 = (Double
15Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
8) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3),
           gridN4 :: Double
gridN4 = (Double
35Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
nDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3
        }
    where
       f :: Double
f = e -> Double
forall e. Ellipsoid e => e -> Double
flattening (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
origin
       n :: Double
n = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
f)  -- Equivalent to (a-b)/(a+b) where b = (1-f)*a


-- | Equation C3 from reference [1].
m :: (Ellipsoid e) => GridTM e -> Double -> Double
m :: forall e. Ellipsoid e => GridTM e -> Double -> Double
m GridTM e
grid Double
lat = Double
bF0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (GridTM e -> Double
forall e. GridTM e -> Double
gridN1 GridTM e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dLat
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridTM e -> Double
forall e. GridTM e -> Double
gridN2 GridTM e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
dLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
sLat
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridTM e -> Double
forall e. GridTM e -> Double
gridN3 GridTM e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dLat) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sLat)
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridTM e -> Double
forall e. GridTM e -> Double
gridN4 GridTM e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dLat) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sLat))
   where
      dLat :: Double
dLat = Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
- Geodetic e -> Double
forall e. Geodetic e -> Double
latitude (GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin GridTM e
grid)
      sLat :: Double
sLat = Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude (GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin GridTM e
grid)
      bF0 :: Double
bF0 = e -> Double
forall e. Ellipsoid e => e -> Double
minorRadius (GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid) Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridTM e -> Double
forall e. GridTM e -> Double
gridScale GridTM e
grid


instance (Ellipsoid e) => GridClass (GridTM e) e where
   fromGrid :: GridPoint (GridTM e) -> Geodetic e
fromGrid GridPoint (GridTM e)
p = -- trace traceMsg $
      Double -> Double -> Double -> e -> Geodetic e
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic
         (Double
lat' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_VII Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_VIII Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_IX)
         (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude (GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin GridTM e
grid)
               Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
east' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_X Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_XI Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_XII Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
east' Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_XIIa)
         (GridPoint (GridTM e) -> Double
forall r. GridPoint r -> Double
altGP GridPoint (GridTM e)
p)
         (GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid)
      where
         GridPoint Double
east' Double
north' Double
_ GridTM e
_ = GridTM e -> GridOffset
forall e. GridTM e -> GridOffset
falseOrigin GridTM e
grid GridOffset -> GridPoint (GridTM e) -> GridPoint (GridTM e)
forall g. GridOffset -> GridPoint g -> GridPoint g
`applyOffset` GridPoint (GridTM e)
p
         lat' :: Double
lat' = (Double, Double) -> Double
forall a b. (a, b) -> a
fst ((Double, Double) -> Double) -> (Double, Double) -> Double
forall a b. (a -> b) -> a -> b
$ Stream (Double, Double) -> (Double, Double)
forall a. Stream a -> a
Stream.head (Stream (Double, Double) -> (Double, Double))
-> Stream (Double, Double) -> (Double, Double)
forall a b. (a -> b) -> a -> b
$ ((Double, Double) -> Bool)
-> Stream (Double, Double) -> Stream (Double, Double)
forall a. (a -> Bool) -> Stream a -> Stream a
Stream.dropWhile ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-5) (Double -> Bool)
-> ((Double, Double) -> Double) -> (Double, Double) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Num a => a -> a
abs (Double -> Double)
-> ((Double, Double) -> Double) -> (Double, Double) -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Double) -> Double
forall a b. (a, b) -> b
snd)
               (Stream (Double, Double) -> Stream (Double, Double))
-> Stream (Double, Double) -> Stream (Double, Double)
forall a b. (a -> b) -> a -> b
$ Stream (Double, Double) -> Stream (Double, Double)
forall a. Stream a -> Stream a
Stream.tail (Stream (Double, Double) -> Stream (Double, Double))
-> Stream (Double, Double) -> Stream (Double, Double)
forall a b. (a -> b) -> a -> b
$ ((Double, Double) -> (Double, Double))
-> (Double, Double) -> Stream (Double, Double)
forall a. (a -> a) -> a -> Stream a
Stream.iterate (Double, Double) -> (Double, Double)
forall {b}. (Double, b) -> (Double, Double)
next (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude (Geodetic e -> Double) -> Geodetic e -> Double
forall a b. (a -> b) -> a -> b
$ GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin GridTM e
grid, Double
1)
            where
               next :: (Double, b) -> (Double, Double)
next (Double
phi, b
_) = let delta :: Double
delta = Double
north' Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridTM e -> Double -> Double
forall e. Ellipsoid e => GridTM e -> Double -> Double
m GridTM e
grid Double
phi in (Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
delta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
aF0, Double
delta)
         -- Terms defined in [1]
         term_VII :: Double
term_VII  = Double
tanLat Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
         term_VIII :: Double
term_VIII = (Double
tanLat Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
24 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3))  Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eta2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta2)
         term_IX :: Double
term_IX   = (Double
tanLat Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
720 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
61 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
90 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
45 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4)
         term_X :: Double
term_X    = Double
1                                                                 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
cosLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
         term_XI :: Double
term_XI   = (Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)                                       Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3)
         term_XII :: Double
term_XII  = ( Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+  Double
28 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+   Double
24 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4)                     Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
120 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5)
         term_XIIa :: Double
term_XIIa = (Double
61 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
662 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1320 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
720 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_6) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
5040 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_7)

         -- Trace message for debugging. Uncomment this code to inspect intermediate values.
         {-
         traceMsg = concat [
            "lat' = ", show lat', "\n",
            "v    = ", show v, "\n",
            "rho  = ", show rho, "\n",
            "eta2 = ", show eta2, "\n",
            "VII  = ", show term_VII, "\n",
            "VIII = ", show term_VIII, "\n",
            "IX   = ", show term_IX, "\n",
            "X    = ", show term_X, "\n",
            "XI   = ", show term_XI, "\n",
            "XII  = ", show term_XII, "\n",
            "XIIa = ", show term_XIIa, "\n"]
         -}
         sinLat :: Double
sinLat = Double -> Double
forall a. Floating a => a -> a
sin Double
lat'
         cosLat :: Double
cosLat = Double -> Double
forall a. Floating a => a -> a
cos Double
lat'
         tanLat :: Double
tanLat = Double -> Double
forall a. Floating a => a -> a
tan Double
lat'
         sinLat2 :: Double
sinLat2 = Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat
         v :: Double
v = Double
aF0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat2)
         rho :: Double
rho = Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat2)
         eta2 :: Double
eta2 = Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1

         aF0 :: Double
aF0 = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius (GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid) Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridTM e -> Double
forall e. GridTM e -> Double
gridScale GridTM e
grid
         e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid
         grid :: GridTM e
grid = GridPoint (GridTM e) -> GridTM e
forall r. GridPoint r -> r
gridBasis GridPoint (GridTM e)
p

   toGrid :: GridTM e -> Geodetic e -> GridPoint (GridTM e)
toGrid GridTM e
grid Geodetic e
geo = -- trace traceMsg $ 
      GridOffset -> GridPoint (GridTM e) -> GridPoint (GridTM e)
forall g. GridOffset -> GridPoint g -> GridPoint g
applyOffset (GridOffset
off  GridOffset -> GridOffset -> GridOffset
forall a. Monoid a => a -> a -> a
`mappend` GridOffset -> GridOffset
offsetNegate (GridTM e -> GridOffset
forall e. GridTM e -> GridOffset
falseOrigin GridTM e
grid)) (GridPoint (GridTM e) -> GridPoint (GridTM e))
-> GridPoint (GridTM e) -> GridPoint (GridTM e)
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> GridTM e -> GridPoint (GridTM e)
forall r. Double -> Double -> Double -> r -> GridPoint r
GridPoint Double
0 Double
0 Double
0 GridTM e
grid
      where
         v :: Double
v = Double
aF0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat2)
         rho :: Double
rho = Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat2)
         eta2 :: Double
eta2 = Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
         off :: GridOffset
off = Double -> Double -> Double -> GridOffset
GridOffset
                  (Double
dLong Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_IV
                   Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dLong Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_V
                   Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dLong Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_VI)
                  (GridTM e -> Double -> Double
forall e. Ellipsoid e => GridTM e -> Double -> Double
m GridTM e
grid Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dLong Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_II
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dLong Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_III
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dLong Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
term_IIIa)
                  Double
0
         -- Terms defined in [1].
         term_II :: Double
term_II   = (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat
         term_III :: Double
term_III  = (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
24) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta2)
         term_IIIa :: Double
term_IIIa = (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
720) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
61 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
58 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4)
         term_IV :: Double
term_IV   = Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat
         term_V :: Double
term_V    = (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)
         term_VI :: Double
term_VI   = (Double
vDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
120) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_5
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
18 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
                              Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta2
                              Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
58 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanLat Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta2)

         -- Trace message for debugging. Uncomment this code to inspect 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 :: Double
lat = Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
         long :: Double
long = Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo
         dLong :: Double
dLong = Double
long Double -> Double -> Double
forall a. Num a => a -> a -> a
- Geodetic e -> Double
forall e. Geodetic e -> Double
longitude (GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin GridTM e
grid)
         sinLat :: Double
sinLat = Double -> Double
forall a. Floating a => a -> a
sin Double
lat
         cosLat :: Double
cosLat = Double -> Double
forall a. Floating a => a -> a
cos Double
lat
         tanLat :: Double
tanLat = Double -> Double
forall a. Floating a => a -> a
tan Double
lat
         sinLat2 :: Double
sinLat2 = Double
sinLat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat
         aF0 :: Double
aF0 = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius (GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid) Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridTM e -> Double
forall e. GridTM e -> Double
gridScale GridTM e
grid
         e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ GridTM e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridTM e
grid

   gridEllipsoid :: GridTM e -> e
gridEllipsoid = Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid (Geodetic e -> e) -> (GridTM e -> Geodetic e) -> GridTM e -> e
forall b c a. (b -> c) -> (a -> b) -> a -> c
. GridTM e -> Geodetic e
forall e. GridTM e -> Geodetic e
trueOrigin