packages feed

geodetic (empty) → 0.1.0

raw patch · 13 files changed

+1214/−0 lines, 13 filesdep +QuickCheckdep +basedep +coordinatebuild-type:Customsetup-changed

Dependencies added: QuickCheck, base, coordinate, directory, doctest, filepath, lens, optional

Files

+ Setup.lhs view
@@ -0,0 +1,44 @@+#!/usr/bin/env runhaskell+\begin{code}+{-# OPTIONS_GHC -Wall #-}+module Main (main) where++import Data.List ( nub )+import Data.Version ( showVersion )+import Distribution.Package ( PackageName(PackageName), PackageId, InstalledPackageId, packageVersion, packageName )+import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) )+import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )+import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose )+import Distribution.Simple.BuildPaths ( autogenModulesDir )+import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), fromFlag )+import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) )+import Distribution.Verbosity ( Verbosity )+import System.FilePath ( (</>) )++main :: IO ()+main = defaultMainWithHooks simpleUserHooks+  { buildHook = \pkg lbi hooks flags -> do+     generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi+     buildHook simpleUserHooks pkg lbi hooks flags+  }++generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO ()+generateBuildModule verbosity pkg lbi = do+  let dir = autogenModulesDir lbi+  createDirectoryIfMissingVerbose verbosity True dir+  withLibLBI pkg lbi $ \_ libcfg -> do+    withTestLBI pkg lbi $ \suite suitecfg -> do+      rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines+        [ "module Build_" ++ testName suite ++ " where"+        , "deps :: [String]"+        , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg))+        ]+  where+    formatdeps = map (formatone . snd)+    formatone p = case packageName p of+      PackageName n -> n ++ "-" ++ showVersion (packageVersion p)++testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]+testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys++\end{code}
+ etc/LICENCE view
@@ -0,0 +1,27 @@+Copyright 2013,2014 Tony Morris++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ geodetic.cabal view
@@ -0,0 +1,76 @@+name:               geodetic+version:            0.1.0+license:            BSD3+license-File:       etc/LICENCE+author:             Tony Morris <ʇǝu˙sıɹɹoɯʇ@ןןǝʞsɐɥ>+maintainer:         Tony Morris+copyright:          Copyright (C) 2013,2014 Tony Morris+synopsis:           Geodetic calculations+category:           Development+description:        Geodetic calculations including Vincenty and Great Circle using a Latitude and Longitude pair+homepage:           https://github.com/tonymorris/geodetic+bug-reports:        https://github.com/tonymorris/geodetic/issues+cabal-version:      >= 1.10+build-type:         Custom++source-repository   head+  type:             git+  location:         git@github.com:tonymorris/geodetic.git++flag                small_base+  description:      Choose the new, split-up base package.++library+  default-language:+                    Haskell2010++  build-depends:+                    base < 5 && >= 3+                    , lens >= 3.10+                    , coordinate >= 0.0.8+                    , optional >= 0.0.1++  ghc-options:+                    -Wall++  default-extensions:+                      NoImplicitPrelude++  hs-source-dirs:+                    src++  exposed-modules:+                    Data.Geo.Geodetic+                    Data.Geo.Geodetic.Azimuth+                    Data.Geo.Geodetic.Bearing+                    Data.Geo.Geodetic.Curve+                    Data.Geo.Geodetic.Ellipsoid+                    Data.Geo.Geodetic.GreatCircle+                    Data.Geo.Geodetic.Haversine+                    Data.Geo.Geodetic.Sphere+                    Data.Geo.Geodetic.Vincenty++test-suite doctests+  type:+                    exitcode-stdio-1.0++  main-is:+                    doctests.hs++  default-language:+                    Haskell2010++  build-depends:+                    base < 5 && >= 3,+                    doctest >= 0.9.7,+                    filepath >= 1.3,+                    directory >= 1.1,+                    QuickCheck >= 2.0++  ghc-options:+                    -Wall+                    -threaded++  hs-source-dirs:+                    test+
+ src/Data/Geo/Geodetic.hs view
@@ -0,0 +1,10 @@+module Data.Geo.Geodetic(module G) where++import Data.Geo.Geodetic.Azimuth as G+import Data.Geo.Geodetic.Bearing as G+import Data.Geo.Geodetic.Curve as G+import Data.Geo.Geodetic.Ellipsoid as G+import Data.Geo.Geodetic.GreatCircle as G+import Data.Geo.Geodetic.Haversine as G+import Data.Geo.Geodetic.Sphere as G+import Data.Geo.Geodetic.Vincenty as G
+ src/Data/Geo/Geodetic/Azimuth.hs view
@@ -0,0 +1,83 @@+-- | An azimuth in degrees between 0 and 360.+module Data.Geo.Geodetic.Azimuth(+  Azimuth+, HasAzimuth(..)+, modAzimuth+, nAzimuth+) where++import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), showParen, showString)+import Data.Maybe(Maybe(..))+import Control.Lens(Prism', Lens', prism')+import Text.Printf(printf)+import Data.Fixed(mod')++-- $setup+-- >>> import Control.Lens((#), (^?))+-- >>> import Data.Foldable(all)+-- >>> import Prelude(Eq(..))++newtype Azimuth =+  Azimuth Double+  deriving (Eq, Ord)++-- | A show instance that prints to 4 decimal places.+-- This is to take floating-point rounding errors into account.+instance Show Azimuth where+  showsPrec n (Azimuth d) =+    showParen (n > 10) (showString ("Azimuth " ++ printf "%0.4f" d))++-- | Construct an azimuth such that if the given value is out of bounds,+-- a modulus is taken to keep it within 0 inclusive and 360 exclusive.+--+-- >>> modAzimuth 7+-- Azimuth 7.0000+--+-- >>> modAzimuth 0+-- Azimuth 0.0000+--+-- >>> modAzimuth 360+-- Azimuth 0.0000+--+-- >>> modAzimuth 361+-- Azimuth 1.0000+--+-- >>> modAzimuth 359.999+-- Azimuth 359.9990+modAzimuth ::+  Double+  -> Azimuth+modAzimuth x =+  Azimuth (x `mod'` 360)++-- | A prism on azimuth to an integer between 0 and 359 inclusive.+--+-- >>> 7 ^? nAzimuth+-- Just (Azimuth 7.0000)+--+-- >>> 0 ^? nAzimuth+-- Just (Azimuth 0.0000)+--+-- >>> 359.999 ^? nAzimuth+-- Just (Azimuth 359.9990)+--+-- >>> 360 ^? nAzimuth+-- Nothing+--+-- prop> all (\m -> nAzimuth # m == n) (n ^? nAzimuth)+nAzimuth ::+  Prism' Double Azimuth+nAzimuth =+  prism'+    (\(Azimuth i) -> i)+    (\i -> case i >= 0 && i < 360 of+             True -> Just (Azimuth i)+             False -> Nothing)++class HasAzimuth t where+  azimuth ::+    Lens' t Azimuth++instance HasAzimuth Azimuth where+  azimuth =+    id
+ src/Data/Geo/Geodetic/Bearing.hs view
@@ -0,0 +1,117 @@+-- | A bearing in degrees between 0 and 360.+module Data.Geo.Geodetic.Bearing(+  Bearing+, HasBearing(..)+, modBearing+, degreeBearing+, radianBearing+) where++import Prelude(Double, Bool(..), Eq, Show(..), Num(..), Fractional(..), Ord(..), id, (&&), (++), (.), showString, showParen, pi)+import Data.Maybe(Maybe(..))+import Control.Lens(Prism', Lens', prism', iso)+import Text.Printf(printf)+import Data.Fixed(mod')++-- $setup+-- >>> import Control.Lens((#), (^?))+-- >>> import Data.Foldable(all)+-- >>> import Prelude(Eq(..))++newtype Bearing =+  Bearing Double+  deriving (Eq, Ord)++-- | A show instance that prints to 4 decimal places.+-- This is to take floating-point rounding errors into account.+instance Show Bearing where+  showsPrec n (Bearing d) =+    showParen (n > 10) (showString ("Bearing " ++ printf "%0.4f" d))++-- | Construct a bearing such that if the given value is out of bounds,+-- a modulus is taken to keep it within 0 inclusive and 360 exclusive.+--+-- >>> modBearing 7+-- Bearing 7.0000+--+-- >>> modBearing 0+-- Bearing 0.0000+--+-- >>> modBearing (-0.0001)+-- Bearing 359.9999+--+-- >>> modBearing 360+-- Bearing 0.0000+--+-- >>> modBearing 359.99999+-- Bearing 360.0000+--+-- >>> modBearing 359.999+-- Bearing 359.9990+modBearing ::+  Double+  -> Bearing+modBearing x =+  Bearing (x `mod'` 360)++-- | A prism on bearing to a double between 0 inclusive and 360 exclusive.+--+-- >>> 7 ^? degreeBearing+-- Just (Bearing 7.0000)+--+-- >>> 0 ^? degreeBearing+-- Just (Bearing 0.0000)+--+-- >>> 359 ^? degreeBearing+-- Just (Bearing 359.0000)+--+-- >>> 359.997 ^? degreeBearing+-- Just (Bearing 359.9970)+--+-- >>> 360 ^? degreeBearing+-- Nothing+--+-- prop> all (\m -> degreeBearing # m == n) (n ^? degreeBearing)+degreeBearing ::+  Prism' Double Bearing+degreeBearing =+  prism'+    (\(Bearing i) -> i)+    (\i -> case i >= 0 && i < 360 of+             True -> Just (Bearing i)+             False -> Nothing)++-- | A prism on bearing to a double between 0 and π exclusive.+--+-- >>> (2 * pi - 0.0000000001) ^? radianBearing+-- Just (Bearing 360.0000)+--+-- >>> 0 ^? radianBearing+-- Just (Bearing 0.0000)+--+-- >>> 0.001 ^? radianBearing+-- Just (Bearing 0.0573)+--+-- >>> 1.78391 ^? radianBearing+-- Just (Bearing 102.2105)+--+-- >>> pi ^? radianBearing+-- Just (Bearing 180.0000)+--+-- >>> (2 * pi) ^? radianBearing+-- Nothing+--+-- >>> (-0.001) ^? radianBearing+-- Nothing+radianBearing ::+  Prism' Double Bearing+radianBearing =+  iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . degreeBearing++class HasBearing t where+  bearing ::+    Lens' t Bearing++instance HasBearing Bearing where+  bearing =+    id
+ src/Data/Geo/Geodetic/Curve.hs view
@@ -0,0 +1,51 @@+-- | A geodetic curve is made of a distance in metres, an azimuth and a reverse azimuth.+module Data.Geo.Geodetic.Curve(+  Curve+, curve+, curveDistance+, curveAzimuth+, curveReverseAzimuth+) where++import Prelude(Eq, Show(..), Ord(..), Double, (.), abs, showString, showParen)+import Data.List(unwords)+import Text.Printf(printf)+import Control.Lens+import Data.Geo.Geodetic.Azimuth++data Curve =+  Curve+    Double -- The ellipsoidal distance.+    Azimuth -- The azimuth.+    Azimuth -- The reverse azimuth.+  deriving (Eq, Ord)++-- | A show instance that prints to 4 decimal places.+-- This is to take floating-point rounding errors into account.+instance Show Curve where+  showsPrec n (Curve d a r) =+    showParen (n > 10) (showString (unwords ["GeodeticCurve", printf "%0.4f" d, show a, show r]))++-- | Construct a geodetic curve with the given parameters.+curve ::+  Double -- ^ The ellipsoidal distance.+  -> Azimuth -- ^ The azimuth.+  -> Azimuth -- ^ The reverse azimuth.+  -> Curve+curve =+  Curve++curveDistance ::+  Lens' Curve Double+curveDistance =+  lens (\(Curve d _ _) -> d) (\(Curve _ a r) d -> Curve d a r)++curveAzimuth ::+  Lens' Curve Azimuth+curveAzimuth =+  lens (\(Curve _ a _) -> a) (\(Curve d _ r) a -> Curve d a r)++curveReverseAzimuth ::+  Lens' Curve Azimuth+curveReverseAzimuth =+  lens (\(Curve _ _ r) -> r) (\(Curve d a _) r -> Curve d a r)
+ src/Data/Geo/Geodetic/Ellipsoid.hs view
@@ -0,0 +1,199 @@+module Data.Geo.Geodetic.Ellipsoid(+  -- * Data type+  Ellipsoid+, HasEllipsoid(..)+  -- * Ellipsoid properties+, HasSemiMajor(..)+, HasSemiMinor(..)+, HasFlattening(..)+, HasInverseFlattening(..)+  -- * Ellipsoid construction+, semiMajorFlattening+, semiMinorFlattening+, semiMajorInverseFlattening+, semiMinorInverseFlattening+  -- * Ellipsoids+, wgs84+, grs80+, grs67+, ans+, wgs72+, au1965+, krasovsky1940+, international1924+, hayford1909+, airy1830+, everest1830+, bessel1841+, clarke1858+, clarke1866+, clarke1880+) where++import Prelude(Eq, Ord, Show, Num(..), Fractional(..), Double, id)+import Control.Lens(Lens', lens)++data Ellipsoid =+  Ellipsoid+    Double -- The semi major axis in metres.+    Double -- The ellipsoidal flattening.+  deriving (Eq, Ord, Show)++semiMajorFlattening ::+  Double -- ^ The semi major axis in metres.+  -> Double -- ^ The ellipsoidal flattening.+  -> Ellipsoid+semiMajorFlattening =+  Ellipsoid++semiMinorFlattening ::+  Double -- ^ The semi minor axis in metres.+  -> Double -- ^ The ellipsoidal flattening.+  -> Ellipsoid+semiMinorFlattening s f =+  Ellipsoid (s * f - 1.0) f++semiMajorInverseFlattening ::+  Double -- ^ The semi major axis in metres.+  -> Double -- ^ The ellipsoidal inverse flattening.+  -> Ellipsoid+semiMajorInverseFlattening s f =+  semiMajorFlattening s (1/f)++semiMinorInverseFlattening ::+  Double -- ^ The semi minor axis in metres.+  -> Double -- ^ The ellipsoidal inverse flattening.+  -> Ellipsoid+semiMinorInverseFlattening s f =+  semiMinorFlattening s (1/f)++class HasEllipsoid t where+  ellipsoid ::+    Lens' t Ellipsoid++instance HasEllipsoid Ellipsoid where+  ellipsoid =+    id++class HasSemiMajor t where+  semiMajor ::+    Lens' t Double++instance HasSemiMajor Double where+  semiMajor =+    id++instance HasSemiMajor Ellipsoid where+  semiMajor =+    lens (\(Ellipsoid s _) -> s) (\(Ellipsoid _ f) s -> Ellipsoid s f)++class HasSemiMinor t where+  semiMinor ::+    Lens' t Double++instance HasSemiMinor Double where+  semiMinor =+    id++instance HasSemiMinor Ellipsoid where+  semiMinor =+    lens (\(Ellipsoid s f) -> (1.0 - f) * s) (\(Ellipsoid _ f) s -> Ellipsoid (s * f - 1.0) f)++class HasFlattening t where+  flattening ::+    Lens' t Double++instance HasFlattening Double where+  flattening =+    id++instance HasFlattening Ellipsoid where+  flattening =+    lens (\(Ellipsoid _ f) -> f) (\(Ellipsoid s _) f -> Ellipsoid s f)++class HasInverseFlattening t where+  inverseFlattening ::+    Lens' t Double++instance HasInverseFlattening Double where+  inverseFlattening =+    id++instance HasInverseFlattening Ellipsoid where+  inverseFlattening =+    lens (\(Ellipsoid _ f) -> 1/f) (\(Ellipsoid s _) f -> Ellipsoid s (1/f))++wgs84 ::+  Ellipsoid+wgs84 =+  semiMajorInverseFlattening 6378137 298.257223563++grs80 ::+  Ellipsoid+grs80 =+  semiMajorInverseFlattening 6378137 298.257222101++grs67 ::+  Ellipsoid+grs67 =+  semiMajorInverseFlattening 6378160 298.25++ans ::+  Ellipsoid+ans =+  semiMajorInverseFlattening 6378160 298.25++wgs72 ::+  Ellipsoid+wgs72 =+  semiMajorInverseFlattening 6378135 298.26++au1965 ::+  Ellipsoid+au1965 =+  semiMajorInverseFlattening 6378160 298.25++krasovsky1940 ::+  Ellipsoid+krasovsky1940 =+  semiMajorInverseFlattening 6378245 298.3++international1924 ::+  Ellipsoid+international1924 =+  semiMajorInverseFlattening 6378388 297++hayford1909 ::+  Ellipsoid+hayford1909 =+  international1924++airy1830 ::+  Ellipsoid+airy1830 =+  semiMajorInverseFlattening 6377563.4 299.32++everest1830 ::+  Ellipsoid+everest1830 =+  semiMajorInverseFlattening 6377276.3 300.8++bessel1841 ::+  Ellipsoid+bessel1841 =+  semiMajorInverseFlattening 6377397.2 299.15++clarke1858 ::+  Ellipsoid+clarke1858 =+  semiMajorInverseFlattening 6378293.645 294.26++clarke1866 ::+  Ellipsoid+clarke1866 =+  semiMajorInverseFlattening 6378206.4 294.98++clarke1880 ::+  Ellipsoid+clarke1880 =+  semiMajorInverseFlattening 6378249.145 293.465
+ src/Data/Geo/Geodetic/GreatCircle.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE FlexibleContexts #-}++-- | Great circle geodetic distance algorithm.+module Data.Geo.Geodetic.GreatCircle(+  sphericalLaw+, sphericalLawD+, sphericalLaw'+) where++import Prelude(Double, Num(..), Fractional(..), pi, sin, cos, acos)+import Control.Lens((#), (^.))+import System.Args.Optional(Optional1(..))+import Data.Geo.Coordinate+import Data.Geo.Geodetic.Sphere++-- $setup+-- >>> import Prelude(Functor(..), Monad(..), String)+-- >>> import Data.Maybe(Maybe)+-- >>> import Text.Printf(printf)++-- | Great circle spherical law algorithm.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw earthMean fr to)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw earthMean fr to)) :: Maybe String+-- Just "17128743.0669"+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw (6350000 ^. nSphere) fr to)) :: Maybe String+-- Just "14959840.4461"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw (6350000 ^. nSphere) fr to)) :: Maybe String+-- Just "17081801.7377"+sphericalLaw ::+  Sphere+  -> Coordinate+  -> Coordinate+  -> Double+sphericalLaw s start end =+  let toRadians n = n * pi / 180+      lat1 = toRadians (fracLatitude # (start ^. latitude))+      lat2 = toRadians (fracLatitude # (end ^. latitude))+      lon1 = toRadians (fracLongitude # (start ^. longitude))+      lon2 = toRadians (fracLongitude # (end ^. longitude))+  in acos (sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)) * nSphere # s++-- | Great circle spherical law algorithm with a default sphere of the earth mean.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLawD fr to)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLawD fr to)) :: Maybe String+-- Just "17128743.0669"+sphericalLawD ::+  Coordinate+  -> Coordinate+  -> Double+sphericalLawD =+  sphericalLaw earthMean++-- | Great circle spherical law algorithm with an optionally applied default sphere of the earth mean.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (sphericalLaw' fr to :: Double)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw' fr to :: Double)) :: Maybe String+-- Just "17128743.0669"+sphericalLaw' ::+  (Optional1+    Sphere+    (+      Coordinate+      -> Coordinate+      -> Double+    ) x) =>+    x+sphericalLaw' =+  optional1 sphericalLaw earthMean
+ src/Data/Geo/Geodetic/Haversine.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE FlexibleContexts #-}++-- | Haversine geodetic distance algorithm.+module Data.Geo.Geodetic.Haversine(+  haversine+, haversineD+, haversine'+) where++import Prelude(Double, Num(..), Fractional(..), (.), pi, atan, sin, atan2, cos, sqrt)+import Control.Lens((#), (^.))+import System.Args.Optional(Optional1(..))+import Data.Geo.Coordinate+import Data.Geo.Geodetic.Sphere++-- $setup+-- >>> import Prelude(Functor(..), Monad(..), String)+-- >>> import Data.Maybe(Maybe)+-- >>> import Text.Printf(printf)++-- | Haversine algorithm.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine earthMean fr to)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine earthMean fr to)) :: Maybe String+-- Just "17128743.0669"+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine (6350000 ^. nSphere) fr to)) :: Maybe String+-- Just "14959840.4461"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine (6350000 ^. nSphere) fr to)) :: Maybe String+-- Just "17081801.7377"+haversine ::+  Sphere+  -> Coordinate+  -> Coordinate+  -> Double+haversine s start end =+  let lat1 = fracLatitude # (start ^. latitude)+      lat2 = fracLatitude # (end ^. latitude)+      toRadians n = n * pi / 180+      dlat = (toRadians (lat1 - lat2)) / 2+      dlon = (toRadians (fracLongitude # (start ^. longitude) - fracLongitude # (end ^. longitude))) / 2+      cosr = cos . toRadians+      square x = x * x+      a = square (sin dlat) + cosr lat1 * cosr lat2 * square (sin (dlon))+      c = 2 * atan2 (sqrt a) (sqrt (1 - a))+  in (nSphere # s) * c++-- | Haversine algorithm with a default sphere of the earth mean.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversineD fr to)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversineD fr to)) :: Maybe String+-- Just "17128743.0669"+haversineD ::+  Coordinate+  -> Coordinate+  -> Double+haversineD =+  haversine earthMean++-- | Haversine algorithm with an optionally applied default sphere of the earth mean.+--+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (haversine' fr to :: Double)) :: Maybe String+-- Just "15000950.5589"+--+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine' fr to :: Double)) :: Maybe String+-- Just "17128743.0669"+haversine' ::+  (Optional1+    Sphere+    (+      Coordinate+      -> Coordinate+      -> Double+    ) x) =>+    x+haversine' =+  optional1 haversine earthMean
+ src/Data/Geo/Geodetic/Sphere.hs view
@@ -0,0 +1,54 @@+-- | A sphere with a radius in metres.+module Data.Geo.Geodetic.Sphere(+  Sphere+, HasSphere(..)+, nSphere+, earthMean+) where++import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), (.), abs, showParen, showString)+import Control.Lens(Iso', Lens', iso)+import Text.Printf(printf)++-- $setup+-- >>> import Control.Lens((#), (^.))+-- >>> import Data.Foldable(all)+-- >>> import Prelude(Eq(..))++newtype Sphere =+  Sphere Double+  deriving (Eq, Ord)++-- | A show instance that prints to 4 decimal places.+-- This is to take floating-point rounding errors into account.+instance Show Sphere where+  showsPrec n (Sphere d) =+    showParen (n > 10) (showString ("Sphere " ++ printf "%0.4f" d))++-- | An isomorphism on sphere to a double.+--+-- >>> 7 ^. nSphere+-- Sphere 7.0000+--+-- >>> 0 ^. nSphere+-- Sphere 0.0000+--+-- >>> (-7) ^. nSphere+-- Sphere -7.0000+nSphere ::+  Iso' Double Sphere+nSphere =+  iso Sphere (\(Sphere d) -> d)++earthMean ::+  Sphere+earthMean =+  Sphere 6367450++class HasSphere t where+  sphere ::+    Lens' t Sphere++instance HasSphere Sphere where+  sphere =+    id
+ src/Data/Geo/Geodetic/Vincenty.hs view
@@ -0,0 +1,361 @@+{-# LANGUAGE FlexibleContexts #-}++-- | An implementation of Thaddeus Vincenty's direct and inverse geodetic algorithms. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>+module Data.Geo.Geodetic.Vincenty(+  Convergence+, convergence+, direct+, directD+, direct'+, VincentyDirectResult+, inverse+, inverseD+, inverse'+) where++import Prelude(Eq(..), Show(..), Ord(..), Num(..), Floating(..), Fractional(..), Double, Int, Bool, Ordering(..), subtract, cos, sin, asin, tan, sqrt, atan, atan2, pi, (.), (++), (&&), ($!), error)+import Control.Lens(lens, (^.), (#), (^?))+import Data.Maybe(fromMaybe)+import System.Args.Optional(Optional2(..))+import Data.Geo.Coordinate+import Data.Geo.Geodetic.Azimuth+import Data.Geo.Geodetic.Bearing+import Data.Geo.Geodetic.Ellipsoid+import Data.Geo.Geodetic.Curve++-- $setup+-- >>> import Prelude++type Convergence =+  Double++-- | A typically acceptable convergence value.+convergence ::+  Convergence+convergence =+  0.0000000000001++data VincentyDirectResult =+  VincentyDirectResult+    Coordinate+    Bearing+  deriving (Eq, Ord, Show)++instance HasCoordinate VincentyDirectResult where+  coordinate =+    lens (\(VincentyDirectResult c _) -> c) (\(VincentyDirectResult _ b) c -> VincentyDirectResult c b)++instance HasBearing VincentyDirectResult where+  bearing =+    lens (\(VincentyDirectResult _ b) -> b) (\(VincentyDirectResult c _) b -> VincentyDirectResult c b)++-- | Vincenty direct algorithm.+--+-- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))+--+-- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))+--+-- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0986)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1464))) (Bearing 165.3451))+--+-- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0662)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4139))) (Bearing 165.3183))+direct ::+  Ellipsoid+  -> Convergence+  -> Coordinate+  -> Bearing+  -> Double+  -> VincentyDirectResult+direct e conv start bear dist =+  let sMnr = e ^. semiMinor+      flat = e ^. flattening+      alpha = radianBearing # bear+      cosAlpha = cos alpha+      sinAlpha = sin alpha+      tanu1 = (1.0 - flat) * tan (radianLatitude # (start ^. latitude))+      cosu1 = 1.0 / sqrt (1.0 + square tanu1)+      sinu1 = tanu1 * cosu1+      sigma1 = atan2 tanu1 cosAlpha+      csa = cosu1 * sinAlpha+      sin2Alpha = square csa+      cos2Alpha = 1 - sin2Alpha+      ab d f g h i = let s = cos2Alpha * (square (e ^. semiMajor / sMnr) - 1)+                     in (s / d) * (f + s * (g + s * (h - i * s)))+      a = 1 + ab 16384 4096 (-768) 320 175+      b = ab 1024 256 (-128) 74 47+      end = let begin = ps (dist / sMnr / a)+                iter p = let tf d = -3 + 4 * d+                             cosSigma'' = cosSigma' p+                             sinSigma'' = sinSigma' p+                             cosSigmaM2'' = cosSigmaM2' sigma1 p+                             cos2SigmaM2'' = cos2SigmaM2' sigma1 p+                             deltaSigma = b * sinSigma'' * (cosSigmaM2'' + b / 4.0 * (cosSigma'' * (-1 + 2 * cos2SigmaM2'') - (b / 6.0) * cosSigmaM2'' * tf (square sinSigma'') * tf cos2SigmaM2''))+                         in transition p deltaSigma+                pred' p = abs (sigma' p - prevSigma' p) >= conv+            in doWhile iter pred' begin+      sigma'' = sigma' end+      sinSigma = sinSigma' end+      cosSigmaM2 = cosSigmaM2' sigma1 end+      cos2SigmaM2 = cos2SigmaM2' sigma1 end+      cosSigma = cos sigma''+      c = flat / 16 * cos2Alpha * (4 + flat * (4 - 3 * cos2Alpha))+      cc = cosu1 * cosSigma+      ccca = cc * cosAlpha+      sss = sinu1 * sinSigma+      latitude' = let r = atan2 (sinu1 * cosSigma + cosu1 * sinSigma * cosAlpha) ((1.0 - flat) * sqrt (sin2Alpha + (sss - ccca) ** 2.0))+                  in fromMaybe (error ("Invariant not met. Latitude in radians not within range " ++ show r)) (r ^? radianLatitude)+      longitude' = let r = fracLongitude # (start ^. longitude) + ((atan2 (sinSigma * sinAlpha) (cc - sss * cosAlpha) - (1 - c) * flat * csa * (sigma'' + c * sinSigma * (cosSigmaM2 + c * cosSigma * (-1 + 2 * cos2SigmaM2)))) * 180 / pi)+                    in fromMaybe (error ("Invariant not met. Longitude in radians not within range " ++ show r)) (r ^? fracLongitude)+  in VincentyDirectResult+       (latitude' .#. longitude')+       (+         let r = atan2 csa (ccca - sss)+         in fromMaybe (error ("Invariant not met. Bearing in radians not within range " ++ show r)) (r ^? radianBearing)+       )++-- | Vincenty direct algorithm with a default ellipsoid of WGS84 and standard convergence.+--+-- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))+--+-- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))+directD ::+  Coordinate+  -> Bearing+  -> Double+  -> VincentyDirectResult+directD =+  direct wgs84 convergence++-- | Vincenty direct algorithm with an optionally applied default ellipsoid of WGS84 and standard convergence.+--+-- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) (27.812 ..#.. 154.295)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude 27) (Minutes 46) (Seconds 30.0981)) (Longitude (DegreesLongitude 154) (Minutes 18) (Seconds 21.1466))) (Bearing 165.3451))+--+-- >>> fmap (\c' -> direct' c' (modBearing 165.34) (4235 :: Double) :: VincentyDirectResult) ((-66.093) ..#.. 12.84)+-- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183))+direct' ::+  (Optional2+    Ellipsoid+    Convergence+    (+      Coordinate+      -> Bearing+      -> Double+      -> VincentyDirectResult+    ) x) =>+    x+direct' =+  optional2 direct wgs84 convergence++-- | Vincenty inverse algorithm.+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse wgs84 convergence fr to)+-- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse wgs84 convergence fr to)+-- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse ans convergence fr to)+-- Just (GeodeticCurve 14998630.4056 Azimuth 180.0000 Azimuth 0.0000)+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse ans convergence fr to)+-- Just (GeodeticCurve 7099229.9126 Azimuth 0.0000 Azimuth 180.0000)+inverse ::+  Ellipsoid+  -> Convergence+  -> Coordinate+  -> Coordinate+  -> Curve+inverse e conv start end =+  let b = e ^. semiMinor+      f = e ^. flattening+      (phi1, phi2) =+        let rl k = radianLatitude # (k ^. latitude)+        in (rl start, rl end)+      a2b2b2 =+        let ss z = square (z e)+        in ss (^. semiMajor) / ss (^. semiMinor) - 1+      omega =+        let rl k = radianLongitude # (k ^. longitude)+        in rl end - rl start+      (u1, u2) =+        let at = atan . ((1 - f) *) . tan+        in (at phi1, at phi2)+      sinu1 = sin u1+      cosu1 = cos u1+      sinu2 = sin u2+      cosu2 = cos u2+      sinu1sinu2 = sinu1 * sinu2+      cosu1sinu2 = cosu1 * sinu2+      sinu1cosu2 = sinu1 * cosu2+      cosu1cosu2 = cosu1 * cosu2+      begin = Q 0 Continue omega 0 0 0+      iter q = let sinlambda = sin (lambda q)+                   coslambda = cos (lambda q)+                   sin2sigma = square cosu2 * square sinlambda + square (cosu1sinu2 - sinu1cosu2 * coslambda)+                   sinsigma = sqrt sin2sigma+                   cossigma = sinu1sinu2 + cosu1cosu2 * coslambda+                   sigma'' = atan2 sinsigma cossigma+                   sinalpha = if sin2sigma == 0.0 then 0.0 else cosu1cosu2 * sinlambda / sinsigma+                   alpha = asin sinalpha+                   cos2alpha = square (cos alpha)+                   cos2sigmam = if cos2alpha == 0.0 then 0.0 else cossigma - 2 * sinu1sinu2 / cos2alpha+                   u2' = cos2alpha * a2b2b2+                   cos2sigmam2 = square cos2sigmam+                   a = 1.0 + u2' / 16384 * (4096 + u2' * (u2' * (320 - 175 * u2') - 768))+                   b' = u2' / 1024 * (256 + u2' * (u2' * (74 - 47 * u2') - 128))+                   deltasigma' = b' * sinsigma * (cos2sigmam + b' / 4 * (cossigma * (2 * cos2sigmam2 - 1) - b' / 6 * cos2sigmam * (4 * sin2sigma - 3) * (cos2sigmam2 * 4 - 3)))+                   c' = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha))+                   l = omega + (1 - c') * f * sinalpha * (sigma'' + c' * sinsigma * (cos2sigmam + c' * cossigma * (2 * cos2sigmam2 - 1)))+                   r = let c = count q+                       in if c == 20+                             then Limit+                             else if c > 1 && cos alpha < conv+                                   then Converge+                                   else Continue+              in Q (count q + 1) r l a sigma'' deltasigma'+      pred' = (== Continue) . result+      ed = whileDo iter pred' begin+      ifi p t a = if p a then t a else a+      (alpha1, alpha2) =+        let alphaNoConverge :: (Fractional b0, Num b0, Ord b0) => Bool -> Ordering -> b0 -> b0 -> (b0, b0)+            alphaNoConverge c cp x y = vmap2 (ifi (>= 360) (subtract 360)) (if c+                                                                              then (x, y)+                                                                              else if cp == GT+                                                                                then (180.0, 0.0)+                                                                                else if cp == LT+                                                                                  then (0.0, 180.0)+                                                                                  else let nan = 0/0+                                                                                        in (nan, nan))+        in alphaNoConverge (result ed == Converge) (compare phi1 phi2) 0 0+  in curve (b * a' ed * (sigma ed - deltasigma ed)) (modAzimuth alpha1) (modAzimuth alpha2)++-- | Vincenty inverse algorithm with a default ellipsoid of WGS84 and standard convergence.+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverseD fr to)+-- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverseD fr to)+-- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)+inverseD ::+  Coordinate+  -> Coordinate+  -> Curve+inverseD =+  inverse wgs84 convergence++-- | Vincenty inverse algorithm with an optionally applied default ellipsoid of WGS84 and standard convergence.+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse' fr to :: Curve)+-- Just (GeodeticCurve 14998576.9860 Azimuth 180.0000 Azimuth 0.0000)+--+-- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse' fr to :: Curve)+-- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000)+inverse' ::+  (Optional2+    Ellipsoid+    Convergence+    (+      Coordinate+      -> Coordinate+      -> Curve+    ) x) =>+    x+inverse' =+  optional2 inverse wgs84 convergence++---- not exported++data P = P {+  origSigma' :: Double+, sigma' :: Double+, prevSigma' :: Double+} deriving Show++vmap2 ::+  (a -> b)+  -> (a, a)+  -> (b, b)+vmap2 f (a1, a2) =+  (f a1, f a2)++ps ::+  Double+  -> P+ps s =+  P s s s++transition ::+  P+  -> Double+  -> P+transition p d =+  P (origSigma' p) (d + origSigma' p) (sigma' p)++sinSigma' ::+  P+  -> Double+sinSigma' =+  sin . sigma'++cosSigma' ::+  P+  -> Double+cosSigma' =+  cos . sigma'++sigmaM2' ::+  Double+  -> P+  -> Double+sigmaM2' s p =+  2.0 * s + sigma' p++cosSigmaM2' ::+  Double+  -> P+  -> Double+cosSigmaM2' s p =+  cos (sigmaM2' s p)++cos2SigmaM2' ::+  Double+  -> P+  -> Double+cos2SigmaM2' s p =+  square (cosSigmaM2' s p)++square ::+  Num a =>+  a+  -> a+square a =+  a * a++doWhile ::+  (a -> a)+  -> (a -> Bool)+  -> a+  -> a+doWhile f p a =+  let x = f a+  in if p x then doWhile f p x else x++whileDo :: (a -> a) -> (a -> Bool) -> a -> a+whileDo f p a = if p a then whileDo f p $! (f a) else a++data InverseResult = Continue | Limit | Converge deriving Eq++data Q = Q {+  count :: Int,+  result :: InverseResult,+  lambda :: Double,+  a' :: Double,+  sigma :: Double,+  deltasigma :: Double+}
+ test/doctests.hs view
@@ -0,0 +1,32 @@+module Main where++import Build_doctests (deps)+import Control.Applicative+import Control.Monad+import Data.List+import System.Directory+import System.FilePath+import Test.DocTest++main ::+  IO ()+main =+  getSources >>= \sources -> doctest $+      "-isrc"+    : "-idist/build/autogen"+    : "-optP-include"+    : "-optPdist/build/autogen/cabal_macros.h"+    : "-hide-all-packages"+    : map ("-package="++) deps ++ sources++getSources :: IO [FilePath]+getSources = filter (isSuffixOf ".hs") <$> go "src"+  where+    go dir = do+      (dirs, files) <- getFilesAndDirectories dir+      (files ++) . concat <$> mapM go dirs++getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])+getFilesAndDirectories dir = do+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c