Use integer powers instead of **

src/Geodetics/Ellipsoids.hs

@@ -14,11 +14,12 @@ specific area.
 -}
 
 module Geodetics.Ellipsoids (
-   -- * Conversion constants
+   -- * Useful constants
    degree,
    arcminute,
    arcsecond,
    kilometer,
+   _2, _3, _4, _5, _6, _7,
    -- * Helmert transform between geodetic reference systems
    Helmert (..),
    inverseHelmert,
@@ -69,11 +70,16 @@ arcsecond = arcminute / 60
 kilometer :: Double
 kilometer = 1000
 
-
--- | Small integers, specialised to @Int@, used for raising to powers.
---
--- If you say @x^2@ then Haskell complains that the @2@ has ambiguous type, so you
--- need to say @x^(2::Int)@ to disambiguate it. This gets tedious in complex formulae.
+-- | Lots of geodetic calculations involve integer powers. Writing e.g. @x ^ 2@ causes
+-- GHC to complain that the @2@ has ambiguous type. @x ** 2@ doesn't complain
+-- but is much slower. So for convenience, here are small integers with type @Int@.
+_2, _3, _4, _5, _6, _7 :: Int
+_2 = 2
+_3 = 3
+_4 = 4
+_5 = 5
+_6 = 6
+_7 = 7
 
 -- | 3d vector as @(X,Y,Z)@.
 type Vec3 a = (a,a,a)
@@ -130,7 +136,8 @@ cross3 :: (Num a) => Vec3 a -> Vec3 a -> Vec3 a
 cross3 (x1,y1,z1) (x2,y2,z2) = (y1*z2 - z1*y2, z1*x2 - x1*z2, x1*y2 - y1*x2)
 
 
--- | The 7 parameter Helmert transformation. The monoid instance allows composition.
+-- | The 7 parameter Helmert transformation. The monoid instance allows composition but
+-- is only accurate for the small values used in practical ellipsoids.
 data Helmert = Helmert {
    cX, cY, cZ :: Double,  -- ^ Offset in meters
    helmertScale :: Double,  -- ^ Parts per million
@@ -246,18 +253,18 @@ minorRadius e = majorRadius e * (1 - flattening e)
 
 -- | The eccentricity squared of an ellipsoid.
 eccentricity2 :: (Ellipsoid e) => e -> Double
-eccentricity2 e = 2 * f - (f * f) where f = flattening e
+eccentricity2 e = 2 * f - f^_2 where f = flattening e
 
 -- | The second eccentricity squared of an ellipsoid.
 eccentricity'2 :: (Ellipsoid e) => e -> Double
-eccentricity'2 e = (f * (2 - f)) / (1 - f * f) where f = flattening e
+eccentricity'2 e = (f * (2 - f)) / (1 - f^_2) where f = flattening e
 
 
 -- | Distance in meters from the surface at the specified latitude to the
 -- axis of the Earth straight down. Also known as the radius of
 -- curvature in the prime vertical, and often denoted @N@.
 normal :: (Ellipsoid e) => e -> Double -> Double
-normal e lat = majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ (2 :: Int))
+normal e lat = majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ _2)
 
 
 -- | Radius of the circle of latitude: the distance from a point
@@ -271,13 +278,13 @@ latitudeRadius e lat = normal e lat * cos lat
 meridianRadius :: (Ellipsoid e) => e -> Double -> Double
 meridianRadius e lat =
    majorRadius e * (1 - eccentricity2 e)
-   / sqrt ((1 - eccentricity2 e * sin lat ** 2) ** 3)
+   / sqrt ((1 - eccentricity2 e * sin lat ^ _2) ^ _3)
 
 
 -- | Radius of curvature of the ellipsoid perpendicular to the meridian at the specified latitude, in meters.
 primeVerticalRadius :: (Ellipsoid e) => e -> Double -> Double
 primeVerticalRadius e lat =
-   majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ (2 :: Int))
+   majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ _2)
 
 
 -- | The isometric latitude. The isometric latitude is conventionally denoted by ψ

src/Geodetics/Geodetic.hs

@@ -103,7 +103,7 @@ showAngle a
    where
       sgn = if a < 0 then "-" else ""
       centisecs :: Integer
-      centisecs = abs $ round $ (a * degree * 360000)  -- hundredths of arcsec per degree.
+      centisecs = abs $ round $ (a / (arcsecond / 100))
       (d, m1) = centisecs `divMod` 360000
       (m, s1) = m1 `divMod` 6000   -- hundredths of arcsec per arcmin
       (s, ds) = s1 `divMod` 100
@@ -153,7 +153,7 @@ geoToEarth geo = (
 -- Journal of Geodesy Volume 76, Number 8 (2002), 451-454. Result is in the form
 -- @(latitude, longitude, altitude)@.
 earthToGeo :: (Ellipsoid e) => e -> ECEF -> (Double, Double, Double)
-earthToGeo e (x,y,z) = (phi, atan2 y x, sqrt (l ** 2 + p2) - norm)
+earthToGeo e (x,y,z) = (phi, atan2 y x, sqrt (l ^ _2 + p2) - norm)
    where
       -- Naming: numeric suffix inicates power. Hence x2 = x * x, x3 = x2 * x, etc.
       p2 = x * x + y * y
@@ -203,7 +203,7 @@ geometricalDistance g1 g2 = sqrt $ geometricalDistanceSq g1 g2
 
 -- | The square of the absolute distance.
 geometricalDistanceSq :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
-geometricalDistanceSq g1 g2 = (x1-x2) ** 2 + (y1-y2) ** 2 + (z1-z2) ** 2
+geometricalDistanceSq g1 g2 = (x1-x2) ^ _2 + (y1-y2) ^ _2 + (z1-z2) ^ _2
    where
       (x1,y1,z1) = geoToEarth g1
       (x2,y2,z2) = geoToEarth g2
@@ -229,13 +229,13 @@ groundDistance p1 p2 = do
      (_, (lambda, (cos2Alpha, delta, sinDelta, cosDelta, cos2DeltaM))) <-
        listToMaybe $ dropWhile converging $ take 100 $ zip lambdas $ drop 1 lambdas
      let
-       uSq = cos2Alpha * (a**2 - b**2) / b**2
+       uSq = cos2Alpha * (a^ _2 - b^ _2) / b^ _2
        bigA = 1 + uSq/16384 * (4096 + uSq * ((-768) + uSq * ((320 - 175*uSq))))
        bigB =     uSq/1024  * (256  + uSq * ((-128) + uSq * ((74 -  47* uSq))))
        deltaDelta =
          bigB * sinDelta * (cos2DeltaM +
-                             bigB/4 * (cosDelta * (2 * cos2DeltaM**2 - 1)
-                                       - bigB/6 * cos2DeltaM * (4 * sinDelta**2 - 3)
+                             bigB/4 * (cosDelta * (2 * cos2DeltaM^ _2 - 1)
+                                       - bigB/6 * cos2DeltaM * (4 * sinDelta^ _2 - 3)
                                           * (4 * cos2DeltaM - 3)))
        s = b * bigA * (delta - deltaDelta)
        alpha1 = atan2(cosU2 * sin lambda) (cosU1 * sinU2 - sinU1 * cosU2 * cos lambda)
@@ -257,19 +257,19 @@ groundDistance p1 p2 = do
       where
         sinLambda = sin lambda
         cosLambda = cos lambda
-        sinDelta = sqrt((cosU2 * sinLambda) ** 2 +
-                        (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ** 2)
+        sinDelta = sqrt((cosU2 * sinLambda) ^ _2 +
+                        (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ _2)
         cosDelta = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
         delta = atan2 sinDelta cosDelta
         sinAlpha = if sinDelta == 0 then 0 else cosU1 * cosU2 * sinLambda / sinDelta
-        cos2Alpha = 1 - sinAlpha ** 2
+        cos2Alpha = 1 - sinAlpha ^ _2
         cos2DeltaM = if cos2Alpha == 0
                      then 0
                      else cosDelta - 2 * sinU1 * sinU2 / cos2Alpha
         c = (f/16) * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha))
         lambda1 = l + (1-c) * f * sinAlpha
                   * (delta + c * sinDelta
-                     * (cos2DeltaM + c * cosDelta *(2 * cos2DeltaM ** 2 - 1)))
+                     * (cos2DeltaM + c * cosDelta *(2 * cos2DeltaM ^ _2 - 1)))
     lambdas = iterate (nextLambda . fst) (l, undefined)
     converging ((l1,_),(l2,_)) = abs (l1 - l2) > 1e-14
 

src/Geodetics/Grid.hs

@@ -23,6 +23,7 @@ module Geodetics.Grid (
 
 import Data.Char
 import Geodetics.Altitude
+import Geodetics.Ellipsoids
 import Geodetics.Geodetic
 
 -- | A Grid is a two-dimensional projection of the ellipsoid onto a plane. Any given type of grid can
@@ -111,7 +112,7 @@ offsetDistance = sqrt . offsetDistanceSq
 -- | The square of the distance represented by an offset.
 offsetDistanceSq :: GridOffset -> Double
 offsetDistanceSq off =
-   deltaEast off ** 2 + deltaNorth off ** 2 + deltaAltitude off ** 2
+   deltaEast off ^ _2 + deltaNorth off ^ _2 + deltaAltitude off ^ _2
 
 
 -- | The direction represented by an offset, as bearing to the right of North.
@@ -142,19 +143,20 @@ unsafeGridCoerce base p = GridPoint (eastings p) (northings p) (altitude p) base
 -- in units of one tenth of the grid square, the second one hundredth, and so on.
 -- The first result is the lower limit of the result, and the second is the size
 -- of the specified offset.
--- So for instance @fromGridDigits (100 *~ kilo meter) "237"@ will return
+-- So for instance @fromGridDigits (100 * kilometer) "237"@ will return
 --
 -- > Just (23700 meters, 100 meters)
 --
 -- If there are any non-digits in the string then the function returns @Nothing@.
 fromGridDigits :: Double -> String -> Maybe (Double, Double)
 fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing
    where
-      n = length ds
+      n :: Integer
+      n = fromIntegral $ length ds
       d = sum $ zipWith (*)
          (map (fromIntegral . digitToInt) ds)
          (drop 1 $ iterate (/ 10) sq)
-      p = sq / (10 ** (fromIntegral n))
+      p = sq / fromIntegral (10 ^ n)
 
 -- | Convert a distance into a digit string suitable for printing as part
 -- of a grid reference. The result is the nearest position to the specified

src/Geodetics/Stereographic.hs

@@ -56,7 +56,7 @@ mkGridStereo tangent origin scale = GridStereo {
       e2 = eccentricity2 ellipse
       e = sqrt e2
       r = sqrt $ meridianRadius ellipse lat0 * primeVerticalRadius ellipse lat0
-      n = sqrt $ 1 + ((e2 * cos lat0 ** 4)/(1 - e2))
+      n = sqrt $ 1 + ((e2 * cos lat0 ^ _4)/(1 - e2))
       s1 = (1 + sinLat0) / (1 - sinLat0)
       s2 = (1 - e * sinLat0) / (1 + e * sinLat0)
       w1 = (s1 * s2 ** e) ** n
@@ -110,7 +110,7 @@ instance (Ellipsoid e) => GridClass (GridStereo e) e where
          long = (longC - long0) / gridN grid + long0
          isoLat = log ((1 + sinLatC) / (gridC grid * (1 - sinLatC))) / (2 * gridN grid)
          lat1 = 2 * atan (exp isoLat) - pi/2
-         next lat = lat - (isoN - isoLat) * cos lat * (1 - e2 * sin lat ** 2) / (1 - e2)
+         next lat = lat - (isoN - isoLat) * cos lat * (1 - e2 * sin lat ^ _2) / (1 - e2)
             where isoN = isometricLatitude (gridEllipsoid grid) lat
                   e2 = eccentricity2 $ gridEllipsoid grid
          lats = Stream.iterate next lat1

src/Geodetics/TransverseMercator.hs

@@ -43,10 +43,10 @@ mkGridTM origin offset sf =
    GridTM {trueOrigin = origin,
            falseOrigin = offset,
            gridScale = sf,
-           gridN1 = 1 + n + (5/4) * n**2 + (5/4) * n**3,
-           gridN2 = 3 * n + 3 * n**2 + (21/8) * n**3,
-           gridN3 = (15/8) * (n**2 + n**3),
-           gridN4 = (35/24) * n**3
+           gridN1 = 1 + n + (5/4) * n^ _2 + (5/4) * n^ _3,
+           gridN2 = 3 * n + 3 * n^ _2 + (21/8) * n^ _3,
+           gridN3 = (15/8) * (n^ _2 + n^ _3),
+           gridN4 = (35/24) * n^ _3
         }
     where
        f = flattening $ ellipsoid origin
@@ -66,75 +66,90 @@ m grid lat = bF0 * (gridN1 grid * dLat
 
 
 instance (Ellipsoid e) => GridClass (GridTM e) e where
-   fromGrid p = Geodetic
-      (lat' - east' ** 2 * tanLat / (2 * rho * v)  -- Term VII
-            + east' ** 4 * (tanLat / (24 * rho * v ** 3))
-                           * (5 + 3 * tanLat ** 2 + eta2 - 9 * tanLat ** 2 * eta2)  -- Term VIII
-            - east' ** 6 * (tanLat / (720 * rho * v ** 5))
-                           * (61 + 90 * tanLat ** 2 + 45 * tanLat ** 4)) -- Term IX
-      (longitude (trueOrigin grid)
-            + east' / (cosLat * v)  -- Term X
-            - (east' ** 3 / (6 * cosLat * v ** 3)) * (v / rho + 2 * tanLat ** 2)  -- Term XI
-            + (east' ** 5 / (120 * cosLat * v ** 5))
-                 * (5 + 28 * tanLat ** 2  + 24 * tanLat ** 4)  -- Term XII
-            - (east' ** 5 * east' ** 2 / (5040 * cosLat * v ** 7))
-                 * (61 + 662 * tanLat ** 2 + 1320 * tanLat ** 4 + 720 * tanLat ** 6)) -- Term XIIa
-     0
-     (gridEllipsoid grid)
-
+   fromGrid p = -- trace traceMsg $
+      Geodetic
+         (lat' - east' ^ _2 * term_VII + east' ^ _4 * term_VIII - east' ^ _6 * term_IX)
+         (longitude (trueOrigin grid)
+               + east' * term_X - east' ^ _3 * term_XI + east' ^ _5 * term_XII - east' ^ _7 * term_XIIa)
+         (altGP p)
+         (gridEllipsoid grid)
       where
          GridPoint east' north' _ _ = falseOrigin grid `applyOffset` p
          lat' = fst $ Stream.head $ Stream.dropWhile ((> 1e-5) . abs . snd)
                $ Stream.tail $ Stream.iterate next (latitude $ trueOrigin grid, 1)
             where
                next (phi, _) = let delta = north' - m grid phi in (phi + delta / aF0, delta)
-
+         -- Terms defined in [1]
+         term_VII  = tanLat / (2 * rho * v)
+         term_VIII = (tanLat / (24 * rho * v ^ _3))  * (5 + 3 * tanLat ^ _2 + eta2 - 9 * tanLat ^ _2 * eta2)
+         term_IX   = (tanLat / (720 * rho * v ^ _5)) * (61 + 90 * tanLat ^ _2 + 45 * tanLat ^ _4)
+         term_X    = 1                                                                 / (cosLat * v)
+         term_XI   = (v / rho + 2 * tanLat ^ _2)                                       / (6 * cosLat * v ^ _3)
+         term_XII  = ( 5 +  28 * tanLat ^ _2 +   24 * tanLat ^ _4)                     / (120 * cosLat * v ^ _5)
+         term_XIIa = (61 + 662 * tanLat ^ _2 + 1320 * tanLat ^ _4 + 720 * tanLat ^ _6) / (5040 * cosLat * v ^ _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 = sin lat'
          cosLat = cos lat'
          tanLat = tan lat'
          sinLat2 = sinLat * sinLat
          v = aF0 / sqrt (1 - e2 * sinLat2)
-         rho = aF0 * (1 - e2) * (1 - e2 * sinLat2) ** (-1.5)
+         rho = v * (1 - e2) / (1 - e2 * sinLat2)
          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
+   toGrid grid geo = -- trace traceMsg $ 
+      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)
+         rho = v * (1 - e2) / (1 - e2 * sinLat2)
          eta2 = v / rho - 1
          off = GridOffset
                   (dLong * term_IV
-                   + dLong ** 3 * term_V
-                   + dLong ** 5 * term_VI)
-                  (m grid lat + dLong ** 2 * term_II
-                     + dLong ** 4 * term_III
-                     + dLong ** 6 * term_IIIa)
+                   + dLong ^ _3 * term_V
+                   + dLong ^ _5 * term_VI)
+                  (m grid lat + dLong ^ _2 * term_II
+                     + dLong ^ _4 * term_III
+                     + dLong ^ _6 * term_IIIa)
                   0
          -- Terms defined in [1].
          term_II   = (v/2) * sinLat * cosLat
-         term_III  = (v/24) * sinLat * cosLat ** 3
-                     * (5 - tanLat ** 2 + 9 * eta2)
-         term_IIIa = (v/720) * sinLat * cosLat ** 5
-                     * (61 - 58 * tanLat ** 2 + tanLat ** 4)
+         term_III  = (v/24) * sinLat * cosLat ^ _3
+                     * (5 - tanLat ^ _2 + 9 * eta2)
+         term_IIIa = (v/720) * sinLat * cosLat ^ _5
+                     * (61 - 58 * tanLat ^ _2 + tanLat ^ _4)
          term_IV   = v * cosLat
-         term_V    = (v/6) * cosLat ** 3 * (v/rho - tanLat ** 2)
-         term_VI   = (v/120) * cosLat ** 5
-                     * (5 - 18 * tanLat ** 2
-                              + tanLat ** 4 + 14 * eta2
-                              - 58 * tanLat ** 2 * eta2)
+         term_V    = (v/6) * cosLat ^ _3 * (v/rho - tanLat ^ _2)
+         term_VI   = (v/120) * cosLat ^ _5
+                     * (5 - 18 * tanLat ^ _2
+                              + tanLat ^ _4 + 14 * eta2
+                              - 58 * tanLat ^ _2 * eta2)
+
+         -- Trace message for debugging. Uncomment this code to inspect intermediate values.
          {-
-         -- 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",
+            "I    = ", show $ m grid lat - deltaNorth (falseOrigin grid), "\n",  -- 
             "II   = ", show term_II, "\n",
             "III  = ", show term_III, "\n",
             "IIIa = ", show term_IIIa, "\n",