packages feed

geodetic 0.1.2 → 0.1.3

raw patch · 10 files changed

+301/−226 lines, 10 filesdep +template-haskelldep ~coordinatedep ~lens

Dependencies added: template-haskell

Dependency ranges changed: coordinate, lens

Files

+ changelog view
@@ -0,0 +1,9 @@+0.1.3++  * Update lens dependency to >= 4.0++  * Update coordinate dependency to >= 0.0.16+  +0.1.2++* Release under NICTA open source.
geodetic.cabal view
@@ -1,5 +1,5 @@ name:               geodetic-version:            0.1.2+version:            0.1.3 license:            BSD3 license-File:       etc/LICENCE author:             Tony Morris <ʇǝu˙sıɹɹoɯʇ@ןןǝʞsɐɥ>@@ -15,6 +15,7 @@ bug-reports:        https://github.com/NICTA/geodetic/issues cabal-version:      >= 1.10 build-type:         Custom+extra-source-files: changelog  source-repository   head   type:             git@@ -29,8 +30,8 @@    build-depends:                     base < 5 && >= 3-                    , lens >= 3.10-                    , coordinate >= 0.0.8+                    , lens >= 4.0+                    , coordinate >= 0.0.16                     , optional >= 0.0.1    ghc-options:@@ -68,7 +69,8 @@                     doctest >= 0.9.7,                     filepath >= 1.3,                     directory >= 1.1,-                    QuickCheck >= 2.0+                    QuickCheck >= 2.0,+                    template-haskell >= 2.8    ghc-options:                     -Wall
src/Data/Geo/Geodetic/Azimuth.hs view
@@ -1,14 +1,18 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}++ -- | An azimuth in degrees between 0 and 360. module Data.Geo.Geodetic.Azimuth(   Azimuth-, HasAzimuth(..)+, AsAzimuth(..) , modAzimuth-, nAzimuth ) where +import Control.Applicative(Applicative) import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), showParen, showString) import Data.Maybe(Maybe(..))-import Control.Lens(Prism', Lens', prism')+import Control.Lens(Choice, Optic', prism') import Text.Printf(printf) import Data.Fixed(mod') @@ -50,34 +54,33 @@ modAzimuth x =   Azimuth (x `mod'` 360) +class AsAzimuth p f s where+  _Azimuth ::+    Optic' p f s Azimuth++instance AsAzimuth p f Azimuth where+  _Azimuth =+    id+ -- | A prism on azimuth to an integer between 0 and 359 inclusive. ----- >>> 7 ^? nAzimuth+-- >> 7 ^? _Azimuth -- Just (Azimuth 7.0000) ----- >>> 0 ^? nAzimuth+-- >> 0 ^? _Azimuth -- Just (Azimuth 0.0000) ----- >>> 359.999 ^? nAzimuth+-- >> 359.999 ^? _Azimuth -- Just (Azimuth 359.9990) ----- >>> 360 ^? nAzimuth+-- >> 360 ^? _Azimuth -- 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+-- prlop> all (\m -> _Azimuth # m == n) (n ^? _Azimuth)+instance (Choice p, Applicative f) => AsAzimuth p f Double where+  _Azimuth =+    prism'+      (\(Azimuth i) -> i)+      (\i -> case i >= 0 && i < 360 of+               True -> Just (Azimuth i)+               False -> Nothing)
src/Data/Geo/Geodetic/Bearing.hs view
@@ -1,15 +1,18 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+ -- | A bearing in degrees between 0 and 360. module Data.Geo.Geodetic.Bearing(   Bearing-, HasBearing(..)+, AsBearing(..) , modBearing-, degreeBearing , radianBearing ) where +import Control.Applicative(Applicative) 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 Control.Lens(Choice, Optic', Prism', prism', iso) import Text.Printf(printf) import Data.Fixed(mod') @@ -54,33 +57,6 @@ 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@@ -106,12 +82,39 @@ radianBearing ::   Prism' Double Bearing radianBearing =-  iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . degreeBearing+  iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Bearing -class HasBearing t where-  bearing ::-    Lens' t Bearing+class AsBearing p f s where+  _Bearing ::+    Optic' p f s Bearing -instance HasBearing Bearing where-  bearing =+instance AsBearing p f Bearing where+  _Bearing =     id++-- | A prism on bearing to a double between 0 inclusive and 360 exclusive.+--+-- >>> (7 :: Double) ^? _Bearing+-- Just (Bearing 7.0000)+--+-- >>> (0 :: Double) ^? _Bearing+-- Just (Bearing 0.0000)+--+-- >>> (359 :: Double) ^? _Bearing+-- Just (Bearing 359.0000)+--+-- >>> (359.997 :: Double) ^? _Bearing+-- Just (Bearing 359.9970)+--+-- >>> (360 :: Double) ^? _Bearing+-- Nothing+--+-- prop> all (\m -> _Bearing # m == n) ((n :: Double) ^? _Bearing)+instance (Choice p, Applicative f) => AsBearing p f Double where+  _Bearing =+    prism'+      (\(Bearing i) -> i)+      (\i -> case i >= 0 && i < 360 of+               True -> Just (Bearing i)+               False -> Nothing)+
src/Data/Geo/Geodetic/Curve.hs view
@@ -1,16 +1,22 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+ -- | A geodetic curve is made of a distance in metres, an azimuth and a reverse azimuth. module Data.Geo.Geodetic.Curve(   Curve+, AsCurve(..) , curve , curveDistance , curveAzimuth , curveReverseAzimuth ) where -import Prelude(Eq, Show(..), Ord(..), Double, (.), abs, showString, showParen)+import Prelude(Eq, Show(..), Ord(..), Double, showString, showParen, id)+import Data.Functor(Functor) import Data.List(unwords) import Text.Printf(printf)-import Control.Lens(Lens', lens)+import Control.Lens(Optic', Lens', Profunctor, lens, iso) import Data.Geo.Geodetic.Azimuth  data Curve =@@ -49,3 +55,17 @@   Lens' Curve Azimuth curveReverseAzimuth =   lens (\(Curve _ _ r) -> r) (\(Curve d a _) r -> Curve d a r)++class AsCurve p f s where+  _Curve ::+    Optic' p f s Curve++instance AsCurve p f Curve where+  _Curve =+    id++instance (Profunctor p, Functor f) => AsCurve p f (Double, Azimuth, Azimuth) where+  _Curve =+    iso+      (\(d, a, r) -> Curve d a r)+      (\(Curve d a r) -> (d, a, r))
src/Data/Geo/Geodetic/Ellipsoid.hs view
@@ -1,12 +1,16 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+ module Data.Geo.Geodetic.Ellipsoid(   -- * Data type   Ellipsoid-, HasEllipsoid(..)+, AsEllipsoid(..)   -- * Ellipsoid properties-, HasSemiMajor(..)-, HasSemiMinor(..)-, HasFlattening(..)-, HasInverseFlattening(..)+, AsSemiMajor(..)+, AsSemiMinor(..)+, AsFlattening(..)+, AsInverseFlattening(..)   -- * Ellipsoid construction , semiMajorFlattening , semiMinorFlattening@@ -30,8 +34,8 @@ , clarke1880 ) where -import Prelude(Eq, Ord, Show, Num(..), Fractional(..), Double, id)-import Control.Lens(Lens', lens)+import Prelude(Eq, Ord, Show, Num(..), Fractional(..), Double, id, Functor)+import Control.Lens(Optic', lens)  data Ellipsoid =   Ellipsoid@@ -67,61 +71,61 @@ semiMinorInverseFlattening s f =   semiMinorFlattening s (1/f) -class HasEllipsoid t where-  ellipsoid ::-    Lens' t Ellipsoid+class AsEllipsoid p f s where+  _Ellipsoid ::+    Optic' p f s Ellipsoid     -instance HasEllipsoid Ellipsoid where-  ellipsoid =+instance AsEllipsoid p f Ellipsoid where+  _Ellipsoid =     id -class HasSemiMajor t where-  semiMajor ::-    Lens' t Double+class AsSemiMajor p f s where+  _SemiMajor ::+    Optic' p f s Double -instance HasSemiMajor Double where-  semiMajor =+instance AsSemiMajor p f Double where+  _SemiMajor =     id -instance HasSemiMajor Ellipsoid where-  semiMajor =-    lens (\(Ellipsoid s _) -> s) (\(Ellipsoid _ f) s -> Ellipsoid s f)+instance (p ~ (->), Functor f) => AsSemiMajor p f Ellipsoid where+  _SemiMajor =+     lens (\(Ellipsoid s _) -> s) (\(Ellipsoid _ f) s -> Ellipsoid s f)    -class HasSemiMinor t where-  semiMinor ::-    Lens' t Double+class AsSemiMinor p f s where+  _SemiMinor ::+    Optic' p f s Double -instance HasSemiMinor Double where-  semiMinor =-    id+instance AsSemiMinor p f 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)+instance (p ~ (->), Functor f) => AsSemiMinor p f 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+class AsFlattening p f s where+  _Flattening ::+    Optic' p f s Double -instance HasFlattening Double where-  flattening =+instance AsFlattening p f Double where+  _Flattening =     id -instance HasFlattening Ellipsoid where-  flattening =+instance (p ~ (->), Functor f) => AsFlattening p f Ellipsoid where+  _Flattening =     lens (\(Ellipsoid _ f) -> f) (\(Ellipsoid s _) f -> Ellipsoid s f) -class HasInverseFlattening t where-  inverseFlattening ::-    Lens' t Double+class AsInverseFlattening p f s where+  _InverseFlattening ::+    Optic' p f s Double -instance HasInverseFlattening Double where-  inverseFlattening =+instance AsInverseFlattening p f Double where+  _InverseFlattening =     id -instance HasInverseFlattening Ellipsoid where-  inverseFlattening =-    lens (\(Ellipsoid _ f) -> 1/f) (\(Ellipsoid s _) f -> Ellipsoid s (1/f))+instance (p ~ (->), Functor f) => AsInverseFlattening p f Ellipsoid where+  _InverseFlattening =+    lens (\(Ellipsoid _ f) -> f) (\(Ellipsoid s _) f -> Ellipsoid s f)  wgs84 ::   Ellipsoid
src/Data/Geo/Geodetic/GreatCircle.hs view
@@ -8,6 +8,7 @@ ) where  import Prelude(Double, Num(..), Fractional(..), pi, sin, cos, acos)+import Control.Applicative(Const) import Control.Lens((#), (^.)) import System.Args.Optional(Optional1(..)) import Data.Geo.Coordinate@@ -20,54 +21,54 @@  -- | 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+-- >>> 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+-- >>> 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+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (sphericalLaw ((6350000 :: Double) ^. _Sphere) 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+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (sphericalLaw ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String -- Just "17081801.7377" sphericalLaw ::-  (HasCoordinate c1, HasCoordinate c2) =>+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>   Sphere -- ^ reference sphere-  -> c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  -> start -- ^ start coordinate+  -> end -- ^ end coordinate   -> Double sphericalLaw s start' end' =-  let start = start' ^. coordinate-      end = end' ^. coordinate+  let start = start' ^. _Coordinate+      end = end' ^. _Coordinate       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+      lat1 = toRadians (_Latitude # (start ^. _Latitude))+      lat2 = toRadians (_Latitude # (end ^. _Latitude))+      lon1 = toRadians (_Longitude # (start ^. _Longitude))+      lon2 = toRadians (_Longitude # (end ^. _Longitude))+  in acos (sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1)) * _Sphere # 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+-- >>> 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+-- >>> 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 ::-  (HasCoordinate c1, HasCoordinate c2) =>-  c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>+  start -- ^ start coordinate+  -> end -- ^ end 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+-- >>> 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+-- >>> 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
src/Data/Geo/Geodetic/Haversine.hs view
@@ -1,4 +1,6 @@ {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}  -- | Haversine geodetic distance algorithm. module Data.Geo.Geodetic.Haversine(@@ -7,71 +9,72 @@ , haversine' ) where -import Prelude(Double, Num(..), Fractional(..), (.), pi, atan, sin, atan2, cos, sqrt)+import Control.Applicative(Const)+import Prelude(Double, Num(..), Fractional(..), (.), pi, 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 Prelude(Functor(..), Monad(..), String, Double) -- >>> 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+-- >>> 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+-- >>> 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+-- >>> fmap (printf "%0.4f") (do fr <- 27.812 <°> 154.295; to <- (-66.093) <°> 12.84; return (haversine ((6350000 :: Double) ^. _Sphere) 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+-- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) <°> 41.935; to <- 6.933 <°> (-162.55); return (haversine ((6350000 :: Double) ^. _Sphere) fr to)) :: Maybe String -- Just "17081801.7377" haversine ::-  (HasCoordinate c1, HasCoordinate c2) =>+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>   Sphere -- ^ reference sphere-  -> c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  -> start -- ^ start coordinate+  -> end -- ^ end coordinate   -> Double haversine s start' end' =-  let start = start' ^. coordinate-      end = end' ^. coordinate-      lat1 = fracLatitude # (start ^. latitude)-      lat2 = fracLatitude # (end ^. latitude)+  let start = start' ^. _Coordinate+      end = end' ^. _Coordinate+      lat1 = _Latitude # (start ^. _Latitude)+      lat2 = _Latitude # (end ^. _Latitude)       toRadians n = n * pi / 180       dlat = (toRadians (lat1 - lat2)) / 2-      dlon = (toRadians (fracLongitude # (start ^. longitude) - fracLongitude # (end ^. longitude))) / 2+      dlon = (toRadians (_Longitude # (start ^. _Longitude) - _Longitude # (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+  in (_Sphere # 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+-- >>> 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+-- >>> 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 ::-  (HasCoordinate c1, HasCoordinate c2) =>-  c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>+  start -- ^ start coordinate+  -> end -- ^ end 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+-- >>> 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+-- >>> 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
src/Data/Geo/Geodetic/Sphere.hs view
@@ -1,13 +1,17 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}++ -- | A sphere with a radius in metres. module Data.Geo.Geodetic.Sphere(   Sphere-, HasSphere(..)-, nSphere+, AsSphere(..) , earthMean ) where -import Prelude(Double, Bool(..), Eq, Show(..), Ord(..), id, (&&), (++), (.), abs, showParen, showString)-import Control.Lens(Iso', Lens', iso)+import Prelude(Double, Eq, Show(..), Ord(..), id, (++), showParen, showString)+import Control.Lens(Optic', Profunctor, iso)+import Data.Functor(Functor) import Text.Printf(printf)  -- $setup@@ -25,30 +29,31 @@   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+class AsSphere p f s where+  _Sphere ::+    Optic' p f s Sphere -instance HasSphere Sphere where-  sphere =+instance AsSphere p f Sphere where+  _Sphere =     id++-- | An isomorphism on sphere to a double.+--+-- >>> (7 :: Double) ^. _Sphere+-- Sphere 7.0000+--+-- >>> (0 :: Double) ^. _Sphere+-- Sphere 0.0000+--+-- >>> (-7 :: Double) ^. _Sphere+-- Sphere -7.0000+instance (Functor f, Profunctor p) => AsSphere p f Double where+  _Sphere =+    iso +      Sphere+      (\(Sphere d) -> d)
src/Data/Geo/Geodetic/Vincenty.hs view
@@ -1,4 +1,7 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}  -- | 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(@@ -13,8 +16,10 @@ , 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 Prelude(Eq(..), Show(..), Ord(..), Num(..), Floating(..), Fractional(..), Double, Int, Bool, Ordering(..), subtract, cos, sin, asin, tan, sqrt, atan, atan2, pi, (.), (++), (&&), ($!), error, id)+import Control.Applicative(Const)+import Control.Lens(Profunctor, Prism', Optic', Iso', (^.), (#), (^?), iso, _1, _2, from)+import Data.Functor(Functor) import Data.Maybe(fromMaybe) import System.Args.Optional(Optional2(..)) import Data.Geo.Coordinate@@ -41,29 +46,43 @@     Bearing   deriving (Eq, Ord, Show) -instance HasCoordinate VincentyDirectResult where-  coordinate =-    lens (\(VincentyDirectResult c _) -> c) (\(VincentyDirectResult _ b) c -> VincentyDirectResult c b)+class AsVincentyDirectResult p f s where+  _VincentyDirectResult ::+    Optic' p f s VincentyDirectResult -instance HasBearing VincentyDirectResult where-  bearing =-    lens (\(VincentyDirectResult _ b) -> b) (\(VincentyDirectResult c _) b -> VincentyDirectResult c b)+instance AsVincentyDirectResult p f VincentyDirectResult where+  _VincentyDirectResult =+    id +instance (Profunctor p, Functor f) => AsVincentyDirectResult p f (Coordinate, Bearing) where+  _VincentyDirectResult =+    iso+      (\(c, b) -> VincentyDirectResult c b)+      (\(VincentyDirectResult c b) -> (c, b))++instance (p ~ (->), Functor f) => AsCoordinate p f VincentyDirectResult where+  _Coordinate =+    from (_VincentyDirectResult :: Iso' (Coordinate, Bearing) VincentyDirectResult) . _1++instance (p ~ (->), Functor f) => AsBearing p f VincentyDirectResult where+  _Bearing =+    from (_VincentyDirectResult :: Iso' (Coordinate, Bearing) VincentyDirectResult) . _2+ -- | Vincenty direct algorithm. ----- >>> fmap (\c' -> direct wgs84 convergence c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)+-- >>> 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)+-- >>> 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)+-- >>> 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)+-- >>> 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 ::-  (HasEllipsoid e, HasCoordinate c, HasBearing b) =>+  (AsCoordinate (->) (Const Coordinate) c, AsBearing (->) (Const Bearing) b, AsEllipsoid (->) (Const Ellipsoid) e) =>   e -- ^ reference ellipsoid   -> Convergence -- ^ convergence point to stop calculating   -> c -- ^ begin coordinate@@ -71,22 +90,24 @@   -> Double -- ^ distance   -> VincentyDirectResult direct e' conv start' bear' dist =-  let e = e' ^. ellipsoid-      start = start' ^. coordinate-      bear = bear' ^. bearing-      sMnr = e ^. semiMinor-      flat = e ^. flattening+  let radianLatitude :: Prism' Double Latitude+      radianLatitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Latitude+      e = e' ^. _Ellipsoid+      start = start' ^. _Coordinate+      bear = bear' ^. _Bearing+      sMnr = e ^. _SemiMinor+      flat = e ^. _Flattening       alpha = radianBearing # bear       cosAlpha = cos alpha       sinAlpha = sin alpha-      tanu1 = (1.0 - flat) * tan (radianLatitude # (start ^. latitude))+      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)+      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@@ -111,8 +132,8 @@       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)+      longitude' = let r = _Longitude # (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 ^? _Longitude)   in VincentyDirectResult        (latitude' .#. longitude')        (@@ -122,13 +143,13 @@  -- | Vincenty direct algorithm with a default ellipsoid of WGS84 and standard convergence. ----- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) (27.812 ..#.. 154.295)+-- >>> 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)+-- >>> 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 ::-  (HasCoordinate c, HasBearing b) =>+  (AsCoordinate (->) (Const Coordinate) c, AsBearing (->) (Const Bearing) b) =>   c -- ^ begin coordinate   -> b -- ^ bearing   -> Double -- ^ distance@@ -138,10 +159,10 @@  -- | 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)+-- >>> 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)+-- >>> 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@@ -159,38 +180,42 @@  -- | Vincenty inverse algorithm. ----- >>> do fr <- 27.812 ..#.. 154.295; to <- (-66.093) ..#.. 12.84; return (inverse wgs84 convergence fr to)+-- >>> 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)+-- >>> 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)+-- >>> 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)+-- >>> 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 ::-  (HasEllipsoid e, HasCoordinate c1, HasCoordinate c2) =>+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end, AsEllipsoid (->) (Const Ellipsoid) e) =>   e -- ^ reference ellipsoid   -> Convergence -- ^ convergence point to stop calculating-  -> c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  -> start -- ^ start coordinate+  -> end -- ^ end coordinate   -> Curve inverse e' conv start' end' =-  let e = e' ^. ellipsoid-      start = start' ^. coordinate-      end = end' ^. coordinate-      b = e ^. semiMinor-      f = e ^. flattening+  let radianLatitude :: Prism' Double Latitude+      radianLatitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Latitude+      radianLongitude :: Prism' Double Longitude+      radianLongitude = iso (\n -> n * 180 / pi) (\n -> n * pi / 180) . _Longitude+      e = e' ^. _Ellipsoid+      start = start' ^. _Coordinate+      end = end' ^. _Coordinate+      b = e ^. _SemiMinor+      f = e ^. _Flattening       (phi1, phi2) =-        let rl k = radianLatitude # (k ^. latitude)+        let rl k = radianLatitude # (k ^. _Latitude)         in (rl start, rl end)       a2b2b2 =         let ss z = square (z e)-        in ss (^. semiMajor) / ss (^. semiMinor) - 1+        in ss (^. _SemiMajor) / ss (^. _SemiMinor) - 1       omega =-        let rl k = radianLongitude # (k ^. longitude)+        let rl k = radianLongitude # (k ^. _Longitude)         in rl end - rl start       (u1, u2) =         let at = atan . ((1 - f) *) . tan@@ -246,25 +271,25 @@  -- | 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)+-- >>> 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)+-- >>> 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 ::-  (HasCoordinate c1, HasCoordinate c2) =>-  c1 -- ^ start coordinate-  -> c2 -- ^ end coordinate+  (AsCoordinate (->) (Const Coordinate) start, AsCoordinate (->) (Const Coordinate) end) =>+  start -- ^ start coordinate+  -> end -- ^ end 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)+-- >>> 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)+-- >>> 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