packages feed

Geodetic (empty) → 0.1

raw patch · 27 files changed

+862/−0 lines, 27 filesdep +basesetup-changed

Dependencies added: base

Files

+ Data/Geo.hs view
@@ -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
+ Data/Geo/Accessor/Azi.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Azi where++import Data.Geo.Azimuth++class Azi a where+  azi :: a -> Azimuth
+ Data/Geo/Accessor/Coordinate.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Coordinate where++import Data.Geo.Coord++class Coordinate a where+  coordinate :: a -> Coord
+ Data/Geo/Accessor/Curve.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Curve where++import Data.Geo.GeodeticCurve++class Curve a where+  curve :: a -> GeodeticCurve
+ Data/Geo/Accessor/Ele.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Ele where++import Data.Geo.Elevation++class Ele a where+  ele :: a -> Elevation
+ Data/Geo/Accessor/EllipsoidalDistance.hs view
@@ -0,0 +1,4 @@+module Data.Geo.Accessor.EllipsoidalDistance where++class EllipsoidalDistance a where+  ellipsoidalDistance:: a -> Double
+ Data/Geo/Accessor/Lat.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Lat where++import Data.Geo.Latitude++class Lat a where+  lat :: a -> Latitude
+ Data/Geo/Accessor/Lon.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.Lon where++import Data.Geo.Longitude++class Lon a where+  lon :: a -> Longitude
+ Data/Geo/Accessor/ReverseAzi.hs view
@@ -0,0 +1,6 @@+module Data.Geo.Accessor.ReverseAzi where++import Data.Geo.Azimuth++class ReverseAzi a where+  reverseAzi :: a -> Azimuth
+ Data/Geo/Azimuth.hs view
@@ -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
+ Data/Geo/Bearing.hs view
@@ -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
+ Data/Geo/Coord.hs view
@@ -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
+ Data/Geo/DMS.hs view
@@ -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)+
+ Data/Geo/ElevatedCurve.hs view
@@ -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))
+ Data/Geo/Elevation.hs view
@@ -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
+ Data/Geo/Ellipsoid.hs view
@@ -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
+ Data/Geo/GeodeticCurve.hs view
@@ -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
+ Data/Geo/GreatCircle.hs view
@@ -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
+ Data/Geo/Haversine.hs view
@@ -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
+ Data/Geo/Latitude.hs view
@@ -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
+ Data/Geo/Longitude.hs view
@@ -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
+ Data/Geo/Position.hs view
@@ -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
+ Data/Geo/Radians.hs view
@@ -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)
+ Data/Geo/Vincenty.hs view
@@ -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+}
+ Geodetic.cabal view
@@ -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
+ LICENSE view
@@ -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.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain