diff --git a/Data/Geo.hs b/Data/Geo.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo.hs
@@ -0,0 +1,51 @@
+module Data.Geo(
+                 module Data.Geo.Azimuth,
+                 module Data.Geo.Bearing,
+                 module Data.Geo.Coord,
+                 module Data.Geo.DMS,
+                 module Data.Geo.ElevatedCurve,
+                 module Data.Geo.Elevation,
+                 module Data.Geo.Ellipsoid,
+                 module Data.Geo.GeodeticCurve,
+                 module Data.Geo.GreatCircle,
+                 module Data.Geo.Haversine,
+                 module Data.Geo.Latitude,
+                 module Data.Geo.Longitude,
+                 module Data.Geo.Position,
+                 module Data.Geo.Radians,
+                 module Data.Geo.Sphere,
+                 module Data.Geo.Vincenty,
+                 module Data.Geo.Accessor.Azi,
+                 module Data.Geo.Accessor.Coordinate,
+                 module Data.Geo.Accessor.Curve,
+                 module Data.Geo.Accessor.Ele,
+                 module Data.Geo.Accessor.EllipsoidalDistance,
+                 module Data.Geo.Accessor.Lat,
+                 module Data.Geo.Accessor.Lon,
+                 module Data.Geo.Accessor.ReverseAzi
+               ) where
+
+import Data.Geo.Azimuth
+import Data.Geo.Bearing
+import Data.Geo.Coord
+import Data.Geo.DMS
+import Data.Geo.ElevatedCurve
+import Data.Geo.Elevation
+import Data.Geo.Ellipsoid
+import Data.Geo.GeodeticCurve
+import Data.Geo.GreatCircle
+import Data.Geo.Haversine
+import Data.Geo.Latitude
+import Data.Geo.Longitude
+import Data.Geo.Position
+import Data.Geo.Radians
+import Data.Geo.Sphere
+import Data.Geo.Vincenty
+import Data.Geo.Accessor.Azi
+import Data.Geo.Accessor.Coordinate
+import Data.Geo.Accessor.Curve
+import Data.Geo.Accessor.Ele
+import Data.Geo.Accessor.EllipsoidalDistance
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Data.Geo.Accessor.ReverseAzi
diff --git a/Data/Geo/Accessor/Azi.hs b/Data/Geo/Accessor/Azi.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Azi.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Azi where
+
+import Data.Geo.Azimuth
+
+class Azi a where
+  azi :: a -> Azimuth
diff --git a/Data/Geo/Accessor/Coordinate.hs b/Data/Geo/Accessor/Coordinate.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Coordinate.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Coordinate where
+
+import Data.Geo.Coord
+
+class Coordinate a where
+  coordinate :: a -> Coord
diff --git a/Data/Geo/Accessor/Curve.hs b/Data/Geo/Accessor/Curve.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Curve.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Curve where
+
+import Data.Geo.GeodeticCurve
+
+class Curve a where
+  curve :: a -> GeodeticCurve
diff --git a/Data/Geo/Accessor/Ele.hs b/Data/Geo/Accessor/Ele.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Ele.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Ele where
+
+import Data.Geo.Elevation
+
+class Ele a where
+  ele :: a -> Elevation
diff --git a/Data/Geo/Accessor/EllipsoidalDistance.hs b/Data/Geo/Accessor/EllipsoidalDistance.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/EllipsoidalDistance.hs
@@ -0,0 +1,4 @@
+module Data.Geo.Accessor.EllipsoidalDistance where
+
+class EllipsoidalDistance a where
+  ellipsoidalDistance:: a -> Double
diff --git a/Data/Geo/Accessor/Lat.hs b/Data/Geo/Accessor/Lat.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Lat.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Lat where
+
+import Data.Geo.Latitude
+
+class Lat a where
+  lat :: a -> Latitude
diff --git a/Data/Geo/Accessor/Lon.hs b/Data/Geo/Accessor/Lon.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/Lon.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.Lon where
+
+import Data.Geo.Longitude
+
+class Lon a where
+  lon :: a -> Longitude
diff --git a/Data/Geo/Accessor/ReverseAzi.hs b/Data/Geo/Accessor/ReverseAzi.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Accessor/ReverseAzi.hs
@@ -0,0 +1,6 @@
+module Data.Geo.Accessor.ReverseAzi where
+
+import Data.Geo.Azimuth
+
+class ReverseAzi a where
+  reverseAzi :: a -> Azimuth
diff --git a/Data/Geo/Azimuth.hs b/Data/Geo/Azimuth.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Azimuth.hs
@@ -0,0 +1,20 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+
+-- | An azimuth in degrees between 0 and 360).
+module Data.Geo.Azimuth(
+                         Azimuth,
+                         azimuth
+                       ) where
+
+import Data.Fixed
+import Data.Geo.Accessor.Value
+
+newtype Azimuth = Azimuth Double
+  deriving (Eq, Ord, Enum, Show, Num, Fractional, Floating)
+
+-- | Construct an azimuth with the number of degrees.
+azimuth :: Double -> Azimuth
+azimuth x = Azimuth (x `mod'` 360)
+
+instance Value Azimuth Double where
+  value (Azimuth x) = x
diff --git a/Data/Geo/Bearing.hs b/Data/Geo/Bearing.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Bearing.hs
@@ -0,0 +1,20 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+
+-- | A bearing in a degrees between 0 and 360.
+module Data.Geo.Bearing(
+                         Bearing,
+                         bearing
+                       ) where
+
+import Data.Geo.Accessor.Value
+import Data.Fixed
+
+newtype Bearing = Bearing Double
+  deriving (Eq, Ord, Enum, Show, Num, Fractional, Floating)
+
+-- | Construct a bearing with the number of degrees.
+bearing :: Double -> Bearing
+bearing x = Bearing (x `mod'` 360)
+
+instance Value Bearing Double where
+  value (Bearing x) = x
diff --git a/Data/Geo/Coord.hs b/Data/Geo/Coord.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Coord.hs
@@ -0,0 +1,51 @@
+-- | A coordinate on earth with a latitude and longitude.
+module Data.Geo.Coord(
+                       Coord,
+                       (|.|),
+                       (!.!),
+                       (|..|),
+                       radians'
+                     ) where
+
+import Data.Geo.Latitude
+import Data.Geo.Longitude
+import Data.Geo.Radians
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Data.Geo.Accessor.Value
+import Control.Applicative
+import Control.Arrow
+
+data Coord = Coord Latitude Longitude
+  deriving Eq
+
+instance Lat Coord where
+  lat (Coord x _) = x
+
+instance Lon Coord where
+  lon (Coord _ x) = x
+
+instance Show Coord where
+  show = show . (value . lat &&& value . lon)
+
+-- | Construct a coordinate with the given latitude and longitude.
+(|.|) :: Latitude -> Longitude -> Coord
+(|.|) = Coord
+
+infixl 8 |.|
+
+-- | -- | Construct a coordinate with the given latitude and longitude in degrees.
+(!.!) :: Double -> Double -> Coord
+(!.!) x y = latitude x |.| longitude y
+
+infixl 8 !.!
+
+-- | -- | Construct a coordinate with the given latitude and longitude in radians.
+(|..|) :: Double -> Double -> Coord
+(|..|) lat' lon' = fromRadians lat' |.| fromRadians lon'
+
+infixl 8 |..|
+
+-- | Convert the latitude and longitude of the given coordinate to radians.
+radians' :: Coord -> (Double, Double)
+radians' = radians <$> lat <*> lon
diff --git a/Data/Geo/DMS.hs b/Data/Geo/DMS.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/DMS.hs
@@ -0,0 +1,80 @@
+-- | Represent a value in degrees, minutes and seconds.
+module Data.Geo.DMS(
+                     DMS,
+                     positive,
+                     degrees,
+                     minutes,
+                     seconds,
+                     dms,
+                     DMSable,
+                     toDMS,
+                     fromDMS,
+                     showNegPos,
+                     showDMS,
+                     coordDMS
+                   ) where
+
+import Data.Geo.Latitude
+import Data.Geo.Longitude
+import Data.Geo.Coord
+import Data.Geo.Accessor.Value
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Control.Arrow
+import Text.Printf
+import Data.List
+
+-- | The structure of a type convertible to degrees, minutes and seconds.
+data DMS = DMS {
+  positive :: Bool,
+  degrees :: Int,
+  minutes :: Int,
+  seconds :: Double
+} deriving (Eq, Ord, Show)
+
+-- | Construct a value of degrees, minutes and seconds.
+dms :: Bool -> Int -> Int -> Double -> DMS
+dms p d m s = DMS p (abs d) (abs m) (abs s)
+
+class DMSable a where
+  toDMS :: a -> DMS
+  fromDMS :: DMS -> a
+  showNegPos :: a -> String
+
+instance DMSable Latitude where
+  toDMS = toDMS' . value
+  fromDMS = latitude . fromDMS'
+  showNegPos x = if value x < 0 then "S" else "N"
+
+instance DMSable Longitude where
+  toDMS = toDMS' . value
+  fromDMS = longitude . fromDMS'
+  showNegPos x = if value x < 0 then "W" else "E"
+
+-- | Show a value in degrees, minutes and seconds.
+showDMS :: (DMSable a) => a -> String
+showDMS x = let DMS _ d m s = toDMS x
+                d' = show d
+                m' = show m
+                s' = printf "%.5f" s
+                suffix = showNegPos x
+            in intercalate "," [d', m', s'] ++ suffix
+
+-- | Show a coordinate in degrees, minutes and seconds.
+coordDMS :: Coord -> String
+coordDMS = show . (showDMS *** showDMS) . (lat &&& lon)
+
+-- not exported
+fromDMS' :: DMS -> Double
+fromDMS' d = let s = seconds d / 3600
+                 m = (/ 60) . fromIntegral . minutes $ d
+                 d' = fromIntegral . degrees $ d
+                 z = if positive d then id else negate
+             in z $ s + m + d'
+
+-- not exported
+toDMS' :: Double -> DMS
+toDMS' x = let (n, f) = second abs . properFraction $ x
+               (m, g) = properFraction . (* 60) $ f
+           in dms (n >= 0) n m (g * 60)
+
diff --git a/Data/Geo/ElevatedCurve.hs b/Data/Geo/ElevatedCurve.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/ElevatedCurve.hs
@@ -0,0 +1,32 @@
+-- | An elevated curve is a geodetic curve with an elevation in metres.
+module Data.Geo.ElevatedCurve(
+                               ElevatedCurve,
+                               elevatedCurve,
+                               curveLength
+                             ) where
+
+import Data.Geo.GeodeticCurve
+import Data.Geo.Elevation
+import Data.Geo.Accessor.Curve
+import Data.Geo.Accessor.Ele
+import Data.Geo.Accessor.EllipsoidalDistance
+import Data.Geo.Accessor.Value
+
+data ElevatedCurve = ElevatedCurve GeodeticCurve Elevation
+  deriving (Eq, Show)
+
+instance Curve ElevatedCurve where
+  curve (ElevatedCurve x _) = x
+
+instance Ele ElevatedCurve where
+  ele (ElevatedCurve _ x) = x
+
+-- | Construct an elevated curve.
+elevatedCurve :: GeodeticCurve -> Elevation -> ElevatedCurve
+elevatedCurve = ElevatedCurve
+
+-- | Compute the length of an elevated curve.
+curveLength :: ElevatedCurve -> Double
+curveLength c = let d = ellipsoidalDistance . curve $ c
+                    square x = x * x
+                in sqrt (square d + square (value . ele $ c))
diff --git a/Data/Geo/Elevation.hs b/Data/Geo/Elevation.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Elevation.hs
@@ -0,0 +1,19 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+
+-- | An elevation in metres, above 0.
+module Data.Geo.Elevation(
+                           Elevation,
+                           elevation
+                         ) where
+
+import Data.Geo.Accessor.Value
+
+newtype Elevation = Elevation Double
+  deriving (Eq, Ord, Enum, Show, Num, Fractional, Floating)
+
+-- | Construct an elevation with the given number of metres.
+elevation :: Double -> Elevation
+elevation = Elevation . abs
+
+instance Value Elevation Double where
+  value (Elevation x) = x
diff --git a/Data/Geo/Ellipsoid.hs b/Data/Geo/Ellipsoid.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Ellipsoid.hs
@@ -0,0 +1,100 @@
+-- | An ellipsoid with parameters such as flattening.
+module Data.Geo.Ellipsoid(
+                           Ellipsoid,
+                           semiMajor,
+                           semiMinor,
+                           flattening,
+                           inverseFlattening,
+                           ellipsoid,
+                           semiMajorInverseF,
+                           semiMajorFlattening,
+                           wgs84,
+                           grs80,
+                           grs67,
+                           ans,
+                           wgs72,
+                           au1965,
+                           krasovsky1940,
+                           international1924,
+                           hayford1909,
+                           airy1830,
+                           everest1830,
+                           bessel1841,
+                           clarke1858,
+                           clarke1866,
+                           clarke1880
+                         ) where
+
+-- | Ellipsoidal parameters. Some are derivable from others.
+data Ellipsoid = Ellipsoid {
+  semiMajor :: Double, -- ^ The semi major axis in metres.
+  semiMinor :: Double, -- ^ The semi minor axis in metres.
+  flattening :: Double, -- ^ The ellipsoidal flattening.
+  inverseFlattening :: Double -- ^ The inverse of the ellipsoidal flattening.
+} deriving (Eq, Show)
+
+-- | Construct an ellipsoid with the given parameters.
+ellipsoid :: Double -- ^ The semi major axis in metres.
+             -> Double -- ^ The semi minor axis in metres.
+             -> Double -- ^ The ellipsoidal flattening.
+             -> Double -- ^ The inverse of the ellipsoidal flattening.
+             -> Ellipsoid
+ellipsoid = Ellipsoid
+
+-- | Construct an ellipsoid using only a semi major axis and inverse flattening. Other parameters are computed.
+semiMajorInverseF :: Double -- ^ The semi major axis.
+                     -> Double -- ^ The inverse of the ellipsoidal flattening.
+                     -> Ellipsoid
+semiMajorInverseF s i = let f = 1 / i
+                        in ellipsoid s ((1.0 - f) * s) f i
+
+-- | Construct an ellipsoid using only a semi major axis and flattening. Other parameters are computed.
+semiMajorFlattening :: Double -- ^ The semi major axis.
+                       -> Double -- ^ The ellipsoidal flattening.
+                       -> Ellipsoid
+semiMajorFlattening s f = ellipsoid s ((1.0 - f) * s) f (1 / f)
+
+wgs84 :: Ellipsoid
+wgs84 = semiMajorInverseF 6378137 298.257223563
+
+grs80 :: Ellipsoid
+grs80 = semiMajorInverseF 6378137 298.257222101
+
+grs67 :: Ellipsoid
+grs67 = semiMajorInverseF 6378160 298.25
+
+ans :: Ellipsoid
+ans = semiMajorInverseF 6378160 298.25
+
+wgs72 :: Ellipsoid
+wgs72 = semiMajorInverseF 6378135 298.26
+
+au1965 :: Ellipsoid
+au1965 = semiMajorInverseF 6378160 298.25
+
+krasovsky1940 :: Ellipsoid
+krasovsky1940 = semiMajorInverseF 6378245 298.3
+
+international1924 :: Ellipsoid
+international1924 = semiMajorInverseF 6378388 297
+
+hayford1909 :: Ellipsoid
+hayford1909 = international1924
+
+airy1830 :: Ellipsoid
+airy1830 = semiMajorInverseF 6377563.4 299.32
+
+everest1830 :: Ellipsoid
+everest1830 = semiMajorInverseF 6377276.3 300.8
+
+bessel1841 :: Ellipsoid
+bessel1841 = semiMajorInverseF 6377397.2 299.15
+
+clarke1858 :: Ellipsoid
+clarke1858 = semiMajorInverseF 6378293.645 294.26
+
+clarke1866 :: Ellipsoid
+clarke1866 = semiMajorInverseF 6378206.4 294.98
+
+clarke1880 :: Ellipsoid
+clarke1880 = semiMajorInverseF 6378249.145 293.465
diff --git a/Data/Geo/GeodeticCurve.hs b/Data/Geo/GeodeticCurve.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/GeodeticCurve.hs
@@ -0,0 +1,29 @@
+-- | A geodetic curve is made of a distance in metres, an azimuth and a reverse azimuth.
+module Data.Geo.GeodeticCurve(
+                               GeodeticCurve,
+                               geodeticCurve
+                             ) where
+
+import Data.Geo.Azimuth
+import Data.Geo.Accessor.EllipsoidalDistance
+import Data.Geo.Accessor.Azi
+import Data.Geo.Accessor.ReverseAzi
+
+data GeodeticCurve = GeodeticCurve Double Azimuth Azimuth
+  deriving (Eq, Show)
+
+-- | Construct a geodetic curve with the given parameters.
+geodeticCurve :: Double -- ^ The ellipsoidal distance.
+                 -> Azimuth -- ^ The azimuth.
+                 -> Azimuth -- ^ The reverse azimuth.
+                 -> GeodeticCurve
+geodeticCurve = GeodeticCurve . abs
+
+instance EllipsoidalDistance GeodeticCurve where
+  ellipsoidalDistance (GeodeticCurve x _ _) = x
+
+instance Azi GeodeticCurve where
+  azi (GeodeticCurve _ x _) = x
+
+instance ReverseAzi GeodeticCurve where
+  reverseAzi (GeodeticCurve _ _ x) = x
diff --git a/Data/Geo/GreatCircle.hs b/Data/Geo/GreatCircle.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/GreatCircle.hs
@@ -0,0 +1,17 @@
+module Data.Geo.GreatCircle(
+                             sphericalLaw
+                           ) where
+
+import Data.Geo.Sphere
+import Data.Geo.Coord
+import Data.Geo.Radians
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Data.Geo.Accessor.Value
+
+sphericalLaw :: Sphere -> Coord -> Coord -> Double
+sphericalLaw s start end = let lat1 = toRadians (lat start)
+                               lat2 = toRadians (lat end)
+                               lon1 = toRadians (lon start)
+                               lon2 = toRadians (lon end)
+                           in acos (sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)) * value s
diff --git a/Data/Geo/Haversine.hs b/Data/Geo/Haversine.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Haversine.hs
@@ -0,0 +1,21 @@
+module Data.Geo.Haversine(
+                           haversine
+                         ) where
+
+import Data.Geo.Sphere
+import Data.Geo.Coord
+import Data.Geo.Radians
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Data.Geo.Accessor.Value
+
+haversine :: Sphere -> Coord -> Coord -> Double
+haversine s start end = let lat1 = lat start
+                            lat2 = lat end
+                            dlat = (toRadians (lat1 - lat2)) / 2
+                            dlon = (toRadians (lon start - lon end)) / 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 value s * c
diff --git a/Data/Geo/Latitude.hs b/Data/Geo/Latitude.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Latitude.hs
@@ -0,0 +1,20 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+
+-- | A latitude on earth in degrees.
+module Data.Geo.Latitude(
+                          Latitude,
+                          latitude
+                        ) where
+
+import Data.Geo.Accessor.Value
+import Data.Fixed
+
+newtype Latitude = Latitude Double
+  deriving (Eq, Ord, Enum, Show, Num, Fractional, Floating)
+
+-- | Construct a latitude using the given number of degrees between -90 and 90.
+latitude :: Double -> Latitude
+latitude x = Latitude ((x + 90) `mod'` 180 - 90)
+
+instance Value Latitude Double where
+  value (Latitude x) = x
diff --git a/Data/Geo/Longitude.hs b/Data/Geo/Longitude.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Longitude.hs
@@ -0,0 +1,20 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+
+-- | A longitude on earth in degrees.
+module Data.Geo.Longitude(
+                           Longitude,
+                           longitude
+                         ) where
+
+import Data.Geo.Accessor.Value
+import Data.Fixed
+
+newtype Longitude = Longitude Double
+  deriving (Eq, Ord, Enum, Show, Num, Fractional, Floating)
+
+-- | Construct a longitude using the given number of degrees between -180 and 180.
+longitude :: Double -> Longitude
+longitude x = Longitude ((x + 180) `mod'` 360 - 180)
+
+instance Value Longitude Double where
+  value (Longitude x) = x
diff --git a/Data/Geo/Position.hs b/Data/Geo/Position.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Position.hs
@@ -0,0 +1,25 @@
+-- | A coordinate on earth with an elevation in metres.
+module Data.Geo.Position(
+                          Position,
+                          (|*|)
+                        ) where
+
+import Data.Geo.Coord
+import Data.Geo.Elevation
+import Data.Geo.Accessor.Coordinate
+import Data.Geo.Accessor.Ele
+
+data Position = Position Coord Elevation
+  deriving (Eq, Show)
+
+instance Coordinate Position where
+  coordinate (Position x _) = x
+
+instance Ele Position where
+  ele (Position _ x) = x
+
+-- | Construct a position with the given coordinate and elevation in metres.
+(|*|) :: Coord -- ^ The coordinate.
+         -> Elevation -- ^ The elevation.
+         -> Position
+(|*|) = Position
diff --git a/Data/Geo/Radians.hs b/Data/Geo/Radians.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Radians.hs
@@ -0,0 +1,37 @@
+-- | Convert types to and from radians/degrees.
+module Data.Geo.Radians(
+                         Radians,
+                         fromRadians,
+                         toRadians,
+                         radians
+                       ) where
+
+import Data.Geo.Latitude
+import Data.Geo.Longitude
+import Data.Geo.Bearing
+import Data.Geo.Accessor.Value
+import Control.Arrow
+
+class Radians a where
+  toRadians :: a -> Double
+  fromRadians :: Double -> a
+
+instance Radians Double where
+  toRadians n = n * pi / 180
+  fromRadians n = n * 180 / pi
+
+instance Radians Latitude where
+  toRadians = toRadians . value
+  fromRadians = latitude . fromRadians
+
+instance Radians Longitude where
+  toRadians = toRadians . value
+  fromRadians = longitude . fromRadians
+
+instance Radians Bearing where
+  toRadians = toRadians . value
+  fromRadians = bearing . fromRadians
+
+-- | Convert two values to radians.
+radians :: (Radians a, Radians b) => a -> b -> (Double, Double)
+radians = curry (toRadians *** toRadians)
diff --git a/Data/Geo/Vincenty.hs b/Data/Geo/Vincenty.hs
new file mode 100644
--- /dev/null
+++ b/Data/Geo/Vincenty.hs
@@ -0,0 +1,196 @@
+{-# LANGUAGE FlexibleInstances #-}
+
+-- | An implementation of Thaddeus Vincenty's direct and inverse geodetic algorithms. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>
+module Data.Geo.Vincenty(
+                          convergence,
+                          VincentyDirect,
+                          direct,
+                          VincentyInverse,
+                          inverse
+                        ) where
+
+import Data.Geo.Ellipsoid
+import Data.Geo.Coord
+import Data.Geo.Bearing
+import Data.Geo.Radians
+import Data.Geo.GeodeticCurve
+import Data.Geo.Accessor.Lat
+import Data.Geo.Accessor.Lon
+import Control.Arrow
+import Control.Monad
+
+-- | An acceptable convergence value.
+convergence :: Double
+convergence = 0.0000000000001
+
+class VincentyDirect a where
+  direct :: a -> Coord -> Bearing -> Double -> (Coord, Bearing)
+
+instance VincentyDirect (Double, Ellipsoid) where
+  direct (conv, e) start bear dist =
+    let sMnr = semiMinor e
+        flat = flattening e
+        alpha = toRadians bear
+        cosAlpha = cos alpha
+        sinAlpha = sin alpha
+        tanu1 = (1.0 - flat) * tan (toRadians . lat $ start)
+        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 (semiMajor e / 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' = fromRadians (atan2 (sinu1 * cosSigma + cosu1 * sinSigma * cosAlpha) ((1.0 - flat) * sqrt (sin2Alpha + (sss - ccca) ** 2.0)))
+        longitude = lon start + fromRadians (atan2 (sinSigma * sinAlpha) (cc - sss * cosAlpha) - (1 - c) * flat * csa * (sigma'' + c * sinSigma * (cosSigmaM2 + c * cosSigma * (-1 + 2 * cos2SigmaM2))))
+    in (latitude' |.| longitude, fromRadians (atan2 csa (ccca - sss)))
+
+instance VincentyDirect Double where
+  direct d = direct (d, wgs84)
+
+instance VincentyDirect Ellipsoid where
+  direct e = direct (convergence, e)
+
+instance VincentyDirect () where
+  direct _ = direct (convergence, wgs84)
+
+class VincentyInverse a where
+  inverse :: a -> Coord -> Coord -> GeodeticCurve
+
+instance VincentyInverse (Double, Ellipsoid) where
+  inverse (conv, e) start end =
+    let b = semiMinor e
+        f = flattening e
+        (phi1, phi2) = let rl k = toRadians . lat $ k
+                      in (rl start, rl end)
+        a2b2b2 = let ss z = square (z e)
+                in ss semiMajor / ss semiMinor - 1
+        omega = let rl k = toRadians . lon $ k
+                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 c cp x y = join (***) (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 geodeticCurve (b * a' ed * (sigma ed - deltasigma ed)) alpha1 alpha2
+
+instance VincentyInverse Double where
+  inverse d = inverse (d, wgs84)
+
+instance VincentyInverse Ellipsoid where
+  inverse e = inverse (convergence, e)
+
+instance VincentyInverse () where
+  inverse _ = inverse (convergence, wgs84)
+
+data P = P {
+  origSigma' :: Double,
+  sigma' :: Double,
+  prevSigma' :: Double
+} deriving Show
+
+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 = join (*)
+
+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/Geodetic.cabal b/Geodetic.cabal
new file mode 100644
--- /dev/null
+++ b/Geodetic.cabal
@@ -0,0 +1,48 @@
+Name:                Geodetic
+Version:             0.1
+License:             BSD3
+License-File:        LICENSE
+Synopsis:            Geodetic calculations
+Description:         Geodetic calculations including Vincenty and Great Circle
+Homepage:            http://code.google.com/p/geodetic/
+Category:            Utils
+Author:              Tony Morris
+Maintainer:          code@tmorris.net
+Copyright:           2009 Tony Morris
+build-type:          Simple
+cabal-version:       >= 1.2
+
+Flag small_base
+  Description:     Choose the new, split-up base package.
+
+Library
+  if flag(small_base)
+    Build-Depends: base < 4 && >= 3
+  else
+    Build-Depends: base < 3
+
+  GHC-Options:    -Wall
+  Exposed-Modules: Data.Geo
+                   Data.Geo.Azimuth
+                   Data.Geo.Bearing
+                   Data.Geo.Coord
+                   Data.Geo.DMS
+                   Data.Geo.ElevatedCurve
+                   Data.Geo.Elevation
+                   Data.Geo.Ellipsoid
+                   Data.Geo.GeodeticCurve
+                   Data.Geo.GreatCircle
+                   Data.Geo.Haversine
+                   Data.Geo.Latitude
+                   Data.Geo.Longitude
+                   Data.Geo.Position
+                   Data.Geo.Radians
+                   Data.Geo.Vincenty
+                   Data.Geo.Accessor.Azi
+                   Data.Geo.Accessor.Coordinate
+                   Data.Geo.Accessor.Curve
+                   Data.Geo.Accessor.Ele
+                   Data.Geo.Accessor.EllipsoidalDistance
+                   Data.Geo.Accessor.Lat
+                   Data.Geo.Accessor.Lon
+                   Data.Geo.Accessor.ReverseAzi
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright 2009 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/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain