diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -20,3 +20,6 @@
 
 Version 0.1.2: Fixed bugs #16 and #17: Unicode PRIME and DOUBLE PRIME now allowed in
    position strings, and the degree symbol is allowed for decimal degrees.
+
+Version 1.0.0: Removed dependency on Dimensional library. This is a breaking change:
+   hence the major version bump. Also fixed bug #18 (and #19).
diff --git a/geodetics.cabal b/geodetics.cabal
--- a/geodetics.cabal
+++ b/geodetics.cabal
@@ -1,20 +1,20 @@
+cabal-version:  3.0
 name:           geodetics
-version:        0.1.2
-cabal-version:  >= 1.10
+version:        1.0.0
 build-type:     Simple
 author:         Paul Johnson <paul@cogito.org.uk>
-data-files:
+extra-doc-files:
                 AddingProjections.txt,
                 LICENSE,
                 README.md,
                 changelog.md,
                 ToDo.txt
-license:        BSD3
-copyright:      Paul Johnson 2018.
+license:        BSD-3-Clause
+copyright:      Paul Johnson 2018,2024
 synopsis:       Terrestrial coordinate systems and geodetic calculations.
 description:    Precise geographical coordinates (latitude & longitude), with conversion between
                 different reference frames and projections.
-                .
+                
                 Certain distinguished reference frames and grids are given distinct
                 types so that coordinates expressed within them cannot be confused with
                 from coordinates in other frames.
@@ -22,7 +22,7 @@
 maintainer:     Paul Johnson <paul@cogito.org.uk>
 homepage:       https://github.com/PaulJohnson/geodetics
 category:       Geography
-tested-with:    GHC==8.6.3
+tested-with:    GHC==9.10.1
 
 source-repository head
   type:     git
@@ -31,10 +31,9 @@
 library
   hs-source-dirs:  src
   build-depends:
-                   base >= 4.6 && < 5,
-                   dimensional >= 1.3,
-                   array >= 0.4,
-                   semigroups >= 0.9
+                   base >= 4.17 && < 5,
+                   array >= 0.1 && < 0.6,
+                   Stream >= 0.4.6 && < 0.5
   ghc-options:     -Wall
   exposed-modules:
                    Geodetics.Altitude,
@@ -55,12 +54,11 @@
   build-depends:   geodetics,
                    base >= 4.6 && < 5,
                    HUnit >= 1.2,
-                   dimensional >= 1.3,
                    QuickCheck >= 2.4,
                    test-framework >= 0.4.1,
                    test-framework-quickcheck2,
                    test-framework-hunit,
-                   array >= 0.4,
+                   array,
                    checkers
   hs-source-dirs:
                    test
diff --git a/src/Geodetics/Altitude.hs b/src/Geodetics/Altitude.hs
--- a/src/Geodetics/Altitude.hs
+++ b/src/Geodetics/Altitude.hs
@@ -2,16 +2,15 @@
   HasAltitude (..)
 ) where
 
-import Numeric.Units.Dimensional.Prelude
 
--- | All geographical coordinate systems need the concept of#
+-- | All geographical coordinate systems need the concept of
 -- altitude above a reference point, usually associated with
 -- local sea level.
 -- 
 -- Minimum definition: altitude, setAltitude.
 class HasAltitude a where
-   altitude :: a -> Length Double
-   setAltitude :: Length Double -> a -> a
+   altitude :: a -> Double
+   setAltitude :: Double -> a -> a
    -- | Set altitude to zero.
    groundPosition :: a -> a
-   groundPosition = setAltitude _0
+   groundPosition = setAltitude 0
diff --git a/src/Geodetics/Ellipsoids.hs b/src/Geodetics/Ellipsoids.hs
--- a/src/Geodetics/Ellipsoids.hs
+++ b/src/Geodetics/Ellipsoids.hs
@@ -1,19 +1,11 @@
 {-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE DataKinds #-}
-{-# LANGUAGE DeriveDataTypeable #-}
-{-# LANGUAGE DeriveGeneric #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE KindSignatures #-}
-{-# LANGUAGE NoImplicitPrelude #-}
-{-# LANGUAGE PatternGuards #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE RoleAnnotations #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE StandaloneDeriving #-}
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeOperators #-}
 
 {- | An Ellipsoid is a reasonable best fit for the surface of the
 Earth over some defined area. WGS84 is the standard used for the whole
@@ -22,12 +14,18 @@
 -}
 
 module Geodetics.Ellipsoids (
-   -- ** Helmert transform between geodetic reference systems
+   -- * Useful constants
+   degree,
+   arcminute,
+   arcsecond,
+   kilometer,
+   _2, _3, _4, _5, _6, _7,
+   -- * Helmert transform between geodetic reference systems
    Helmert (..),
    inverseHelmert,
    ECEF,
    applyHelmert,
-   -- ** Ellipsoid models of the Geoid
+   -- * Ellipsoid models of the Geoid
    Ellipsoid (..),
    WGS84 (..),
    LocalEllipsoid (..),
@@ -35,13 +33,13 @@
    minorRadius,
    eccentricity2,
    eccentricity'2,
-   -- ** Auxiliary latitudes and related Values
+   -- * Auxiliary latitudes and related Values
    normal,
    latitudeRadius,
    meridianRadius,
    primeVerticalRadius,
    isometricLatitude,
-   -- ** Tiny linear algebra library for 3D vectors
+   -- * Tiny linear algebra library for 3D vectors
    Vec3,
    Matrix3,
    add3,
@@ -54,14 +52,35 @@
    cross3
 ) where
 
-import Data.Monoid (Monoid)
-import Data.Semigroup (Semigroup, (<>))
-import Numeric.Units.Dimensional
-import Numeric.Units.Dimensional.Prelude
-import qualified Numeric.Units.Dimensional.Dimensions.TypeLevel as T
--- import Prelude ()  -- Numeric instances.
 
+-- | All angles in this library are in radians. This is one degree in radians.
+degree :: Double
+degree = pi/180
 
+-- | One arc-minute in radians.
+arcminute :: Double
+arcminute = degree / 60
+
+-- | One arc-second in radians.
+arcsecond :: Double
+arcsecond = arcminute / 60
+
+
+-- | All distances in this library are in meters. This is one kilometer in meters.
+kilometer :: Double
+kilometer = 1000
+
+-- | 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)
 
@@ -70,31 +89,29 @@
 
 
 -- | Multiply a vector by a scalar.
-scale3 :: (Num a) =>
-   Vec3 (Quantity d a) -> Quantity d' a -> Vec3 (Quantity (d T.* d') a)
+scale3 :: (Num a) =>  Vec3 a -> a -> Vec3 a
 scale3 (x,y,z) s = (x*s, y*s, z*s)
 
 
 -- | Negation of a vector.
-negate3 :: (Num a) => Vec3 (Quantity d a) -> Vec3 (Quantity d a)
+negate3 :: (Num a) => Vec3 a -> Vec3 a
 negate3 (x,y,z) = (negate x, negate y, negate z)
 
 -- | Add two vectors
-add3 :: (Num a) => Vec3 (Quantity d a) -> Vec3 (Quantity d a) -> Vec3 (Quantity d a)
+add3 :: (Num a) => Vec3 a -> Vec3 a -> Vec3 a
 add3 (x1,y1,z1) (x2,y2,z2) = (x1+x2, y1+y2, z1+z2)
 
 
 -- | Multiply a matrix by a vector in the Dimensional type system.
 transform3 :: (Num a) =>
-   Matrix3 (Quantity d a) -> Vec3 (Quantity d' a) -> Vec3 (Quantity (d T.* d') a)
+   Matrix3 a -> Vec3 a -> Vec3 a
 transform3 (tx,ty,tz) v = (t tx v, t ty v, t tz v)
    where
       t (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2
 
 
 -- | Inverse of a 3x3 matrix.
-invert3 :: (Fractional a) =>
-   Matrix3 (Quantity d a) -> Matrix3 (Quantity ((d T.* d)/(d T.* d T.* d)) a)
+invert3 :: (Fractional a) => Matrix3 a -> Matrix3 a
 invert3 ((x1,y1,z1),
          (x2,y2,z2),
          (x3,y3,z3)) =
@@ -111,21 +128,21 @@
 
 
 -- | Dot product of two vectors
-dot3 :: (Num a) =>
-   Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Quantity (d1 T.* d2) a
+dot3 :: (Num a) => Vec3 a -> Vec3 a -> a
 dot3 (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2
 
 -- | Cross product of two vectors
-cross3 :: (Num a) =>
-   Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Vec3 (Quantity (d1 T.* d2) a)
+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 :: Length Double,
-   helmertScale :: Dimensionless Double,  -- ^ Parts per million
-   rX, rY, rZ :: Dimensionless Double } deriving (Eq, Show)
+   cX, cY, cZ :: Double,  -- ^ Offset in meters
+   helmertScale :: Double,  -- ^ Parts per million
+   rX, rY, rZ :: Double  -- ^ Rotation around each axis in radians.
+} deriving (Eq, Show)
 
 instance Semigroup Helmert where
     h1 <> h2 = Helmert (cX h1 + cX h2) (cY h1 + cY h2) (cZ h1 + cZ h2)
@@ -133,7 +150,7 @@
                        (rX h1 + rX h2) (rY h1 + rY h2) (rZ h1 + rZ h2)
 
 instance Monoid Helmert where
-   mempty = Helmert (0 *~ meter) (0 *~ meter) (0 *~ meter) _0 _0 _0 _0
+   mempty = Helmert 0 0 0 0 0 0 0
    mappend = (<>)
 
 -- | The inverse of a Helmert transformation.
@@ -145,7 +162,7 @@
 
 -- | Earth-centred, Earth-fixed coordinates as a vector. The origin and axes are
 -- not defined: use with caution.
-type ECEF = Vec3 (Length Double)
+type ECEF = Vec3 Double
 
 -- | Apply a Helmert transformation to earth-centered coordinates.
 applyHelmert:: Helmert -> ECEF -> ECEF
@@ -154,7 +171,7 @@
       cY h + s * (        rZ h  * x +        y - rX h * z),
       cZ h + s * (negate (rY h) * x + rX h * y +        z))
    where
-      s = _1 + helmertScale h * (1e-6 *~ one)
+      s = 1 + helmertScale h * 1e-6
 
 
 -- | An Ellipsoid is defined by the major radius and the inverse flattening (which define its shape),
@@ -170,17 +187,17 @@
 -- > helmertToWGS84 = applyHelmert . helmert
 -- > helmertFromWGS84 e . helmertToWGS84 e = id
 class (Show a, Eq a) => Ellipsoid a where
-   majorRadius :: a -> Length Double
-   flatR :: a -> Dimensionless Double
+   majorRadius :: a -> Double
+   flatR :: a -> Double
       -- ^ Inverse of the flattening.
-   helmert :: a -> Helmert
-   helmertToWSG84 :: a -> ECEF -> ECEF
+   helmert :: a -> Helmert  -- ^ The Helmert parameters relative to WGS84,
+   helmertToWGS84 :: a -> ECEF -> ECEF
       -- ^ The Helmert transform that will convert a position wrt
       -- this ellipsoid into a position wrt WGS84.
-   helmertToWSG84 e = applyHelmert (helmert e)
-   helmertFromWSG84 :: a -> ECEF -> ECEF
+   helmertToWGS84 e = applyHelmert (helmert e)
+   helmertFromWGS84 :: a -> ECEF -> ECEF
       -- ^ And its inverse.
-   helmertFromWSG84 e = applyHelmert (inverseHelmert $ helmert e)
+   helmertFromWGS84 e = applyHelmert (inverseHelmert $ helmert e)
 
 
 -- | The WGS84 geoid, major radius 6378137.0 meters, flattening = 1 / 298.257223563
@@ -197,11 +214,11 @@
    show _ = "WGS84"
 
 instance Ellipsoid WGS84 where
-   majorRadius _ = 6378137.0 *~ meter
-   flatR _ = 298.257223563 *~ one
+   majorRadius _ = 6378137.0
+   flatR _ = 298.257223563
    helmert _ = mempty
-   helmertToWSG84 _ = id
-   helmertFromWSG84 _ = id
+   helmertToWGS84 _ = id
+   helmertFromWGS84 _ = id
 
 
 -- | Ellipsoids other than WGS84, used within a defined geographical area where
@@ -211,9 +228,10 @@
 -- Creating two different local ellipsoids with the same name is a Bad Thing.
 data LocalEllipsoid = LocalEllipsoid {
    nameLocal :: String,
-   majorRadiusLocal :: Length Double,
-   flatRLocal :: Dimensionless Double,
-   helmertLocal :: Helmert } deriving (Eq)
+   majorRadiusLocal :: Double,
+   flatRLocal :: Double,
+   helmertLocal :: Helmert
+} deriving (Eq)
 
 instance Show LocalEllipsoid where
     show = nameLocal
@@ -225,48 +243,48 @@
 
 
 -- | Flattening (f) of an ellipsoid.
-flattening :: (Ellipsoid e) => e -> Dimensionless Double
-flattening e = _1 / flatR e
+flattening :: (Ellipsoid e) => e -> Double
+flattening e = 1 / flatR e
 
--- | The minor radius of an ellipsoid.
-minorRadius :: (Ellipsoid e) => e -> Length Double
-minorRadius e = majorRadius e * (_1 - flattening e)
+-- | The minor radius of an ellipsoid in meters.
+minorRadius :: (Ellipsoid e) => e -> Double
+minorRadius e = majorRadius e * (1 - flattening e)
 
 
 -- | The eccentricity squared of an ellipsoid.
-eccentricity2 :: (Ellipsoid e) => e -> Dimensionless Double
-eccentricity2 e = _2 * f - (f * f) where f = flattening e
+eccentricity2 :: (Ellipsoid e) => e -> Double
+eccentricity2 e = 2 * f - f^_2 where f = flattening e
 
 -- | The second eccentricity squared of an ellipsoid.
-eccentricity'2 :: (Ellipsoid e) => e -> Dimensionless Double
-eccentricity'2 e = (f * (_2 - f)) / (_1 - f * f) where f = flattening e
+eccentricity'2 :: (Ellipsoid e) => e -> Double
+eccentricity'2 e = (f * (2 - f)) / (1 - f^_2) where f = flattening e
 
 
--- | Distance from the surface at the specified latitude to the
+-- | 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 -> Angle Double -> Length Double
-normal e lat = majorRadius e / sqrt (_1 - eccentricity2 e * sin lat ^ pos2)
+normal :: (Ellipsoid e) => e -> Double -> Double
+normal e lat = majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ _2)
 
 
 -- | Radius of the circle of latitude: the distance from a point
--- at that latitude to the axis of the Earth.
-latitudeRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double
+-- at that latitude to the axis of the Earth, in meters.
+latitudeRadius :: (Ellipsoid e) => e -> Double -> Double
 latitudeRadius e lat = normal e lat * cos lat
 
 
--- | Radius of curvature in the meridian at the specified latitude.
+-- | Radius of curvature in the meridian at the specified latitude, in meters
 -- Often denoted @M@.
-meridianRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double
+meridianRadius :: (Ellipsoid e) => e -> Double -> Double
 meridianRadius e lat =
-   majorRadius e * (_1 - eccentricity2 e)
-   / sqrt ((_1 - eccentricity2 e * sin lat ^ pos2) ^ pos3)
+   majorRadius e * (1 - eccentricity2 e)
+   / sqrt ((1 - eccentricity2 e * sin lat ^ _2) ^ _3)
 
 
--- | Radius of curvature of the ellipsoid perpendicular to the meridian at the specified latitude.
-primeVerticalRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double
+-- | 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 ^ pos2)
+   majorRadius e / sqrt (1 - eccentricity2 e * sin lat ^ _2)
 
 
 -- | The isometric latitude. The isometric latitude is conventionally denoted by ψ
@@ -275,7 +293,7 @@
 -- Mercator projection. The name "isometric" arises from the fact that at any point
 -- on the ellipsoid equal increments of ψ and longitude λ give rise to equal distance
 -- displacements along the meridians and parallels respectively.
-isometricLatitude :: (Ellipsoid e) => e -> Angle Double -> Angle Double
+isometricLatitude :: (Ellipsoid e) => e -> Double -> Double
 isometricLatitude ellipse lat = atanh sinLat - e * atanh (e * sinLat)
    where
       sinLat = sin lat
diff --git a/src/Geodetics/Geodetic.hs b/src/Geodetics/Geodetic.hs
--- a/src/Geodetics/Geodetic.hs
+++ b/src/Geodetics/Geodetic.hs
@@ -1,5 +1,5 @@
 module Geodetics.Geodetic (
-   -- ** Geodetic Coordinates
+   -- * Geodetic Coordinates
    Geodetic (..),
    readGroundPosition,
    toLocal,
@@ -10,25 +10,21 @@
    groundDistance,
    properAngle,
    showAngle,
-   -- ** Earth Centred Earth Fixed Coordinates
+   -- * Earth Centred Earth Fixed Coordinates
    ECEF,
    geoToEarth,
    earthToGeo,
-   -- ** Re-exported for convenience
+   -- * Re-exported for convenience
    WGS84 (..)
 ) where
 
 
 import Data.Char (chr)
-import Data.Function
 import Data.Maybe
-import Data.Monoid
 import Geodetics.Altitude
 import Geodetics.Ellipsoids
 import Geodetics.LatLongParser
-import Numeric.Units.Dimensional.Prelude hiding ((.))
 import Text.ParserCombinators.ReadP
-import qualified Prelude as P
 
 -- | Defines a three-D position on or around the Earth using latitude,
 -- longitude and altitude with respect to a specified ellipsoid, with
@@ -60,8 +56,8 @@
 -- the same to within a given tolerance then use "geometricDistance"
 -- (or its squared variant to avoid an extra @sqrt@ operation).
 data Geodetic e = Geodetic {
-   latitude, longitude :: Angle Double,
-   geoAlt :: Length Double,
+   latitude, longitude :: Double,  -- ^ In radians.
+   geoAlt :: Double,  -- ^ In meters.
    ellipsoid :: e
 }
 
@@ -70,8 +66,7 @@
       showAngle (abs $ latitude g),  " ", letter "SN" (latitude g),  ", ",
       showAngle (abs $ longitude g), " ", letter "WE" (longitude g), ", ",
       show (altitude g), " ", show (ellipsoid g)]
-      where letter s n = [s !! (if n < _0 then 0 else 1)]
-
+      where letter s n = [s !! (if n < 0 then 0 else 1)]
 
 
 -- | Read the latitude and longitude of a ground position and
@@ -94,25 +89,24 @@
 readGroundPosition e str =
    case map fst $ filter (null . snd) $ readP_to_S latLong str of
       [] -> Nothing
-      (lat,long) : _ -> Just $ groundPosition $ Geodetic (lat *~ degree) (long *~ degree) undefined e
+      (lat,long) : _ -> Just $ groundPosition $ Geodetic (lat * degree) (long * degree) undefined e
 
 
 -- | Show an angle as degrees, minutes and seconds to two decimal places.
-showAngle :: Angle Double -> String
+showAngle :: Double -> String
 showAngle a
-   | isNaN a1       = "NaN"  -- Not a Nangle
-   | isInfinite a1  = sgn ++ "Infinity"
+   | isNaN a        = "NaN"  -- Not a Nangle
+   | isInfinite a   = sgn ++ "Infinity"
    | otherwise      = concat [sgn, show d, [chr 0xB0, ' '],
                               show m, "\8242 ",
                               show s, ".", dstr, "\8243" ]
    where
-      a1 = a /~ one
-      sgn = if a < _0 then "-" else ""
+      sgn = if a < 0 then "-" else ""
       centisecs :: Integer
-      centisecs = P.abs $ P.round $ (a /~ degree) P.* 360000  -- hundredths of arcsec per degree.
-      (d, m1) = centisecs `P.divMod` 360000
-      (m, s1) = m1 `P.divMod` 6000   -- hundredths of arcsec per arcmin
-      (s, ds) = s1 `P.divMod` 100
+      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
       dstr = reverse $ take 2 $ reverse (show ds) ++ "00" -- Decimal fraction with zero padding.
 
 
@@ -128,8 +122,8 @@
 antipode g = Geodetic lat long (geoAlt g) (ellipsoid g)
    where
       lat = negate $ latitude g
-      long' = longitude g - 180 *~ degree
-      long | long' < _0  = long' + 360 *~ degree
+      long' = longitude g - 180 * degree
+      long | long' < 0  = long' + 360 * degree
            | otherwise  = long'
 
 
@@ -140,7 +134,7 @@
 geoToEarth geo = (
       (n + h) * coslat * coslong,
       (n + h) * coslat * sinlong,
-      (n * (_1 - eccentricity2 e) + h) * sinlat)
+      (n * (1 - eccentricity2 e) + h) * sinlat)
    where
       n = normal e $ latitude geo
       e = ellipsoid geo
@@ -156,26 +150,27 @@
 --
 -- Uses the closed form solution of H. Vermeille: Direct
 -- transformation from geocentric coordinates to geodetic coordinates.
--- Journal of Geodesy Volume 76, Number 8 (2002), 451-454
-earthToGeo :: (Ellipsoid e) => e -> ECEF -> (Angle Double, Angle Double, Length Double)
-earthToGeo e (x,y,z) = (phi, atan2 y x, sqrt (l ^ pos2 + p2) - norm)
+-- 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)
    where
       -- Naming: numeric suffix inicates power. Hence x2 = x * x, x3 = x2 * x, etc.
-      p2 = x ^ pos2 + y ^ pos2
+      p2 = x * x + y * y
       a = majorRadius e
-      a2 = a ^ pos2
+      a2 = a * a
       e2 = eccentricity2 e
-      e4 = e2 ^ pos2
-      zeta = (_1-e2) * (z ^ pos2 / a2)
-      rho = (p2 / a2 + zeta - e4) / _6
-      rho2 = rho ^ pos2
+      e4 = e2 * e2
+      zeta = (1-e2) * (z * z / a2)
+      rho = (p2 / a2 + zeta - e4) / 6
+      rho2 = rho * rho
       rho3 = rho * rho2
-      s = e4 * zeta * p2 / (_4 * a2)
-      t = cbrt (s + rho3 + sqrt (s * (s + _2 * rho3)))
+      s = e4 * zeta * p2 / (4 * a2)
+      t = (s + rho3 + sqrt (s * (s + 2 * rho3))) ** (1/3) -- Cube root
       u = rho + t + rho2 / t
-      v = sqrt (u ^ pos2 + e4 * zeta)
-      w = e2 * (u + v - zeta) / (_2 * v)
-      kappa = _1 + e2 * (sqrt (u + v + w ^ pos2) + w) / (u + v)
+      v = sqrt (u * u + e4 * zeta)
+      w = e2 * (u + v - zeta) / (2 * v)
+      kappa = 1 + e2 * (sqrt (u + v + w * w) + w) / (u + v)
       phi = atan (kappa * z / sqrt p2)
       norm = normal e phi
       l = z + e2 * norm * sin phi
@@ -203,12 +198,12 @@
 -- points. They must be on the same ellipsoid.
 -- Note that this is not the geodetic distance taken by following
 -- the curvature of the earth.
-geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Length Double
+geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
 geometricalDistance g1 g2 = sqrt $ geometricalDistanceSq g1 g2
 
--- | The square of the absolute distance. Comes out as "Area" type of course.
-geometricalDistanceSq :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Area Double
-geometricalDistanceSq g1 g2 = (x1-x2) ^ pos2 + (y1-y2) ^ pos2 + (z1-z2) ^ pos2
+-- | 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
    where
       (x1,y1,z1) = geoToEarth g1
       (x2,y2,z2) = geoToEarth g2
@@ -229,23 +224,19 @@
 -- equations\". T. Vincenty. Survey Review XXII 176, April
 -- 1975. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>
 groundDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e ->
-                  Maybe (Length Double, Dimensionless Double, Dimensionless Double)
+                  Maybe (Double, Double, Double)
 groundDistance p1 p2 = do
      (_, (lambda, (cos2Alpha, delta, sinDelta, cosDelta, cos2DeltaM))) <-
-       listToMaybe $ dropWhile converging $ take 100 $ zip lambdas $ tail lambdas
+       listToMaybe $ dropWhile converging $ take 100 $ zip lambdas $ drop 1 lambdas
      let
-       uSq = cos2Alpha * (a^pos2 - b^pos2) / b^pos2
-       bigA = _1 + uSq/(16384*~one) * ((4096*~one) + uSq *
-                                      (((-768)*~one) + uSq * ((320*~one)
-                                                            - (175*~one)*uSq)))
-       bigB = uSq/(1024*~one) * ((256*~one) +
-                                 uSq * (((-128)*~one) +
-                                        uSq * ((74*~one) - (47*~one)*uSq)))
+       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^pos2 - _1)
-                                        - bigB/_6 * cos2DeltaM * (_4 * sinDelta^pos2 - _3)
-                                          * (_4 * cos2DeltaM - _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)
        alpha2 = atan2(cosU1 * sin lambda) (cosU1 * sinU2 * cos lambda - sinU1 * cosU2)
@@ -255,8 +246,8 @@
     a = majorRadius $ ellipsoid p1
     b = minorRadius $ ellipsoid p1
     l = abs $ longitude p1 - longitude p2
-    u1 = atan ((_1-f) * tan (latitude p1))
-    u2 = atan ((_1-f) * tan (latitude p2))
+    u1 = atan ((1-f) * tan (latitude p1))
+    u2 = atan ((1-f) * tan (latitude p2))
     sinU1 = sin u1
     cosU1 = cos u1
     sinU2 = sin u2
@@ -266,25 +257,25 @@
       where
         sinLambda = sin lambda
         cosLambda = cos lambda
-        sinDelta = sqrt((cosU2 * sinLambda) ^ pos2 +
-                        (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ pos2)
+        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 ^ pos2
-        cos2DeltaM = if cos2Alpha == _0
-                     then _0
-                     else cosDelta - _2 * sinU1 * sinU2 / cos2Alpha
-        c = f/(16 *~ one) * cos2Alpha * (_4 + f * (_4 - _3 * cos2Alpha))
-        lambda1 = l + (_1-c) * f * sinAlpha
+        sinAlpha = if sinDelta == 0 then 0 else cosU1 * cosU2 * sinLambda / sinDelta
+        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 ^ pos2 - _1)))
+                     * (cos2DeltaM + c * cosDelta *(2 * cos2DeltaM ^ _2 - 1)))
     lambdas = iterate (nextLambda . fst) (l, undefined)
-    converging ((l1,_),(l2,_)) = abs (l1 - l2) > (1e-14 *~ one)
+    converging ((l1,_),(l2,_)) = abs (l1 - l2) > 1e-14
 
 
 -- | Add or subtract multiples of 2*pi so that for all @t@, @-pi < properAngle t < pi@.
-properAngle :: Angle Double -> Angle Double
+properAngle :: Double -> Double
 properAngle t
    | r1 <= negate pi    = r1 + pi2
    | r1 > pi            = r1 - pi2
@@ -292,6 +283,6 @@
    where
       pf :: Double -> (Int, Double)
       pf = properFraction  -- Shut up GHC warning about defaulting to Integer.
-      (_,r) = pf (t/pi2 /~ one)
-      r1 = (r *~ one) * pi2
-      pi2 = pi * _2
+      (_,r) = pf (t/pi2)
+      r1 = r * pi2
+      pi2 = pi * 2
diff --git a/src/Geodetics/Grid.hs b/src/Geodetics/Grid.hs
--- a/src/Geodetics/Grid.hs
+++ b/src/Geodetics/Grid.hs
@@ -1,11 +1,11 @@
-{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
+{-# LANGUAGE FunctionalDependencies #-}
 
 module Geodetics.Grid (
-   -- ** Grid types
+   -- * Grid types
    GridClass (..),
    GridPoint (..),
    GridOffset (..),
-   -- ** Grid operations
+   -- * Grid operations
    polarOffset,
    offsetScale,
    offsetNegate,
@@ -14,25 +14,21 @@
    offsetDistanceSq,
    offsetBearing,
    gridOffset,
-   -- ** Unsafe conversion
+   -- * Unsafe conversion
    unsafeGridCoerce,
-   -- ** Utility functions for grid references
+   -- * Utility functions for grid references
    fromGridDigits,
    toGridDigits
 ) where
 
 import Data.Char
-import Data.Function
-import Data.Monoid (Monoid)
-import Data.Semigroup (Semigroup, (<>))
 import Geodetics.Altitude
+import Geodetics.Ellipsoids
 import Geodetics.Geodetic
-import Numeric.Units.Dimensional.Prelude hiding ((.))
-import qualified Prelude as P
 
 -- | A Grid is a two-dimensional projection of the ellipsoid onto a plane. Any given type of grid can
 -- usually be instantiated with parameters such as a tangential point or line, and these parameters
--- will include the terrestrial reference frame ("Ellipsoid" in this library) used as a foundation. 
+-- will include the terrestrial reference frame ("Ellipsoid" in this library) used as a foundation.
 -- Hence conversion from a geodetic to a grid point requires the \"basis\" for the grid in question,
 -- and grid points carry that basis with them because without it there is no defined relationship
 -- between the grid points and terrestrial positions.
@@ -42,17 +38,17 @@
    gridEllipsoid :: r -> e
 
 
--- | A point on the specified grid. 
+-- | A point on the specified grid.
 data GridPoint r = GridPoint {
-   eastings, northings, altGP :: Length Double,
+   eastings, northings, altGP :: Double,
    gridBasis :: r
 } deriving (Show)
 
 
 instance Eq (GridPoint r) where
-   p1 == p2  = 
-      eastings p1 == eastings p2 && 
-      northings p1 == northings p2 && 
+   p1 == p2  =
+      eastings p1 == eastings p2 &&
+      northings p1 == northings p2 &&
       altGP p1 == altGP p2
 
 instance HasAltitude (GridPoint g) where
@@ -60,14 +56,13 @@
    setAltitude h gp = gp{altGP = h}
 
 
-
--- | A vector relative to a point on a grid.
+-- | A vector relative to a point on a grid. All distances are in meters.
 -- Operations that use offsets will only give
 -- meaningful results if all the points come from the same grid.
--- 
+--
 -- The monoid instance is the sum of offsets.
 data GridOffset = GridOffset {
-   deltaEast, deltaNorth, deltaAltitude :: Length Double
+   deltaEast, deltaNorth, deltaAltitude :: Double
 } deriving (Eq, Show)
 
 instance Semigroup GridOffset where
@@ -76,20 +71,20 @@
                         (deltaAltitude g1 + deltaAltitude g2)
 
 instance Monoid GridOffset where
-   mempty = GridOffset _0 _0 _0
+   mempty = GridOffset 0 0 0
    mappend = (<>)
 
--- | An offset defined by a distance and a bearing to the right of North.
+-- | An offset defined by a distance (m) and a bearing (radians) to the right of North.
 --
 -- There is no elevation parameter because we are using a plane to approximate an ellipsoid,
 -- so elevation would not provide a useful result.  If you want to work with elevations
 -- then "rayPath" will give meaningful results.
-polarOffset :: Length Double -> Angle Double -> GridOffset
-polarOffset r d = GridOffset (r * sin d) (r * cos d) _0
+polarOffset :: Double -> Double -> GridOffset
+polarOffset r d = GridOffset (r * sin d) (r * cos d) 0
 
 
 -- | Scale an offset by a scalar.
-offsetScale :: Dimensionless Double -> GridOffset -> GridOffset
+offsetScale :: Double -> GridOffset -> GridOffset
 offsetScale s off = GridOffset (deltaEast off * s)
                                (deltaNorth off * s)
                                (deltaAltitude off * s)
@@ -103,36 +98,36 @@
 
 -- Add an offset on to a point to get another point.
 applyOffset :: GridOffset -> GridPoint g -> GridPoint g
-applyOffset off p = GridPoint (eastings p + deltaEast off) 
+applyOffset off p = GridPoint (eastings p + deltaEast off)
                            (northings p + deltaNorth off)
                            (altitude p + deltaAltitude off)
                            (gridBasis p)
 
 
 -- | The distance represented by an offset.
-offsetDistance :: GridOffset -> Length Double
+offsetDistance :: GridOffset -> Double
 offsetDistance = sqrt . offsetDistanceSq
 
 
 -- | The square of the distance represented by an offset.
-offsetDistanceSq :: GridOffset -> Area Double
-offsetDistanceSq off = 
-   deltaEast off ^ pos2 + deltaNorth off ^ pos2 + deltaAltitude off ^ pos2
+offsetDistanceSq :: GridOffset -> Double
+offsetDistanceSq off =
+   deltaEast off ^ _2 + deltaNorth off ^ _2 + deltaAltitude off ^ _2
 
-              
+
 -- | The direction represented by an offset, as bearing to the right of North.
-offsetBearing :: GridOffset -> Angle Double
+offsetBearing :: GridOffset -> Double
 offsetBearing off = atan2 (deltaEast off) (deltaNorth off)
 
 
--- | The offset required to move from p1 to p2.             
+-- | The offset required to move from p1 to p2.
 gridOffset :: GridPoint g -> GridPoint g -> GridOffset
 gridOffset p1 p2 = GridOffset (eastings p2 - eastings p1)
                               (northings p2 - northings p1)
                               (altitude p2 - altitude p1)
 
 
--- | Coerce a grid point of one type into a grid point of a different type, 
+-- | Coerce a grid point of one type into a grid point of a different type,
 -- but with the same easting, northing and altitude. This is unsafe because it
 -- will produce a different position unless the two grids are actually equal.
 --
@@ -148,41 +143,41 @@
 -- 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 :: Length Double -> String -> Maybe (Length Double, Length Double)
+fromGridDigits :: Double -> String -> Maybe (Double, Double)
 fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing
    where
-      n = length ds
-      d = sum $ zipWith (*) 
-         (map ((*~ one) . fromIntegral . digitToInt) ds) 
-         (tail $ iterate (/ (10 *~ one)) sq)
-      p = sq / ((10 *~ one) ** (fromIntegral n *~ one))
-      
+      n :: Integer
+      n = fromIntegral $ length ds
+      d = sum $ zipWith (*)
+         (map (fromIntegral . digitToInt) ds)
+         (drop 1 $ iterate (/ 10) sq)
+      p = sq / fromIntegral ((10 :: Int) ^ 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
 -- number of digits, expressed as an integer count of squares and a string of digits.
 -- If any arguments are invalid then @Nothing@ is returned.
 toGridDigits ::
-   Length Double    -- ^ Size of enclosing grid square. Must be at least 1 km.
+   Double           -- ^ Size of enclosing grid square. Must be at least 1000m.
    -> Int           -- ^ Number of digits to return. Must be positive.
-   -> Length Double -- ^ Offset to convert into grid.
+   -> Double        -- ^ Offset to convert into grid (m).
    -> Maybe (Integer, String)
 toGridDigits sq n d =
-   if sq < (1 *~ kilo meter) || n < 0 || d < _0 
+   if sq < 1000 || n < 0 || d < 0
    then Nothing
    else
       Just (sqs, pad)
    where
       p :: Integer
-      p = 10 P.^ n
-      unit :: Length Double
-      unit = sq / (fromIntegral p *~ one)
-      u = round ((d / unit) /~ one)
+      p = 10 ^ n
+      unit :: Double
+      unit = sq / fromIntegral p
+      u = round (d / unit)
       (sqs, d1) = u `divMod` p
       s = show d1
-      pad = if n == 0 then "" else replicate (n P.- length s) '0' ++ s
-      
+      pad = if n == 0 then "" else replicate (n - length s) '0' ++ s
diff --git a/src/Geodetics/LatLongParser.hs b/src/Geodetics/LatLongParser.hs
--- a/src/Geodetics/LatLongParser.hs
+++ b/src/Geodetics/LatLongParser.hs
@@ -1,9 +1,10 @@
 -- | The default reader for Geodetic ground positions is flexible but slow. If you are
 -- going to read positions in a known format and performance matters then use one of
 -- the more specialised parsers here.
+--
+-- All angles are returned in degrees.
 
 module Geodetics.LatLongParser (
-
    degreesMinutesSeconds,
    degreesMinutesSecondsUnits,
    degreesDecimalMinutes,
@@ -78,7 +79,8 @@
    return $ d + ms
 
 
--- | Parse an unsigned angle written using degrees, minutes and seconds with units (° ' \"). At least one component must be specified.
+-- | Parse an unsigned angle written using degrees, minutes and seconds with units (° ' \").
+-- At least one component must be specified.
 degreesMinutesSecondsUnits :: ReadP Double
 degreesMinutesSecondsUnits = do
    (s, a) <- gather $ do
diff --git a/src/Geodetics/Path.hs b/src/Geodetics/Path.hs
--- a/src/Geodetics/Path.hs
+++ b/src/Geodetics/Path.hs
@@ -6,12 +6,10 @@
 import Control.Monad
 import Geodetics.Ellipsoids
 import Geodetics.Geodetic
-import Numeric.Units.Dimensional.Prelude
-import Prelude ()
 
 
--- | Lower and upper exclusive bounds within which a path is valid. 
-type PathValidity = (Length Double, Length Double)
+-- | Lower and upper exclusive distance bounds within which a path is valid. 
+type PathValidity = (Double, Double)
 
 -- | A path is a parametric function of distance along the path. The result is the
 -- position, and the direction of the path at that point as heading and elevation angles.
@@ -24,18 +22,19 @@
 -- Outside its validity the path function may
 -- return anything or bottom.
 data Path e = Path {
-      pathFunc :: Length Double -> (Geodetic e, Angle Double, Angle Double),
+      pathFunc :: Double -> (Geodetic e, Double, Double),
+         -- ^ Takes a length and returns a position, and direction as heading and elevation angles.
       pathValidity :: PathValidity
    }
    
 -- | Convenience value for paths that are valid for all distances.
 alwaysValid :: PathValidity
 alwaysValid = (negate inf, inf) where
-   inf = (1.0 *~ meter) / (0 *~ one)  -- Assumes IEE arithmetic.
+   inf = 1.0 / 0  -- Assumes IEE arithmetic.
 
 
 -- | True if the path is valid at that distance.
-pathValidAt :: Path e -> Length Double -> Bool
+pathValidAt :: Path e -> Double -> Bool
 pathValidAt path d = d > x1 && d < x2
    where (x1,x2) = pathValidity path
 
@@ -50,9 +49,9 @@
 bisect :: 
    Path e 
    -> (Geodetic e -> Ordering)        -- ^ Evaluation function.
-   -> Length Double                   -- ^ Required accuracy in terms of distance along the path.
-   -> Length Double -> Length Double  -- ^ Initial bounds.
-   -> Maybe (Length Double)
+   -> Double                          -- ^ Required accuracy in terms of distance along the path.
+   -> Double -> Double                -- ^ Initial bounds.
+   -> Maybe Double
 bisect path f t b1 b2 = do
       guard $ pathValidAt path b1
       guard $ pathValidAt path b2
@@ -65,7 +64,7 @@
       hasRoot (v1, v2) = snd v1 <= EQ && EQ <= snd v2
       sortPair (v1, v2) = if snd v1 <= snd v2 then (v1, v2) else (v2, v1)
       bisect1 ((d1, r1), (d2, r2)) =
-         let d3 = (d1 + d2) / _2
+         let d3 = (d1 + d2) / 2
              r3 = f' d3
              c1 = ((d1, r1), (d3, r3))
              c2 = ((d3, r3), (d2, r2))
@@ -85,16 +84,16 @@
 --
 -- If either estimate departs from its path validity then @Nothing@ is returned.
 intersect :: (Ellipsoid e) =>
-   Length Double -> Length Double     -- ^ Starting estimates.
-   -> Length Double                   -- ^ Required accuracy.
+   Double -> Double                   -- ^ Starting estimates.
+   -> Double                          -- ^ Required accuracy.
    -> Int                             -- ^ Iteration limit. Returns @Nothing@ if this is reached.  
    -> Path e -> Path e                -- ^ Paths to intersect.
-   -> Maybe (Length Double, Length Double)
+   -> Maybe (Double, Double)
 intersect d1 d2 accuracy n path1 path2
    | not $ pathValidAt path1 d1     = Nothing
    | not $ pathValidAt path2 d2     = Nothing
    | n <= 0                         = Nothing
-   | mag < (1e-15 *~ one)           = Nothing
+   | mag < 1e-15                    = Nothing
    | mag3 (nv1 `cross3` nv2) * r <= accuracy = Just (d1, d2)
        -- Assumes that sin (accuracy/r) == accuracy/r
    | otherwise = 
@@ -104,8 +103,8 @@
    where
       (pt1, h1, _) = pathFunc path1 d1
       (pt2, h2, _) = pathFunc path2 d2
-      vectors :: Angle Double -> Angle Double -> Angle Double 
-                 -> (Vec3 (Dimensionless Double), Vec3 (Dimensionless Double))
+      vectors :: Double -> Double -> Double 
+                 -> (Vec3 Double, Vec3 Double)
       vectors lat lon b = (
           -- Unit vector of normal to surface at (lat,lon)
          (cosLat*cosLon, cosLat*sinLon, sinLat),
@@ -125,7 +124,7 @@
       (nv2, gc2) = vectors (latitude pt2) (longitude pt2) h2
       nv3 = gc1 `cross3` gc2         -- Intersection of the great circles
       mag = mag3 nv3
-      nv3a = scale3 nv3 (_1 / mag)   -- Scale to unit. See outer function for case when mag3 == 0
+      nv3a = scale3 nv3 (1 / mag)   -- Scale to unit. See outer function for case when mag3 == 0
       nv3b = negate3 nv3a            -- Antipodal result. Take the closest.
       -- Find "nearest" intersection, defined as smaller of sum of distances to current points.
       d1a = gcDist gc1 nv1 nv3a * r
@@ -135,7 +134,7 @@
       -- Signed angle between v1 and v2, 
       gcDist norm v1 v2 = 
          let c = v1 `cross3` v2 
-         in (if c `dot3` norm < _0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) 
+         in (if c `dot3` norm < 0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) 
       r = majorRadius $ ellipsoid pt1
           
 {- Note on derivation of "intersect"
@@ -170,8 +169,8 @@
 -- | A ray from a point heading in a straight line in 3 dimensions. 
 rayPath :: (Ellipsoid e) => 
    Geodetic e          -- ^ Start point.
-   -> Angle Double     -- ^ Bearing.
-   -> Angle Double     -- ^ Elevation.
+   -> Double           -- ^ Bearing.
+   -> Double           -- ^ Elevation.
    -> Path e
 rayPath pt1 bearing elevation = Path ray alwaysValid
    where
@@ -181,14 +180,14 @@
             (lat,long,alt) = earthToGeo (ellipsoid pt1) pt2'  -- Geodetic of result point.
             (dE,dN,dU) = transform3 (trans3 $ ecefMatrix lat long) delta  -- Direction of ray at result point.
             elevation2 = asin dU
-            bearing2 = if dE == _0 && dN == _0 then bearing else atan2 dE dN  -- Allow for vertical elevation.
+            bearing2 = if dE == 0 && dN == 0 then bearing else atan2 dE dN  -- Allow for vertical elevation.
             
       ecefMatrix lat long =   -- Transform matrix for vectors from (East, North, Up) to (X,Y,Z).
          ((negate sinLong, negate cosLong*sinLat, cosLong*cosLat),
               --    East X      North X               Up X
           (       cosLong, negate sinLong*sinLat, sinLong*cosLat),
               --    East Y      North Y               Up Y
-          (  _0           ,      cosLat         , sinLat))
+          (             0,      cosLat         , sinLat))
               --    East Z      North Z               Up Z
          where
             sinLong = sin long
@@ -215,24 +214,24 @@
 -- the approximation is accurate to within a few meters over 1000km.
 rhumbPath :: (Ellipsoid e) =>
    Geodetic e            -- ^ Start point.
-   -> Angle Double       -- ^ Course.
+   -> Double             -- ^ Course.
    -> Path e
 rhumbPath pt course = Path rhumb validity
    where
-      rhumb distance = (Geodetic lat (properAngle lon) _0 (ellipsoid pt), course, _0)
+      rhumb distance = (Geodetic lat (properAngle lon) 0 (ellipsoid pt), course, 0)
          where
             lat' = lat0 + distance * cosC / m0   -- Kaplan Eq 13.
-            lat = lat0 + (m0 / (a*(_1-e2))) * ((_1-_3*e2/_4)*(lat'-lat0)
-                                              + (_3*e2/_8)*(sin (_2*lat') - sin (_2*lat0)))
-            lon | abs cosC > 1e-7 *~ one 
+            lat = lat0 + (m0 / (a*(1-e2))) * ((1-3*e2/4)*(lat'-lat0)
+                                            + (3*e2/8)*(sin (2*lat') - sin (2*lat0)))
+            lon | abs cosC > 1e-7
                      = lon0 + tanC * (q lat - q0)     -- Kaplan Eq 16.
                 | otherwise
-                     = lon0 + distance * sinC / latitudeRadius (ellipsoid pt) ((lat0 + lat')/_2)
+                     = lon0 + distance * sinC / latitudeRadius (ellipsoid pt) ((lat0 + lat')/2)
       validity
-         | cosC > _0  = ((negate pi/_2 - latitude pt) * b / cosC, (pi/_2 - latitude pt) * b / cosC)
-         | otherwise  = ((pi/_2 - latitude pt) * b / cosC, (negate pi/_2 - latitude pt) * b / cosC)
+         | cosC > 0  = ((negate pi/2 - latitude pt) * b / cosC, (pi/2 - latitude pt) * b / cosC)
+         | otherwise  = ((pi/2 - latitude pt) * b / cosC, (negate pi/2 - latitude pt) * b / cosC)
       q0 = q lat0
-      q phi = log (tan (pi/_4+phi/_2)) + e * log ((_1-eSinPhi)/(_1+eSinPhi)) / _2
+      q phi = log (tan (pi/4+phi/2)) + e * log ((1-eSinPhi)/(1+eSinPhi)) / 2
          where                                -- Factor out expression from Eq 16 of Kaplan
             eSinPhi = e * sin phi
       sinC = sin course
@@ -255,11 +254,11 @@
    -> Path e
 latitudePath pt = Path line alwaysValid
    where
-      line distance = (pt2, pi/_2, _0) 
+      line distance = (pt2, pi/2, 0)
          where
             pt2 = Geodetic 
                (latitude pt) (longitude pt + distance / r)
-               _0 (ellipsoid pt)
+               0 (ellipsoid pt)
       r = latitudeRadius (ellipsoid pt) (latitude pt)
 
 
@@ -270,4 +269,4 @@
 longitudePath :: (Ellipsoid e) =>
    Geodetic e    -- ^ Start point.
    -> Path e
-longitudePath pt = rhumbPath pt _0
+longitudePath pt = rhumbPath pt 0
diff --git a/src/Geodetics/Stereographic.hs b/src/Geodetics/Stereographic.hs
--- a/src/Geodetics/Stereographic.hs
+++ b/src/Geodetics/Stereographic.hs
@@ -10,30 +10,28 @@
    mkGridStereo
 ) where
 
-
 import Geodetics.Ellipsoids
 import Geodetics.Geodetic
 import Geodetics.Grid
-import Numeric.Units.Dimensional.Prelude
-import Prelude ()
 
+import qualified Data.Stream as Stream
 
 -- | A stereographic projection with its origin at an arbitrary point on Earth, other than the poles.
 data GridStereo e = GridStereo {
       gridTangent :: Geodetic e, -- ^ Point where the plane of projection touches the ellipsoid. Often known as the Natural Origin.
       gridOrigin :: GridOffset,  -- ^ Grid position of the tangent point. Often known as the False Origin.
-      gridScale :: Dimensionless Double, -- ^ Scaling factor that balances the distortion between the center and the edges. 
+      gridScale :: Double, -- ^ Scaling factor that balances the distortion between the center and the edges.
                                          -- Should be slightly less than unity.
       
       -- Memoised parameters derived from the tangent point.
-      gridR :: Length Double,
-      gridN, gridC, gridSin, gridCos :: Dimensionless Double,
-      gridLatC :: Angle Double,
-      gridG, gridH :: Length Double
+      gridR :: Double,
+      gridN, gridC, gridSin, gridCos :: Double,
+      gridLatC :: Double,
+      gridG, gridH :: Double
    } deriving (Show)
    
 -- | Create a stereographic projection. The tangency point must not be one of the poles.  
-mkGridStereo :: (Ellipsoid e) => Geodetic e -> GridOffset -> Dimensionless Double -> GridStereo e
+mkGridStereo :: (Ellipsoid e) => Geodetic e -> GridOffset -> Double -> GridStereo e
 mkGridStereo tangent origin scale = GridStereo {
       gridTangent = tangent,
       gridOrigin = origin,
@@ -42,7 +40,7 @@
       gridN = n,
       gridC = c,
       gridSin = sinLatC1,
-      gridCos = sqrt $ _1 - sinLatC1 * sinLatC1,
+      gridCos = sqrt $ 1 - sinLatC1 * sinLatC1,
       gridLatC = asin sinLatC1,
       gridG = g,
       gridH = h
@@ -51,43 +49,43 @@
       -- The reference seems to use χO to refer to two slightly different values. 
       -- Here these will be called LatC0 and LatC1.
       ellipse = ellipsoid tangent
-      op :: Num a => Quantity d a -> Quantity d a    -- Values of longitude, tangent longitude, E and N
-      op = if latitude tangent < _0 then negate else id  -- must be negated in the southern hemisphere.
+      op :: Num a => a -> a    -- Values of longitude, tangent longitude, E and N
+      op = if latitude tangent < 0 then negate else id  -- must be negated in the southern hemisphere.
       lat0 = op $ latitude tangent
       sinLat0 = sin lat0
       e2 = eccentricity2 ellipse
       e = sqrt e2
       r = sqrt $ meridianRadius ellipse lat0 * primeVerticalRadius ellipse lat0
-      n = sqrt $ _1 + ((e2 * cos lat0 ^ pos4)/(_1 - e2))
-      s1 = (_1 + sinLat0) / (_1 - sinLat0)
-      s2 = (_1 - e * sinLat0) / (_1 + e * sinLat0)
+      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
-      sinLatC0 = (w1 - _1)/(w1 + _1)
-      c = ((n + sin lat0) * (_1 - sinLatC0)) / ((n - sin lat0) * (_1 + sinLatC0))
+      sinLatC0 = (w1 - 1)/(w1 + 1)
+      c = ((n + sin lat0) * (1 - sinLatC0)) / ((n - sin lat0) * (1 + sinLatC0))
       w2 = c * w1
-      sinLatC1 = (w2 - _1)/(w2 + _1)
-      g = _2 * r * scale * tan (pi/_4 - latC1/_2)
-      h = _4 * r * scale * tan latC1 + g
+      sinLatC1 = (w2 - 1)/(w2 + 1)
+      g = 2 * r * scale * tan (pi/4 - latC1/2)
+      h = 4 * r * scale * tan latC1 + g
       latC1 = asin sinLatC1
       
 
 instance (Ellipsoid e) => GridClass (GridStereo e) e where
    toGrid grid geo = applyOffset (gridOrigin grid) $ GridPoint east north (geoAlt geo) grid
       where
-         op :: Num a => Quantity d a -> Quantity d a    -- Values of longitude, tangent longitude, E and N
-         op = if latitude (gridTangent grid) < _0 then negate else id  -- must be negated in the southern hemisphere.
-         sinLatC = (w - _1)/(w + _1)
-         cosLatC = sqrt $ _1 - sinLatC * sinLatC
+         op :: Num a => a -> a    -- Values of longitude, tangent longitude, E and N
+         op = if latitude (gridTangent grid) < 0 then negate else id  -- must be negated in the southern hemisphere.
+         sinLatC = (w - 1)/(w + 1)
+         cosLatC = sqrt $ 1 - sinLatC * sinLatC
          longC = gridN grid * (op (longitude geo) - long0) + long0
          w = gridC grid * (sA * sB ** e) ** gridN grid
-         sA = (_1+sinLat) / (_1 - sinLat)
-         sB = (_1 - e*sinLat) / (_1 + e*sinLat)
+         sA = (1+sinLat) / (1 - sinLat)
+         sB = (1 - e*sinLat) / (1 + e*sinLat)
          sinLat = sin $ op $ latitude geo
          e = sqrt $ eccentricity2 $ ellipsoid geo
          long0 = op $ longitude $ gridTangent grid
-         b = _1 + sinLatC * gridSin grid + cosLatC * gridCos grid * cos (longC - long0)
-         east = _2 * gridR grid * gridScale grid * cosLatC * sin (longC - long0) / b
-         north = _2 * gridR grid * gridScale grid * (sinLatC * gridCos grid - cosLatC * gridSin grid * cos (longC - long0)) / b
+         b = 1 + sinLatC * gridSin grid + cosLatC * gridCos grid * cos (longC - long0)
+         east = 2 * gridR grid * gridScale grid * cosLatC * sin (longC - long0) / b
+         north = 2 * gridR grid * gridScale grid * (sinLatC * gridCos grid - cosLatC * gridSin grid * cos (longC - long0)) / b
    
    fromGrid gp = 
       {- trace (    -- Remove comment brackets for debugging.
@@ -97,8 +95,8 @@
          "\n   lat1 = " ++ show lat1 ++ "\n   latN = " ++ show latN ) $ -}
          Geodetic (op latN) (op long) height $ gridEllipsoid grid
       where
-         op :: Num a => Quantity d a -> Quantity d a                   -- Values of longitude, tangent longitude, E and N
-         op = if latitude (gridTangent grid) < _0 then negate else id  -- must be negated in the southern hemisphere.
+         op :: Num a => a -> a                   -- Values of longitude, tangent longitude, E and N
+         op = if latitude (gridTangent grid) < 0 then negate else id  -- must be negated in the southern hemisphere.
          GridPoint east north height _ = applyOffset (offsetNegate $ gridOrigin grid) gp
          east' = east
          north' = north
@@ -106,16 +104,15 @@
          long0 = op $ longitude $ gridTangent grid
          i = atan2 east' (gridH grid + north')
          j = atan2 east' (gridG grid - north') - i
-         latC = gridLatC grid + _2 * atan2 (north' - east' * tan (j/_2)) (_2 * gridR grid * gridScale grid)
-         longC = j + _2 * i + long0
+         latC = gridLatC grid + 2 * atan2 (north' - east' * tan (j/2)) (2 * gridR grid * gridScale grid)
+         longC = j + 2 * i + long0
          sinLatC = sin latC
          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 ^ pos2) / (_1 - e2)
+         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)
             where isoN = isometricLatitude (gridEllipsoid grid) lat
                   e2 = eccentricity2 $ gridEllipsoid grid
-         lats = iterate next lat1
-         latN = snd $ head $ dropWhile (\(v1, v2) -> abs (v1-v2) > 0.01 *~ arcsecond) $ zip lats $ tail lats 
-            
+         lats = Stream.iterate next lat1
+         latN = snd $ Stream.head $ Stream.dropWhile (\(v1, v2) -> abs (v1-v2) > 0.01 * arcsecond) $ Stream.zip lats $ Stream.drop 1 lats
    gridEllipsoid = ellipsoid . gridTangent
diff --git a/src/Geodetics/TransverseMercator.hs b/src/Geodetics/TransverseMercator.hs
--- a/src/Geodetics/TransverseMercator.hs
+++ b/src/Geodetics/TransverseMercator.hs
@@ -5,60 +5,60 @@
    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 ()
 
+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 {
    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 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 
+   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.
-   gridN1, gridN2, gridN3, gridN4 :: Dimensionless Double
+   gridN1, gridN2, gridN3, gridN4 :: Double
 } deriving (Show)
 
 
 -- | Create a Transverse Mercator grid.
-mkGridTM :: (Ellipsoid e) => 
+mkGridTM :: (Ellipsoid e) =>
    Geodetic e               -- ^ True origin.
    -> GridOffset            -- ^ Vector from true origin to false origin.
-   -> Dimensionless Double  -- ^ Scale factor.
+   -> 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
+           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 
+    where
        f = flattening $ ellipsoid origin
-       n = f / (_2-f)  -- Equivalent to (a-b)/(a+b) where b = (1-f)*a
-
-
+       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 
+m :: (Ellipsoid e) => GridTM e -> Double -> 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))
+                    + 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)
@@ -66,77 +66,90 @@
 
 
 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)
-            
-            
+   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 $ head $ dropWhile ((> 0.01 *~ milli meter) . snd) 
-               $ tail $ iterate next (latitude $ trueOrigin grid, 1 *~ meter) 
+         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) 
-               -- head and tail are safe because iterate returns an infinite list.
-          
+               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 ^ pos2
-         v = aF0 / sqrt (_1 - e2 * sinLat2)
-         rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one)
-         eta2 = v / rho - _1
-               
-               
+         sinLat2 = sinLat * sinLat
+         v = aF0 / sqrt (1 - e2 * sinLat2)
+         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) *~ one)
-         eta2 = v / rho - _1
+         v = aF0 / sqrt (1 - e2 * sinLat2)
+         rho = v * (1 - e2) / (1 - e2 * sinLat2)
+         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)
+                   + 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*~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_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_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.
+         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.
+         {-
          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",
@@ -151,8 +164,8 @@
          sinLat = sin lat
          cosLat = cos lat
          tanLat = tan lat
-         sinLat2 = sinLat ^ pos2
-         aF0 = (majorRadius $ gridEllipsoid grid) * gridScale grid
+         sinLat2 = sinLat * sinLat
+         aF0 = majorRadius (gridEllipsoid grid) * gridScale grid
          e2 = eccentricity2 $ gridEllipsoid grid
-         
+
    gridEllipsoid = ellipsoid . trueOrigin
diff --git a/src/Geodetics/UK.hs b/src/Geodetics/UK.hs
--- a/src/Geodetics/UK.hs
+++ b/src/Geodetics/UK.hs
@@ -9,38 +9,33 @@
    toUkGridReference
 ) where
 
-import Control.Applicative
 import Control.Monad
 import Data.Array
 import Data.Char
-import Data.Monoid
 import Geodetics.Geodetic
 import Geodetics.Grid
 import Geodetics.Ellipsoids
 import Geodetics.TransverseMercator
-import Numeric.Units.Dimensional.Prelude
-import qualified Prelude as P
 
 
 
--- | Ellipsoid definition for Great Britain. Airy 1830 offset from the centre of the Earth 
+-- | Ellipsoid definition for Great Britain. Airy 1830 offset from the centre of the Earth
 -- and rotated slightly.
--- 
--- The Helmert parameters are from the Ordnance Survey document 
+--
+-- The Helmert parameters are from the Ordnance Survey document
 -- \"A Guide to Coordinate Systems in Great Britain\", which notes that it
 -- can be in error by as much as 5 meters and should not be used in applications
--- requiring greater accuracy.  A more precise conversion requires a large table 
+-- requiring greater accuracy.  A more precise conversion requires a large table
 -- of corrections for historical inaccuracies in the triangulation of the UK.
 data OSGB36 = OSGB36 deriving (Eq, Show)
 
 instance Ellipsoid OSGB36 where
-   majorRadius _ = 6377563.396 *~ meter
-   flatR _ = 299.3249646 *~ one
+   majorRadius _ = 6377563.396
+   flatR _ = 299.3249646
    helmert _ = Helmert {
-      cX = 446.448 *~ meter, cY = (-125.157) *~ meter, cZ = 542.06 *~ meter,
-      helmertScale = (-20.4894) *~ one,
-      rX = 0.1502 *~ arcsecond, rY = 0.247 *~ arcsecond, rZ = 0.8421 *~ arcsecond }
-
+      cX = 446.448, cY = (-125.157), cZ = 542.06,
+      helmertScale = (-20.4894),
+      rX = 0.1502 * arcsecond, rY = 0.247 * arcsecond, rZ = 0.8421 * arcsecond }
 
 -- | The UK National Grid is a Transverse Mercator projection with a true origin at
 -- 49 degrees North, 2 degrees West on OSGB36, and a false origin 400km West and 100 km North of
@@ -56,73 +51,74 @@
 
 ukTrueOrigin :: Geodetic OSGB36
 ukTrueOrigin = Geodetic {
-   latitude = 49 *~ degree,
-   longitude = (-2) *~ degree,
-   geoAlt = 0 *~ meter,
+   latitude = 49 * degree,
+   longitude = (-2) * degree,
+   geoAlt = 0,
    ellipsoid = OSGB36
 }
 
-ukFalseOrigin :: GridOffset 
-ukFalseOrigin = GridOffset ((-400) *~ kilo meter) (100 *~ kilo meter) (0 *~ meter)
+ukFalseOrigin :: GridOffset
+ukFalseOrigin = GridOffset ((-400) * kilometer) (100 * kilometer) (0 * kilometer)
 
 
 -- | Numerical definition of the UK national grid.
 ukGrid :: GridTM OSGB36
-ukGrid = mkGridTM ukTrueOrigin ukFalseOrigin 
-   ((10 *~ one) ** (0.9998268 *~ one - _1))
+ukGrid = mkGridTM ukTrueOrigin ukFalseOrigin
+   (10 ** (0.9998268 - 1))
 
 
 -- | Size of a UK letter-pair grid square.
-ukGridSquare :: Length Double
-ukGridSquare = 100 *~ kilo meter
+ukGridSquare :: Double
+ukGridSquare = 100 * kilometer
 
 
--- | Convert a grid reference to a position, if the reference is valid. 
--- This returns the position of the south-west corner of the nominated 
+-- | Convert a grid reference to a position, if the reference is valid.
+-- This returns the position of the south-west corner of the nominated
 -- grid square and an offset to its centre. Altitude is set to zero.
 fromUkGridReference :: String -> Maybe (GridPoint UkNationalGrid, GridOffset)
-fromUkGridReference str = if length str < 2 then Nothing else do
-      let 
-         c1:c2:ds = str
-         n = length ds
-      guard $ even n
-      let (dsE, dsN) = splitAt (n `div` 2) ds
-      (east, sq) <- fromGridDigits ukGridSquare dsE
-      (north, _) <- fromGridDigits ukGridSquare dsN
-      base <- fromUkGridLetters c1 c2
-      let half = sq / (2 *~ one)
-      return (applyOffset (GridOffset east north (0 *~ meter)) base,
-              GridOffset half half (0 *~ meter))
+fromUkGridReference str =
+      case str of
+         c1:c2:ds -> do
+            let n = length ds
+            guard $ even n
+            let (dsE, dsN) = splitAt (n `div` 2) ds
+            (east, sq) <- fromGridDigits ukGridSquare dsE
+            (north, _) <- fromGridDigits ukGridSquare dsN
+            base <- fromUkGridLetters c1 c2
+            let half = sq / 2
+            return (applyOffset (GridOffset east north 0) base,
+                  GridOffset half half 0)
+         _ -> Nothing
 
-      
 
 
+
 -- | The south west corner of the nominated grid square, if it is a legal square.
 -- This function works for all pairs of letters except 'I' (as that is not used).
 -- In practice only those pairs covering the UK are actually considered meaningful.
 fromUkGridLetters :: Char -> Char -> Maybe (GridPoint UkNationalGrid)
 fromUkGridLetters c1 c2 = applyOffset <$> (mappend <$> g1 <*> g2) <*> letterOrigin
    where
-      letterOrigin = Just $ GridPoint ((-1000) *~ kilo meter) ((-500) *~ kilo meter) m0 UkNationalGrid
-      gridIndex c = 
-         if inRange ('A', 'H') c then Just $ ord c P.- ord 'A'  -- 'I' is not used.
-         else if inRange ('J', 'Z') c then Just $ ord c P.- ord 'B'
+      letterOrigin = Just $ GridPoint ((-1000) * kilometer) ((-500) * kilometer) m0 UkNationalGrid
+      gridIndex c =
+         if inRange ('A', 'H') c then Just $ ord c - ord 'A'  -- 'I' is not used.
+         else if inRange ('J', 'Z') c then Just $ ord c - ord 'B'
          else Nothing
       gridSquare c = do -- Maybe monad
          g <- gridIndex c
-         let (y,x) = g `divMod` 5 
-         return (fromIntegral x *~ one, _4 - fromIntegral y *~ one)
+         let (y,x) = g `divMod` 5
+         return (fromIntegral x, 4 - fromIntegral y)
       g1 = do
          (x,y) <- gridSquare c1
-         return $ GridOffset (x * (500 *~ kilo meter)) (y * (500 *~ kilo meter)) m0
+         return $ GridOffset (x * (500 * kilometer)) (y * (500 * kilometer)) m0
       g2 = do
          (x,y) <- gridSquare c2
-         return $ GridOffset (x * (100 *~ kilo meter)) (y * (100 *~ kilo meter)) m0
-      m0 = 0 *~ meter
+         return $ GridOffset (x * (100 * kilometer)) (y * (100 * kilometer)) m0
+      m0 = 0
 
 
 -- | Find the nearest UK grid reference point to a specified position. The Int argument is the number of
--- digits precision, so 2 for a 4-figure reference and 3 for a 6-figure reference, although any value 
+-- digits precision, so 2 for a 4-figure reference and 3 for a 6-figure reference, although any value
 -- between 0 and 5 can be used (giving a 1 meter precision).
 -- Altitude is ignored. If the result is outside the area defined by the two letter grid codes then
 -- @Nothing@ is returned.
@@ -130,16 +126,15 @@
 toUkGridReference n p
    | n < 0         = error "toUkGridReference: precision argument must not be negative."
    | otherwise     = do
-      (gx, strEast) <- toGridDigits ukGridSquare n $ eastings p + 1000 *~ kilo meter
-      (gy, strNorth) <- toGridDigits ukGridSquare n $ northings p + 500 *~ kilo meter
-      let (gx1, gx2) = (fromIntegral gx) `divMod` 5
-          (gy1, gy2) = (fromIntegral gy) `divMod` 5
+      (gx, strEast) <- toGridDigits ukGridSquare n $ eastings p + 1000 * kilometer
+      (gy, strNorth) <- toGridDigits ukGridSquare n $ northings p + 500 * kilometer
+      let (gx1, gx2) = fromIntegral gx `divMod` 5
+          (gy1, gy2) = fromIntegral gy `divMod` 5
       guard (gx1 < 5 && gy1 < 5)
       let c1 = gridSquare gx1 gy1
           c2 = gridSquare gx2 gy2
       return $ c1 : c2 : strEast ++ strNorth
    where
-      gridSquare x y = letters ! (4 P.- y, x)
+      gridSquare x y = letters ! (4 - y, x)
       letters :: Array (Int, Int) Char
       letters = listArray ((0,0),(4,4)) $ ['A'..'H'] ++ ['J'..'Z']
-   
diff --git a/test/ArbitraryInstances.hs b/test/ArbitraryInstances.hs
--- a/test/ArbitraryInstances.hs
+++ b/test/ArbitraryInstances.hs
@@ -5,7 +5,6 @@
 
 module ArbitraryInstances where
 
-import Control.Applicative
 import Control.Monad
 import Geodetics.Altitude
 import Geodetics.Geodetic
@@ -14,67 +13,67 @@
 import Geodetics.Path
 import Geodetics.Stereographic as SG
 import Geodetics.TransverseMercator as TM
-import Numeric.Units.Dimensional.Prelude
-import qualified Prelude ()
 import Test.QuickCheck
 
 
+-- | Shrink an angle so that shrunk values are round numbers of degrees.
+shrinkAngle :: Double -> [Double]
+shrinkAngle v = (degree *) <$> shrink (v/degree)
 
--- | Shrink using a dimension, so that shrunk values are round numbers in that dimension.
-shrinkDimension :: forall a (d :: Dimension) (m :: Metricality) .
-                   (Fractional a, Arbitrary a) => Unit m d a -> Quantity d a -> [Quantity d a]
-shrinkDimension u v = (*~ u) <$> shrink (v /~ u)
+-- | Shrink a distance so that shrunk values are round numbers of kilometers.
+shrinkDistance :: Double -> [Double]
+shrinkDistance v = (kilometer *) <$> shrink (v/kilometer)
 
 -- | Wrapper for arbitrary angles.
-newtype Bearing = Bearing (Dimensionless Double)
+newtype Bearing = Bearing Double
 
 instance Show Bearing where
    show (Bearing b) = "Bearing " ++ showAngle b
 
 instance Arbitrary Bearing where
-   arbitrary = Bearing <$> (*~ degree) <$> choose (-180,180)
-   shrink (Bearing b) = Bearing <$> shrinkDimension degree b
+   arbitrary = Bearing <$> (* degree) <$> choose (-180,180)
+   shrink (Bearing b) = Bearing <$> shrinkAngle b
    
    
-newtype Azimuth = Azimuth (Dimensionless Double)
+newtype Azimuth = Azimuth Double
 
 instance Show Azimuth where
    show (Azimuth a) = "Azimuth " ++ showAngle a
    
 instance Arbitrary Azimuth where
-   arbitrary = Azimuth <$> (*~ degree) <$> choose (0,90)
-   shrink (Azimuth a) = Azimuth <$> shrinkDimension degree a
+   arbitrary = Azimuth <$> (* degree) <$> choose (0,90)
+   shrink (Azimuth a) = Azimuth <$> shrinkAngle a
    
    
 -- | Wrapper for arbitrary distances up to 10,000 km
-newtype Distance = Distance (Length Double) deriving (Show)
+newtype Distance = Distance Double deriving (Show)
 
 instance Arbitrary Distance where
-   arbitrary = Distance <$> (*~ kilo meter) <$> choose (0,10000)
-   shrink (Distance d) = Distance <$> shrinkDimension (kilo meter) d
+   arbitrary = Distance <$> (* kilometer) <$> choose (0,10000)
+   shrink (Distance d) = Distance <$> shrinkDistance d
    
 
 -- | Wrapper for arbitrary distances up to 1,000 km
-newtype Distance2 = Distance2 (Length Double) deriving (Show)
+newtype Distance2 = Distance2 Double deriving (Show)
 
 instance Arbitrary Distance2 where
-   arbitrary = Distance2 <$> (*~ kilo meter) <$> choose (0,1000)
-   shrink (Distance2 d) = Distance2 <$> shrinkDimension (kilo meter) d
+   arbitrary = Distance2 <$> (* kilometer) <$> choose (0,1000)
+   shrink (Distance2 d) = Distance2 <$> shrinkDistance d
 
 -- | Wrapper for arbitrary altitudes up to 10 km
-newtype Altitude = Altitude (Length Double) deriving (Show)
+newtype Altitude = Altitude Double deriving (Show)
 
 instance Arbitrary Altitude where
-   arbitrary = Altitude <$> (*~ kilo meter) <$> choose (0,10)
-   shrink (Altitude h) = Altitude <$> shrinkDimension (kilo meter) h
+   arbitrary = Altitude <$> (* kilometer) <$> choose (0,10)
+   shrink (Altitude h) = Altitude <$> shrinkDistance h
 
 
 -- | Wrapper for arbitrary dimensionless numbers (-10 .. 10)
-newtype Scalar = Scalar (Dimensionless Double) deriving (Show)
+newtype Scalar = Scalar Double deriving (Show)
 
 instance Arbitrary Scalar where
-   arbitrary = Scalar <$> (*~ one) <$> choose (-10,10)
-   shrink (Scalar s) = Scalar <$> shrinkDimension one s
+   arbitrary = Scalar <$> choose (-10,10)
+   shrink (Scalar s) = Scalar <$> shrink s
 
 
 -- | Wrapper for arbitrary grid references.
@@ -92,21 +91,21 @@
 
 
 -- | Generate in range +/- <arg> m.
-genOffset :: Double -> Gen (Length Double)
-genOffset d = (*~ meter) <$> choose (-d, d)
+genOffset :: Double -> Gen Double
+genOffset d = choose (-d, d)
 
-genAlt :: Gen (Length Double)
-genAlt = (*~ meter) <$> choose (0,10000)
+genAlt :: Gen Double
+genAlt = choose (0,10000)
 
 
-genLatitude :: Gen (Dimensionless Double)
-genLatitude = (*~ degree) <$> choose (-90,90)
+genLatitude :: Gen Double
+genLatitude = (* degree) <$> choose (-90,90)
 
-genLongitude :: Gen (Dimensionless Double)
-genLongitude = (*~ degree) <$> choose (-180,180)
+genLongitude :: Gen Double
+genLongitude = (* degree) <$> choose (-180,180)
 
-genSeconds :: Gen (Dimensionless Double)
-genSeconds = (*~ arcsecond) <$> choose (-10,10)
+genSeconds :: Gen Double
+genSeconds = (* arcsecond) <$> choose (-10,10)
     
 
 -- | Shrinking with the original value preserved. Used for shrinking records.  See 
@@ -114,30 +113,18 @@
 shrink' :: (Arbitrary a) => a -> [a]
 shrink' x = x : shrink x
 
--- | Shrink a quantity in the given units.
-shrinkQuantity :: forall a (d :: Dimension) (m :: Metricality).
-                  (Arbitrary a, Fractional a) => Unit m d a -> Quantity d a -> [Quantity d a]
-shrinkQuantity u q = map (*~ u) $ shrink' $ q /~ u
-
-shrinkLength :: (Arbitrary a, Fractional a) => Length a -> [Length a]
-shrinkLength = shrinkQuantity meter
-
-shrinkUnit :: (Arbitrary a, Fractional a) => Dimensionless a -> [Dimensionless a]
-shrinkUnit = shrinkQuantity one
-
-shrinkAngle :: (Arbitrary a, Floating a) => Dimensionless a -> [Dimensionless a]
-shrinkAngle = shrinkQuantity degree
+shrinkAngle' :: Double -> [Double]
+shrinkAngle' a = a : shrinkAngle a
 
 
 instance Arbitrary Helmert where
    arbitrary = 
       Helmert <$> genOffset 300 <*> genOffset 300 <*> genOffset 300 <*> 
-         ((*~ one) <$> choose (-5,10)) <*>
-         genSeconds <*> genSeconds <*> genSeconds
+         (choose (-5,10)) <*> genSeconds <*> genSeconds <*> genSeconds
    shrink h = 
-      tail $ Helmert <$> shrinkLength (cX h) <*> shrinkLength (cY h) <*> shrinkLength (cZ h) <*>
-         shrinkUnit (helmertScale h) <*>
-         shrinkUnit (rX h) <*> shrinkUnit (rY h) <*> shrinkUnit (rZ h)      
+      drop 1 $ Helmert <$> shrink' (cX h) <*> shrink' (cY h) <*> shrink' (cZ h) <*>
+         shrink' (helmertScale h) <*>
+         shrink' (rX h) <*> shrink' (rY h) <*> shrink' (rZ h)      
 
 
 instance Arbitrary WGS84 where
@@ -148,55 +135,55 @@
 instance Arbitrary LocalEllipsoid where
    arbitrary =
       LocalEllipsoid <$> (("Local_" ++) <$> replicateM 3 (choose ('A','Z'))) <*>  -- name
-         ((*~ meter) <$> choose (6378100, 6378400)) <*>                  -- majorRadius
-         ((*~ one) <$> choose (297,300)) <*>                             -- flatR
-         arbitrary                                                       -- helmert
-   shrink e = tail $ LocalEllipsoid (nameLocal e) (majorRadius e) (flatR e) <$> shrink' (helmert e)
+         (choose (6378100, 6378400)) <*>                    -- majorRadius
+         (choose (297,300)) <*>                             -- flatR
+         arbitrary                                          -- helmert
+   shrink e = drop 1 $ LocalEllipsoid (nameLocal e) (majorRadius e) (flatR e) <$> shrink' (helmert e)
 
         
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (Geodetic e) where
    arbitrary = 
       Geodetic <$> genLatitude <*> genLongitude <*> genOffset 1 <*> arbitrary
    shrink g = 
-      tail $ Geodetic <$> shrinkAngle (latitude g) <*> shrinkAngle (longitude g) <*> 
-         shrinkLength (altitude g) <*> shrink' (ellipsoid g)
+      drop 1 $ Geodetic <$> shrinkAngle' (latitude g) <*> shrinkAngle' (longitude g) <*> 
+         shrink' (altitude g) <*> shrink' (ellipsoid g)
 
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridTM e)) where
    arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary
-   shrink p = tail $ GridPoint <$> 
-      shrinkLength (eastings p) <*> 
-      shrinkLength (northings p) <*> 
-      shrinkLength (altitude p) <*> 
+   shrink p = drop 1 $ GridPoint <$> 
+      shrink' (eastings p) <*> 
+      shrink' (northings p) <*> 
+      shrink' (altitude p) <*> 
       shrink' (gridBasis p)
 
 
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridStereo e)) where
    arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary
-   shrink p = tail $ GridPoint <$> 
-      shrinkLength (eastings p) <*> 
-      shrinkLength (northings p) <*> 
-      shrinkLength (altitude p) <*> 
+   shrink p = drop 1 $ GridPoint <$> 
+      shrink' (eastings p) <*> 
+      shrink' (northings p) <*> 
+      shrink' (altitude p) <*> 
       shrink' (gridBasis p)
 
 
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridTM e) where
-   arbitrary = mkGridTM <$> arbitrary <*> arbitrary <*> ((*~ one) <$> choose (0.95,1.0))
-   shrink tm = tail $ mkGridTM <$> shrink' (trueOrigin tm) <*> shrink' (falseOrigin tm) <*> [TM.gridScale tm]
+   arbitrary = mkGridTM <$> arbitrary <*> arbitrary <*> choose (0.95,1.0)
+   shrink tm = drop 1 $ mkGridTM <$> shrink' (trueOrigin tm) <*> shrink' (falseOrigin tm) <*> [TM.gridScale tm]
    
    
 instance Arbitrary GridOffset where
    arbitrary = GridOffset <$> genOffset 100000 <*> genOffset 100000 <*> genAlt
-   shrink d = tail $ GridOffset <$> 
-      shrinkLength (deltaEast d) <*> shrinkLength (deltaNorth d) <*> shrinkLength (deltaAltitude d)
+   shrink d = drop 1 $ GridOffset <$> 
+      shrink' (deltaEast d) <*> shrink' (deltaNorth d) <*> shrink' (deltaAltitude d)
 
 
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridStereo e) where
-   arbitrary = mkGridStereo <$> arbitrary <*> arbitrary <*> ((*~ one) <$> choose (0.95,1.0))
-   shrink sg = tail $ mkGridStereo <$> shrink' (gridTangent sg) <*> shrink' (gridOrigin sg) <*> [SG.gridScale sg]
+   arbitrary = mkGridStereo <$> arbitrary <*> arbitrary <*> choose (0.95,1.0)
+   shrink sg = drop 1 $ mkGridStereo <$> shrink' (gridTangent sg) <*> shrink' (gridOrigin sg) <*> [SG.gridScale sg]
    
 
 -- | Wrapper for arbitrary rays, along with creation parameters for printing and shrinking.
-data Ray e = Ray (Geodetic e) (Angle Double) (Angle Double)
+data Ray e = Ray (Geodetic e) Double Double
 
 instance (Ellipsoid e) => Show (Ray e) where
    show (Ray p0 b e ) = "(Ray " ++ show p0 ++ ", " ++ showAngle b ++ ", " ++ showAngle e ++ ")"
@@ -207,13 +194,13 @@
 instance (Ellipsoid e, Arbitrary e) => Arbitrary (Ray e) where
    arbitrary = do
       p0 <- arbitrary
-      b <- (*~ degree) <$> choose (-180,180)
-      e <- (*~ degree) <$> choose (0,90)
+      b <- (* degree) <$> choose (-180,180)
+      e <- (* degree) <$> choose (0,90)
       return $ Ray p0 b e
-   shrink (Ray p0 b e) = tail $ do
+   shrink (Ray p0 b e) = drop 1 $ do
       p0' <- shrink' p0
-      b' <- shrinkAngle b
-      e' <- shrinkAngle e
+      b' <- shrinkAngle' b
+      e' <- shrinkAngle' e
       return $ Ray p0' b' e'
       
      
@@ -230,15 +217,15 @@
    show rp2 = show (pt1, Bearing b1) ++ show (pt2, Bearing b2)
       where 
          (p1, p2) = mk2RhumbPaths rp2
-         (pt1, b1, _) = pathFunc p1 (0 *~ meter)
-         (pt2, b2, _) = pathFunc p2 (0 *~ meter)
+         (pt1, b1, _) = pathFunc p1 0
+         (pt2, b2, _) = pathFunc p2 0
           
 instance Arbitrary RhumbPaths2 where
    arbitrary = RP2 
-      <$> arbitrary `suchThat` ((< 70 *~ degree) . abs . latitude)
+      <$> arbitrary `suchThat` ((< (70 * degree)) . abs . latitude)
       <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
    shrink rp = 
-      tail $ RP2 <$> 
+      drop 1 $ RP2 <$> 
          shrink' (rp2Point0 rp) <*> 
          shrink' (rp2Bearing0 rp) <*> 
          shrink' (rp2Distance rp) <*> 
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -3,9 +3,6 @@
 module Main where
 
 import Data.Maybe
-import Data.Monoid
-import Numeric.Units.Dimensional.Prelude
-import qualified Prelude as P
 import Test.Framework (Test, defaultMainWithOpts, testGroup)
 import Test.Framework.Options (TestOptions, TestOptions'(..))
 import Test.Framework.Runners.Options (RunnerOptions, RunnerOptions'(..))
@@ -44,7 +41,7 @@
 instance EqProp GridOffset where
   (GridOffset a b c) =-= (GridOffset a' b' c') =
     eq True $ a ≈ a' && b ≈ b' && c ≈ c'
-    where x ≈ y = abs (x - y) < 0.00001 *~ meter
+    where x ≈ y = abs (x - y) < 0.00001
 
 instance EqProp Helmert where
   (Helmert cX' cY' cZ' s rX' rY' rZ') =-= (Helmert cX'' cY'' cZ'' s' rX'' rY'' rZ'') =
@@ -52,8 +49,8 @@
                    s ≈- s',
                    rX' ≈- rX'', rY' ≈- rY'', rZ' ≈- rZ'']
 
-    where x ≈ y = abs (x - y) < 0.00001 *~ meter
-          x ≈- y = abs (x - y) < (_1 / (_5 * _2) ** (_5))
+    where x ≈ y = abs (x - y) < 0.00001
+          x ≈- y = abs (x - y) < (1 / (5 * 2) ^ _5)
 
 tests :: [Test]
 tests = [
@@ -67,6 +64,9 @@
       testProperty "Grid Offset 2" prop_offset2,
       testProperty "Grid Offset 3" prop_offset3,
       testProperty "Grid 1" prop_grid1 ],
+   testGroup "TransverseMercator" [
+      testCase "fromGrid . toGrid == id" $ HU.assertBool "" prop_tmGridInverse
+      ],
    testGroup "UK" [
       testProperty "UK Grid 1" prop_ukGrid1,
       testGroup "UK Grid 2" $ map ukGridTest2 ukSampleGrid,
@@ -95,38 +95,38 @@
 
 -- | The positions are within 30 cm.
 samePlace :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool
-samePlace p1 p2 = geometricalDistance p1 p2 < 0.3 *~ meter
+samePlace p1 p2 = geometricalDistance p1 p2 < 0.3
 
 
 -- | The positions are within 10 m.
 closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool
-closeEnough p1 p2 = geometricalDistance p1 p2 < 10 *~ meter
+closeEnough p1 p2 = geometricalDistance p1 p2 < 10
 
 
 -- | The angles are within 0.01 arcsec
-sameAngle :: Angle Double -> Angle Double -> Bool
-sameAngle v1 v2 = abs (properAngle (v1 - v2)) < 0.01 *~ arcsecond
+sameAngle :: Double -> Double -> Bool
+sameAngle v1 v2 = abs (properAngle (v1 - v2)) < 0.01 * arcsecond
 
 -- | The grid positions are within 1mm
 sameGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool
 sameGrid p1 p2 = check eastings && check northings && check altitude
-   where check f = f p1 - f p2 < 1 *~ milli meter
+   where check f = f p1 - f p2 < 1e-3
 
 
 -- | Grid offsets are within 1mm.
 sameOffset :: GridOffset -> GridOffset -> Bool
 sameOffset go1 go2 = check deltaNorth && check deltaEast && check deltaAltitude
-   where check f = f go1 - f go2 < 1 *~ milli meter
+   where check f = f go1 - f go2 < 1e-3
 
 
 -- | The grid X and Y are both within 1 meter
 closeGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool
 closeGrid p1 p2 = check eastings && check northings && check altitude
-   where check f = f p1 - f p2 < 1 *~ meter
+   where check f = f p1 - f p2 < 1
 
 -- | Degrees, minutes and seconds into radians.
-dms :: Int -> Int -> Double -> Dimensionless Double
-dms d m s = fromIntegral d *~ degree + fromIntegral m *~ arcminute + s *~ arcsecond
+dms :: Int -> Int -> Double -> Double
+dms d m s = fromIntegral d * degree + fromIntegral m * arcminute + s * arcsecond
 
 -- | Round-trip from local to WGS84 and back is identity (approximately)
 prop_WGS84_and_back :: Geodetic LocalEllipsoid -> Bool
@@ -138,46 +138,46 @@
 prop_zero_ground p =
    case groundDistance p p of
       Nothing -> False
-      Just (d, _, _) -> abs d < 1 *~ milli meter
+      Just (d, _, _) -> abs d < 1e-3
 
 
 -- | Sample pairs of points with bearings and distances.
 -- The Oracle for these values is the @FORWARD@ program from
 --  <http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html>
-worldLines :: [(String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double)]
+worldLines :: [(String, Geodetic WGS84, Geodetic WGS84, {-Length-} Double, {-Angle-} Double, {-Angle-} Double)]
 worldLines = [
-   ("Ordinary", Geodetic (40*~degree) (30*~degree) _0 WGS84, Geodetic (30*~degree) (50*~degree) _0 WGS84,
-      2128852.999*~meter, 115.19596706*~degree, 126.79044315*~degree),
-   ("Over Pole", Geodetic (60*~degree) (0*~degree) _0 WGS84, Geodetic (60*~degree) (180*~degree) _0 WGS84,
-      6695785.820*~meter, 0*~degree, 180*~degree),
-   ("Equator to Pole", Geodetic (0*~degree) (0*~degree) _0 WGS84, Geodetic (90*~degree) (180*~degree) _0 WGS84,
-      10001965.729*~meter, 0*~degree, 180*~degree)]
+   ("Ordinary", Geodetic (40 * degree) (30 * degree) 0 WGS84, Geodetic (30 * degree) (50 * degree) 0 WGS84,
+      2128852.999, 115.19596706 * degree, 126.79044315 * degree),
+   ("Over Pole", Geodetic (60 * degree) (0 * degree) 0 WGS84, Geodetic (60 * degree) (180 * degree) 0 WGS84,
+      6695785.820, 0 * degree, 180 * degree),
+   ("Equator to Pole", Geodetic (0 * degree) (0 * degree) 0 WGS84, Geodetic (90 * degree) (180 * degree) 0 WGS84,
+      10001965.729, 0 * degree, 180 * degree)]
 
 
-worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double) -> Test
+worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Double, Double, Double) -> Test
 worldLineTests (str, g1, g2, d, a, b) = testCase str $ HU.assertBool "" $ ok $ groundDistance g1 g2
    where
       ok Nothing = False
       ok (Just (d1, a1, b1)) =
-         abs (d - d1) < 0.01 *~ meter
-         && abs (a - a1) < 0.01 *~ arcsecond
-         && abs (b - b1) < 0.01 *~ arcsecond
+         abs (d - d1) < 0.01
+         && abs (a - a1) < 0.01 * arcsecond
+         && abs (b - b1) < 0.01 * arcsecond
 
 -- | Sample points for UK tests. The oracle for these values is the script at
 -- <http://www.movable-type.co.uk/scripts/latlong-convert-coords.html>, which uses
 -- the same Helmert transform as this library. Hence the results should match to within 30 cm.
 ukPoints :: [(String, Geodetic WGS84, Geodetic OSGB36)]
 ukPoints = [
-   ("Greenwich",        Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) _0 WGS84,
-                        Geodetic (dms 51 28 39.00) (dms 0 0 0) _0 OSGB36),
-   ("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) _0 WGS84,
-                        Geodetic (dms 55 56 56.51) (dms (-3) (-11) (-57.61)) _0 OSGB36),
-   ("Lands End",        Geodetic (dms 50 03 56.68) (dms (-5) (-42) (-51.20)) _0 WGS84,
-                        Geodetic (dms 50 03 54.51) (dms (-5) (-42) (-47.87)) _0 OSGB36),
-   ("Gt. Yarmouth Pier",Geodetic (dms 52 36 29.33) (dms 1 44 27.79) _0 WGS84,
-                        Geodetic (dms 52 36 27.84) (dms 1 44 34.52) _0 OSGB36),
-   ("Stanhope",         Geodetic (dms 54 44 49.08) (dms (-2) 0 (-19.89)) _0 WGS84,
-                        Geodetic (dms 54 44 48.71) (dms (-2) 0 (-14.41)) _0 OSGB36) ]
+   ("Greenwich",        Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) 0 WGS84,
+                        Geodetic (dms 51 28 39.00) (dms 0 0 0) 0 OSGB36),
+   ("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) 0 WGS84,
+                        Geodetic (dms 55 56 56.51) (dms (-3) (-11) (-57.61)) 0 OSGB36),
+   ("Lands End",        Geodetic (dms 50 03 56.68) (dms (-5) (-42) (-51.20)) 0 WGS84,
+                        Geodetic (dms 50 03 54.51) (dms (-5) (-42) (-47.87)) 0 OSGB36),
+   ("Gt. Yarmouth Pier",Geodetic (dms 52 36 29.33) (dms 1 44 27.79) 0 WGS84,
+                        Geodetic (dms 52 36 27.84) (dms 1 44 34.52) 0 OSGB36),
+   ("Stanhope",         Geodetic (dms 54 44 49.08) (dms (-2) 0 (-19.89)) 0 WGS84,
+                        Geodetic (dms 54 44 48.71) (dms (-2) 0 (-14.41)) 0 OSGB36) ]
 
 
 
@@ -201,19 +201,34 @@
 prop_offset3 :: GridOffset -> Bool
 prop_offset3 delta = sameOffset delta0
                                 (polarOffset (offsetDistance delta0) (offsetBearing delta))
-   where delta0 = delta {deltaAltitude = 0 *~ meter}
+   where delta0 = delta {deltaAltitude = 0}
 
 -- | Given a grid point and an offset, applying the offset to the point gives a new point which
 -- is offset from the first point by the argument offset.
 prop_grid1 :: GridPoint (GridTM LocalEllipsoid) -> GridOffset -> Bool
 prop_grid1 p d = sameOffset d $ p `gridOffset` applyOffset d p
 
-
+-- | Check that using toGrid/fromGrid for TransverseMercator projection are inverses
+-- | for negative latitudes near the coordinates 0,0
+prop_tmGridInverse :: Bool
+prop_tmGridInverse = 
+   let origin = Geodetic 
+         { latitude = 0 * degree
+         , longitude = 0 * degree
+         , geoAlt = 0
+         , ellipsoid = WGS84
+         }
+       g = mkGridTM origin mempty 1
+       testPoint = origin { latitude = (-1) * arcminute }
+       tp1 = toGrid g testPoint
+       tp2 = fromGrid tp1
+   in tp2 `closeEnough` testPoint
+   
 -- | Converting a UK grid reference to a GridPoint and back is a null operation.
 prop_ukGrid1 :: GridRef -> Bool
 prop_ukGrid1 (GridRef str) =
    str ==
-   (fromJust $ toUkGridReference ((length str P.- 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str)
+   (fromJust $ toUkGridReference ((length str - 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str)
 
 -- | UK Grid Reference points. The oracle for these points was the
 -- UK Grid Reference Finder (gridreferencefinder.com), retrieved on 26 Jan 2013.
@@ -228,11 +243,14 @@
    ("ND3804872787", 338048, 972787, 58.638518, -3.0688688,   "John O Groats"),
    ("SC3915875189", 239158, 475189, 54.147275, -4.4641148,   "Douglas Harbour"),
    ("ST1922474591", 319224, 174591, 51.464505, -3.1641741,   "Torchwood HQ"),
-   ("SK3520736502", 435207, 336502, 52.924784, -1.4777486,   "Derby Cathedral")]
+   ("SK3520736502", 435207, 336502, 52.924784, -1.4777486,   "Derby Cathedral"),
+   ("TG5141013177", 651410, 313177, 52.657979 , 1.7160519,   "Caister Water Tower"),
+   ("TG2623802646", 626238, 302646, 52.574548 , 1.3373749,   "Framingham")]
+   -- Caister and Framingham are taken from Ordnance Survey worked examples.
    where
       convert (grid, x, y, lat, long, desc) =
-         (grid, GridPoint (x *~ meter) (y *~ meter) (0 *~ meter) UkNationalGrid,
-          Geodetic (lat *~ degree) (long *~ degree) (0 *~ meter) WGS84, desc)
+         (grid, GridPoint (x) (y) (0) UkNationalGrid,
+          Geodetic (lat * degree) (long * degree) (0) WGS84, desc)
 
 type GridPointTest = (String, GridPoint UkNationalGrid, Geodetic WGS84, String) -> Test
 
@@ -254,12 +272,12 @@
 -- | Check that WGS84 to grid point works close enough for sample points.
 ukGridTest5 :: GridPointTest
 ukGridTest5 (_, gp, geo, testName) = testCase testName $ HU.assertBool ""
-   $ offsetDistance (gridOffset gp $ toGrid UkNationalGrid $ toLocal OSGB36 geo) < 1 *~ meter
+   $ offsetDistance (gridOffset gp $ toGrid UkNationalGrid $ toLocal OSGB36 geo) < 1
 
 
 -- | Worked example for UK Geodetic to GridPoint, taken from "A Guide to Coordinate Systems in Great Britain" [1]
 ukTest :: Geodetic OSGB36
-ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) (0 *~ meter) OSGB36
+ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) (0) OSGB36
 
 {-
    v = 6.3885023333E+06
@@ -280,22 +298,22 @@
 
 -- | Standard stereographic grid for point tests in the Northern Hemisphere.
 stereoGridN :: GridStereo LocalEllipsoid
-stereoGridN = mkGridStereo tangent origin (0.9999079 *~ one)
+stereoGridN = mkGridStereo tangent origin (0.9999079)
    where
-      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155 *~ metre) (299.15281 *~ one) mempty
-      tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) (0 *~ meter) ellipse
-      origin = GridOffset (155000 *~ metre) (463000 *~ metre) (0 *~ meter)
+      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155) (299.15281) mempty
+      tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) (0) ellipse
+      origin = GridOffset (155000) (463000) (0)
 
 
 -- | Standard steregraphic grid for point tests in the Southern Hemisphere.
 --
 -- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated.
 stereoGridS :: GridStereo LocalEllipsoid
-stereoGridS = mkGridStereo tangent origin (0.9999079 *~ one)
+stereoGridS = mkGridStereo tangent origin (0.9999079)
    where
-      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155 *~ metre) (299.15281 *~ one) mempty
-      tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) (0 *~ meter) ellipse
-      origin = GridOffset ((-155000) *~ metre) (463000 *~ metre) (0 *~ meter)
+      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155) (299.15281) mempty
+      tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) (0) ellipse
+      origin = GridOffset ((-155000)) (463000) (0)
 
 
 -- | Data for the stereographic tests taken from
@@ -303,31 +321,31 @@
 stereographicToGridN :: Bool
 stereographicToGridN = sameGrid g1 g1'
    where
-      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN
-      g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN
+      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridN
+      g1 = GridPoint (196105.283) (557057.739) (0) stereoGridN
       g1' = toGrid stereoGridN p1
 
 stereographicFromGridN :: Bool
 stereographicFromGridN = samePlace p1 p1'
    where
-      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN
-      g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN
+      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridN
+      g1 = GridPoint (196105.283) (557057.739) (0) stereoGridN
       p1' = fromGrid g1
 
 
 stereographicToGridS :: Bool
 stereographicToGridS = sameGrid g1 g1'
    where
-      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridS
-      g1 = GridPoint ((-196105.283) *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridS
+      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridS
+      g1 = GridPoint ((-196105.283)) (557057.739) (0) stereoGridS
       g1' = toGrid stereoGridS p1
 
 
 stereographicFromGridS :: Bool
 stereographicFromGridS = samePlace p1 p1'
    where
-      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridS
-      g1 = GridPoint ((-196105.283) *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridS
+      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridS
+      g1 = GridPoint ((-196105.283)) (557057.739) (0) stereoGridS
       p1' = fromGrid g1
 
 
@@ -345,7 +363,7 @@
 prop_rayPath1 :: Ray WGS84 -> Bool
 prop_rayPath1 r@(Ray pt b e) =
       samePlace pt pt1 && sameAngle b b1 && sameAngle e e1
-   where (pt1,b1,e1) = pathFunc (getRay r) _0
+   where (pt1,b1,e1) = pathFunc (getRay r) 0
 
 
 type ContinuityTest e = Geodetic e -> Bearing -> Azimuth -> Distance -> Distance -> Property
@@ -356,7 +374,7 @@
 -- and have the property that any (point,bearing,azimuth) triple on
 -- the path will specify the same path with a distance offset.
 prop_pathContinuity :: (Ellipsoid e) =>
-   (Geodetic e -> Angle Double -> Angle Double -> Path e) -> ContinuityTest e
+   (Geodetic e -> Double -> Double -> Path e) -> ContinuityTest e
 prop_pathContinuity pf pt0 (Bearing b0) (Azimuth a0) (Distance d1) (Distance d2) =
    counterexample (show ((pt2, Bearing b2, Azimuth a2), (pt3, Bearing b3, Azimuth a3))) $
       pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>
@@ -371,7 +389,7 @@
 
 -- | For continuity testing of ground-based paths (azimuth & altitude always zero)
 -- where lower accuracy is required.
-prop_pathContinuity1 :: (Ellipsoid e) => (Geodetic e -> Angle Double -> Path e) -> ContinuityTest1 e
+prop_pathContinuity1 :: (Ellipsoid e) => (Geodetic e -> Double -> Path e) -> ContinuityTest1 e
 prop_pathContinuity1 pf pt0 (Bearing b0) (Distance2 d1) (Distance2 d2) =
    counterexample (show ((pt2, Bearing b2), (pt3, Bearing b3))) $
       pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>
@@ -393,9 +411,9 @@
 -- This is a test of bisection rather than rays.
 prop_rayBisect :: Ray WGS84 -> Altitude -> Bool
 prop_rayBisect r (Altitude height) =
-   case bisect ray0 f (1 *~ centi meter) (0 *~ meter) (1000 *~ kilo meter) of
+   case bisect ray0 f (1e-2) (0) (1000 * kilometer) of
       Nothing -> False
-      Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1 *~ centi meter
+      Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1e-2
    where
       f g = compare (altitude g) height
       ray0 = getRay r
@@ -409,7 +427,7 @@
 -- | Two rhumb paths intersect at the same place.
 prop_rhumbIntersect :: RhumbPaths2 -> Property
 prop_rhumbIntersect rp =
-   case intersect _0 _0 (10.0 *~ centi meter) 100 path1 path2 of
+   case intersect 0 0 (0.1) 100 path1 path2 of
       Just (d1, d2) ->
          let (pt1, _, _) = pathFunc path1 d1
              (pt2, _, _) = pathFunc path2 d2