diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -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}
diff --git a/etc/LICENCE b/etc/LICENCE
new file mode 100644
--- /dev/null
+++ b/etc/LICENCE
@@ -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.
diff --git a/geodetic.cabal b/geodetic.cabal
new file mode 100644
--- /dev/null
+++ b/geodetic.cabal
@@ -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
+
diff --git a/src/Data/Geo/Geodetic.hs b/src/Data/Geo/Geodetic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Azimuth.hs b/src/Data/Geo/Geodetic/Azimuth.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Azimuth.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Bearing.hs b/src/Data/Geo/Geodetic/Bearing.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Bearing.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Curve.hs b/src/Data/Geo/Geodetic/Curve.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Curve.hs
@@ -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)
diff --git a/src/Data/Geo/Geodetic/Ellipsoid.hs b/src/Data/Geo/Geodetic/Ellipsoid.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Ellipsoid.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/GreatCircle.hs b/src/Data/Geo/Geodetic/GreatCircle.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/GreatCircle.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Haversine.hs b/src/Data/Geo/Geodetic/Haversine.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Haversine.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Sphere.hs b/src/Data/Geo/Geodetic/Sphere.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Sphere.hs
@@ -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
diff --git a/src/Data/Geo/Geodetic/Vincenty.hs b/src/Data/Geo/Geodetic/Vincenty.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Geo/Geodetic/Vincenty.hs
@@ -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
+}
diff --git a/test/doctests.hs b/test/doctests.hs
new file mode 100644
--- /dev/null
+++ b/test/doctests.hs
@@ -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