geodetics 0.1.2 → 1.0.0

raw patch · 13 files changed

+585/−570 lines, 13 filesdep +Streamdep −dimensionaldep −semigroupsdep ~arraydep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: Stream

Dependencies removed: dimensional, semigroups

Dependency ranges changed: array, base

API changes (from Hackage documentation)

Files

• changelog.md +3/−0
• geodetics.cabal +11/−13
• src/Geodetics/Altitude.hs +4/−5
• src/Geodetics/Ellipsoids.hs +90/−72
• src/Geodetics/Geodetic.hs +62/−71
• src/Geodetics/Grid.hs +46/−51
• src/Geodetics/LatLongParser.hs +4/−2
• src/Geodetics/Path.hs +34/−35
• src/Geodetics/Stereographic.hs +36/−39
• src/Geodetics/TransverseMercator.hs +91/−78
• src/Geodetics/UK.hs +49/−54
• test/ArbitraryInstances.hs +70/−83
• test/Main.hs +85/−67

@@ -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).

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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']-   

@@ -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) <*> 

@@ -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