geodetic 0.1.0 → 0.1.1
raw patch · 4 files changed
+64/−46 lines, 4 files
Files
- geodetic.cabal +1/−1
- src/Data/Geo/Geodetic/GreatCircle.hs +12/−8
- src/Data/Geo/Geodetic/Haversine.hs +12/−8
- src/Data/Geo/Geodetic/Vincenty.hs +39/−29
geodetic.cabal view
@@ -1,5 +1,5 @@ name: geodetic-version: 0.1.0+version: 0.1.1 license: BSD3 license-File: etc/LICENCE author: Tony Morris <ʇǝu˙sıɹɹoɯʇ@ןןǝʞsɐɥ>
src/Data/Geo/Geodetic/GreatCircle.hs view
@@ -32,12 +32,15 @@ -- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLaw (6350000 ^. nSphere) fr to)) :: Maybe String -- Just "17081801.7377" sphericalLaw ::- Sphere- -> Coordinate- -> Coordinate+ (HasCoordinate c1, HasCoordinate c2) =>+ Sphere -- ^ reference sphere+ -> c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Double-sphericalLaw s start end =- let toRadians n = n * pi / 180+sphericalLaw s start' end' =+ 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))@@ -52,8 +55,9 @@ -- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (sphericalLawD fr to)) :: Maybe String -- Just "17128743.0669" sphericalLawD ::- Coordinate- -> Coordinate+ (HasCoordinate c1, HasCoordinate c2) =>+ c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Double sphericalLawD = sphericalLaw earthMean@@ -75,4 +79,4 @@ ) x) => x sphericalLaw' =- optional1 sphericalLaw earthMean+ optional1 (sphericalLaw :: Sphere -> Coordinate -> Coordinate -> Double) earthMean
src/Data/Geo/Geodetic/Haversine.hs view
@@ -32,12 +32,15 @@ -- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversine (6350000 ^. nSphere) fr to)) :: Maybe String -- Just "17081801.7377" haversine ::- Sphere- -> Coordinate- -> Coordinate+ (HasCoordinate c1, HasCoordinate c2) =>+ Sphere -- ^ reference sphere+ -> c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Double-haversine s start end =- let lat1 = fracLatitude # (start ^. latitude)+haversine s start' end' =+ let start = start' ^. coordinate+ end = end' ^. coordinate+ lat1 = fracLatitude # (start ^. latitude) lat2 = fracLatitude # (end ^. latitude) toRadians n = n * pi / 180 dlat = (toRadians (lat1 - lat2)) / 2@@ -56,8 +59,9 @@ -- >>> fmap (printf "%0.4f") (do fr <- (-16.7889) ..#.. 41.935; to <- 6.933 ..#.. (-162.55); return (haversineD fr to)) :: Maybe String -- Just "17128743.0669" haversineD ::- Coordinate- -> Coordinate+ (HasCoordinate c1, HasCoordinate c2) =>+ c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Double haversineD = haversine earthMean@@ -79,4 +83,4 @@ ) x) => x haversine' =- optional1 haversine earthMean+ optional1 (haversine :: Sphere -> Coordinate -> Coordinate -> Double) earthMean
src/Data/Geo/Geodetic/Vincenty.hs view
@@ -63,14 +63,18 @@ -- >>> fmap (\c' -> direct ans convergence c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84) -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0662)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4139))) (Bearing 165.3183)) direct ::- Ellipsoid- -> Convergence- -> Coordinate- -> Bearing- -> Double+ (HasEllipsoid e, HasCoordinate c, HasBearing b) =>+ e -- ^ reference ellipsoid+ -> Convergence -- ^ convergence point to stop calculating+ -> c -- ^ begin coordinate+ -> b -- ^ bearing+ -> Double -- ^ distance -> VincentyDirectResult-direct e conv start bear dist =- let sMnr = e ^. semiMinor+direct e' conv start' bear' dist =+ let e = e' ^. ellipsoid+ start = start' ^. coordinate+ bear = bear' ^. bearing+ sMnr = e ^. semiMinor flat = e ^. flattening alpha = radianBearing # bear cosAlpha = cos alpha@@ -124,9 +128,10 @@ -- >>> fmap (\c' -> directD c' (modBearing 165.34) 4235) ((-66.093) ..#.. 12.84) -- Just (VincentyDirectResult (Coordinate (Latitude (DegreesLatitude (-66)) (Minutes 7) (Seconds 47.0667)) (Longitude (DegreesLongitude 12) (Minutes 51) (Seconds 49.4142))) (Bearing 165.3183)) directD ::- Coordinate- -> Bearing- -> Double+ (HasCoordinate c, HasBearing b) =>+ c -- ^ begin coordinate+ -> b -- ^ bearing+ -> Double -- ^ distance -> VincentyDirectResult directD = direct wgs84 convergence@@ -150,7 +155,7 @@ ) x) => x direct' =- optional2 direct wgs84 convergence+ optional2 (direct :: Ellipsoid -> Convergence -> Coordinate -> Bearing -> Double -> VincentyDirectResult) wgs84 convergence -- | Vincenty inverse algorithm. --@@ -166,13 +171,17 @@ -- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverse ans convergence fr to) -- Just (GeodeticCurve 7099229.9126 Azimuth 0.0000 Azimuth 180.0000) inverse ::- Ellipsoid- -> Convergence- -> Coordinate- -> Coordinate+ (HasEllipsoid e, HasCoordinate c1, HasCoordinate c2) =>+ e -- ^ reference ellipsoid+ -> Convergence -- ^ convergence point to stop calculating+ -> c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Curve-inverse e conv start end =- let b = e ^. semiMinor+inverse e' conv start' end' =+ let e = e' ^. ellipsoid+ start = start' ^. coordinate+ end = end' ^. coordinate+ b = e ^. semiMinor f = e ^. flattening (phi1, phi2) = let rl k = radianLatitude # (k ^. latitude)@@ -223,15 +232,15 @@ ed = whileDo iter pred' begin ifi p t a = if p a then t a else a (alpha1, alpha2) =- let alphaNoConverge :: (Fractional b0, Num b0, Ord b0) => Bool -> Ordering -> b0 -> b0 -> (b0, b0)- alphaNoConverge c cp x y = vmap2 (ifi (>= 360) (subtract 360)) (if c- then (x, y)- else if cp == GT- then (180.0, 0.0)- else if cp == LT- then (0.0, 180.0)- else let nan = 0/0- in (nan, nan))+ let alphaNoConverge c cp x y =+ vmap2 (ifi (>= 360) (subtract 360)) (if c+ then (x, y)+ else if cp == GT+ then (180.0, 0.0)+ else if cp == LT+ then (0.0, 180.0)+ else let nan = 0/0+ in (nan, nan)) in alphaNoConverge (result ed == Converge) (compare phi1 phi2) 0 0 in curve (b * a' ed * (sigma ed - deltasigma ed)) (modAzimuth alpha1) (modAzimuth alpha2) @@ -243,8 +252,9 @@ -- >>> do fr <- 27.812 ..#.. 154.295; to <- 87.7769 ..#.. 19.944; return (inverseD fr to) -- Just (GeodeticCurve 7099204.2589 Azimuth 0.0000 Azimuth 180.0000) inverseD ::- Coordinate- -> Coordinate+ (HasCoordinate c1, HasCoordinate c2) =>+ c1 -- ^ start coordinate+ -> c2 -- ^ end coordinate -> Curve inverseD = inverse wgs84 convergence@@ -267,7 +277,7 @@ ) x) => x inverse' =- optional2 inverse wgs84 convergence+ optional2 (inverse :: Ellipsoid -> Convergence -> Coordinate -> Coordinate -> Curve) wgs84 convergence ---- not exported