jord 0.3.1.0 → 0.4.0.0
raw patch · 34 files changed
+3355/−1878 lines, 34 files
Files
- ChangeLog.md +7/−0
- README.md +47/−10
- app/Eval.hs +427/−92
- app/Main.hs +183/−170
- jord.cabal +96/−90
- src/Data/Geo/Jord.hs +9/−3
- src/Data/Geo/Jord/Angle.hs +224/−224
- src/Data/Geo/Jord/AngularPosition.hs +33/−33
- src/Data/Geo/Jord/Duration.hs +144/−0
- src/Data/Geo/Jord/Earth.hs +89/−89
- src/Data/Geo/Jord/EcefPosition.hs +37/−37
- src/Data/Geo/Jord/Frames.hs +54/−7
- src/Data/Geo/Jord/Geodetics.hs +267/−226
- src/Data/Geo/Jord/Kinematics.hs +444/−0
- src/Data/Geo/Jord/Length.hs +104/−104
- src/Data/Geo/Jord/NVector.hs +25/−25
- src/Data/Geo/Jord/Parse.hs +52/−52
- src/Data/Geo/Jord/Quantity.hs +7/−7
- src/Data/Geo/Jord/Rotation.hs +2/−2
- src/Data/Geo/Jord/Speed.hs +134/−0
- src/Data/Geo/Jord/Transform.hs +0/−192
- src/Data/Geo/Jord/Transformation.hs +192/−0
- src/Data/Geo/Jord/Vector3d.hs +2/−2
- test/Data/Geo/Jord/AngleSpec.hs +92/−92
- test/Data/Geo/Jord/DurationSpec.hs +26/−0
- test/Data/Geo/Jord/EarthSpec.hs +27/−27
- test/Data/Geo/Jord/FramesSpec.hs +14/−14
- test/Data/Geo/Jord/GeodeticsSpec.hs +213/−213
- test/Data/Geo/Jord/KinematicsSpec.hs +191/−0
- test/Data/Geo/Jord/LatLongSpec.hs +63/−63
- test/Data/Geo/Jord/LengthSpec.hs +40/−35
- test/Data/Geo/Jord/SpeedSpec.hs +41/−0
- test/Data/Geo/Jord/TransformSpec.hs +0/−69
- test/Data/Geo/Jord/TransformationSpec.hs +69/−0
ChangeLog.md view
@@ -1,3 +1,10 @@+### 0.4.0.0 + +- Added ECEF, frames and delta to REPL +- Added Speed +- Added Duration +- Added Kinematics: course, position, CPA and intercept + ### 0.3.1.0 - Added ECEF position
README.md view
@@ -8,12 +8,14 @@ ## What is this? -Jord is a [Haskell](https://www.haskell.org) library that implements various geographical position calculations using the algorithms described in [Gade, K. (2010). A Non-singular Horizontal Position Representation](http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf). +Jord is a [Haskell](https://www.haskell.org) library that implements various geographical position calculations using the algorithms described in [Gade, K. (2010). A Non-singular Horizontal Position Representation](http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf) and in +[Shudde, Rex H. (1986). Some tactical algorithms for spherical geometry](https://calhoun.nps.edu/bitstream/handle/10945/29516/sometacticalalgo00shud.pdf) -- Transformation between ECEF (earth-centred, earth-fixed), Latitude/Longitude and N-Vector positions for spherical and ellipsoidal earth model -- Transformation between Latitude/Longitude and N-Vector positions -- Local, Body and North, East, Down Frames: delta between positions, target position from reference position and delta -- surface distance, initial & final bearing, interpolated position, great circle intersections, cross track distance, ... +- Transformation between ECEF (earth-centred, earth-fixed), Latitude/Longitude and N-Vector positions for spherical and ellipsoidal earth model. +- Transformation between Latitude/Longitude and N-Vector positions. +- Local, Body and North, East, Down Frames: delta between positions, target position from reference position and delta. +- Geodetics: surface distance, initial & final bearing, interpolated position, great circle intersections, cross track distance, ... +- Kinematics: position from p0, bearing and speed, closest point of approach between tracks, intercept (time, speed, min speed). ## How do I build it? @@ -34,16 +36,34 @@ let p1 = decimalLatLongHeight 1 2 (metres (-3)) let p2 = decimalLatLongHeight 4 5 (metres (-6)) let w = decimalDegrees 5 -- wander azimuth -deltaBetween p1 p2 (frameL w) wgs84 -- = deltaMetres 359490.579 302818.523 17404.272 +deltaBetween p1 p2 (frameL w) wgs84 -- deltaMetres 359490.579 302818.523 17404.272 -- destination position from 531914N0014347W having travelled 500Nm on a heading of 96.0217° --- using mean earth radius derived from the WG84 ellipsoid -destination (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) r84 +-- using mean earth radius derived from the WGS84 ellipsoid +destination84 (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) +-- using mean earth radius derived from the GRS80 ellipsoid +destination (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) r80 -- surface distance between 54°N,154°E and its antipodal position --- using mean earth radius derived from the WG84 ellipsoid let p = decimalLatLong 54 154 -surfaceDistance p (antipode p) r84 +-- using mean earth radius derived from the WGS84 ellipsoid +surfaceDistance84 p (antipode p) +-- using mean earth radius derived from the GRS80 ellipsoid +surfaceDistance p (antipode p) r80 + +-- closest point of approach between tracks +let p1 = decimalLatLong 20 (-60) +let b1 = decimalDegrees 10 +let s1 = knots 15 +let p2 = decimalLatLong 34 (-50) +let b2 = decimalDegrees 220 +let s2 = knots 300 +let t1 = Track p1 b1 s1 +let t2 = Track p2 b2 s2 +-- using mean earth radius derived from the WGS84 ellipsoid +cpa84 t1 t2 +-- using mean earth radius derived from the WGS72 ellipsoid +cpa t1 t2 r72 ``` Jord comes with a REPL (built with [haskeline](https://github.com/judah/haskeline)): @@ -52,4 +72,21 @@ $ jord-exe jord> finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E jord> angle: 126°0'0.0" (126.0) +jord> f = frameB 10d 20d 30d +jord> Body (vehicle) frame: + yaw : 10°0'0.000" (10.0) + pitch: 20°0'0.000" (20.0) + roll : 30°0'0.000" (30.0) +jord> d = delta 3000 2000 100 +jord> Delta: + x: 3000.0m <-> 3.0km <-> 1.6198704103671706nm <-> 9842.51968503937ft + y: 2000.0m <-> 2.0km <-> 1.079913606911447nm <-> 6561.679790026246ft + z: 100.0m <-> 0.1km <-> 5.399568034557235e-2nm <-> 328.0839895013123ft +jord> p0 = geo 49.66618 3.45063 0 +jord> latlong: 49°39'58.248"N,3°27'2.268"E (49.66618, 3.45063) + height : 0.0m <-> 0.0km <-> 0.0nm <-> 0.0ft +jord> target p0 f d wgs84 +jord> latlong: 49°41'30.486"N,3°28'52.561"E (49.69180166666667, 3.4812669444444446) + height : 6.0077m <-> 6.0077e-3km <-> 3.24389848812095e-3nm <-> 19.71030183727034ft +jord> ```
app/Eval.hs view
@@ -26,14 +26,8 @@ import Control.Monad.Fail import Data.Bifunctor-import Data.Geo.Jord.Angle-import Data.Geo.Jord.AngularPosition-import Data.Geo.Jord.Geodetics-import Data.Geo.Jord.LatLong-import Data.Geo.Jord.Length-import Data.Geo.Jord.NVector-import Data.Geo.Jord.Quantity-import Data.Geo.Jord.Transform+import Data.Either (rights)+import Data.Geo.Jord import Data.List hiding (delete, insert, lookup) import Data.Maybe import Prelude hiding (fail, lookup)@@ -44,14 +38,125 @@ data Value = Ang Angle -- ^ angle | Bool Bool -- ^ boolean + | Cpa (Cpa (AngularPosition NVector)) -- ^ CPA + | Dlt Delta -- ^ delta + | Dur Duration -- ^ duration | Double Double -- ^ double - | Len Length -- ^ length + | Ep EcefPosition -- ^ ECEF position + | Em Earth -- ^ earth model + | FrmB Angle+ Angle+ Angle -- ^ yaw, pitch and roll of Body frame + | FrmL Angle -- ^ wander azimuth of Local frame + | FrmN -- ^ North, east, down frame | Gc GreatCircle -- ^ great circle - | Geo (AngularPosition LatLong) -- ^ latitude, longitude and height - | NVec (AngularPosition NVector) -- ^ n-vector and height + | Gp (AngularPosition LatLong) -- ^ latitude, longitude and height + | Intp (Intercept (AngularPosition NVector)) -- ^ Intercept + | Len Length -- ^ length + | Ned Ned -- ^ north east down + | Np (AngularPosition NVector) -- ^ n-vector and height + | Spd Speed -- ^ speed + | Trk (Track (AngularPosition NVector)) -- ^ track | Vals [Value] -- array of values - deriving (Eq, Show) +-- | show value. +instance Show Value where+ show (Ang a) = "angle: " ++ showAng a+ show (Bool b) = show b+ show (Cpa c) =+ "closest point of approach:" +++ "\n time : " +++ show (cpaTime c) +++ "\n distance: " +++ showLen (cpaDistance c) +++ "\n pos1 : " +++ showLl (fromNVector . cpaPosition1 $ c :: LatLong) +++ "\n pos2 : " ++ showLl (fromNVector . cpaPosition2 $ c :: LatLong)+ show (Dlt d) =+ "Delta:" +++ "\n x: " +++ showLen (dx d) ++ "\n y: " ++ showLen (dy d) ++ "\n z: " ++ showLen (dz d)+ show (Dur d) = "duration: " ++ show d+ show (Double d) = show d+ show (Em m) = "Earth model: " ++ show m+ show (Ep p) =+ "ECEF:" +++ "\n x: " +++ showLen (ex p) ++ "\n y: " ++ showLen (ey p) ++ "\n z: " ++ showLen (ez p)+ show (FrmB y p r) =+ "Body (vehicle) frame:" +++ "\n yaw : " +++ showAng y ++ "\n pitch: " ++ showAng p ++ "\n roll : " ++ showAng r+ show (FrmL w) = "Local frame:" ++ "\n wander azimuth: " ++ showAng w+ show FrmN = "North, East, Down frame"+ show (Gc gc) = "great circle: " ++ show gc+ show (Gp g) = "latlong: " ++ showLl ll ++ "\n height : " ++ showLen h+ where+ ll = pos g+ h = height g+ show (Intp i) =+ "intercept:" +++ "\n time : " +++ show (interceptTime i) +++ "\n distance : " +++ showLen (interceptDistance i) +++ "\n pos : " +++ showLl (fromNVector . interceptPosition $ i :: LatLong) +++ "\n interceptor speed : " +++ showSpd (interceptorSpeed i) +++ "\n interceptor bearing: " ++ showAng (interceptorBearing i)+ show (Len l) = "length: " ++ showLen l+ show (Ned d) =+ "NED:" +++ "\n north: " +++ showLen (north d) +++ "\n east : " ++ showLen (east d) ++ "\n down : " ++ showLen (down d)+ show (Np nv) =+ "n-vector: " +++ show x ++ ", " ++ show y ++ ", " ++ show z ++ "\n height : " ++ showLen h+ where+ v = vec (pos nv)+ x = vx v+ y = vy v+ z = vz v+ h = height nv+ show (Trk t) =+ "track:" +++ "\n position: " +++ showLl (fromNVector . trackPos $ t :: LatLong) +++ "\n height : " +++ showLen (height . trackPos $ t) +++ "\n bearing : " +++ showAng (trackBearing t) ++ "\n speed : " ++ showSpd (trackSpeed t)+ show (Spd s) = "speed: " ++ showSpd s+ show (Vals []) = "empty"+ show (Vals vs) = "\n " ++ intercalate "\n\n " (map show vs)++showAng :: Angle -> String+showAng a = show a ++ " (" ++ show (toDecimalDegrees a) ++ ")"++showLl :: LatLong -> String+showLl ll =+ show ll +++ " (" +++ show (toDecimalDegrees (latitude ll)) ++ ", " ++ show (toDecimalDegrees (longitude ll)) ++ ")"++showLen :: Length -> String+showLen l =+ show (toMetres l) +++ "m <-> " +++ show (toKilometres l) +++ "km <-> " ++ show (toNauticalMiles l) ++ "nm <-> " ++ show (toFeet l) ++ "ft"++showSpd :: Speed -> String+showSpd s =+ show (toKilometresPerHour s) +++ "km/h <-> " +++ show (toMetresPerSecond s) +++ "m/s <-> " +++ show (toKnots s) +++ "kt <-> " ++ show (toMilesPerHour s) ++ "mph <-> " ++ show (toFeetPerSecond s) ++ "ft/s"+ -- | 'Either' an error or a 'Value'. type Result = Either String Value @@ -93,7 +198,7 @@ -- -- @ -- dest = eval "destination 54°N,154°E 54° 1000m" -- Right Ll --- dest = eval "toNVector (destination 54°N,154°E 54° 1000m)" -- Right NVec +-- dest = eval "toNVector (destination 54°N,154°E 54° 1000m)" -- Right Np -- @ -- -- Every function call must be wrapped between parentheses, however they can be ommitted for the top level call. @@ -115,12 +220,12 @@ convert r True = r convert r False = case r of- Right v@(NVec _) -> Right (toGeo v)+ Right v@(Np _) -> Right (toGeo v) Right (Vals vs) -> Right (Vals (map toGeo vs)) oth -> oth toGeo :: Value -> Value-toGeo (NVec v) = Geo (fromNVector v)+toGeo (Np v) = Gp (fromNVector v) toGeo val = val -- | All supported functions. @@ -128,19 +233,34 @@ functions = [ "antipode" , "crossTrackDistance"+ , "cpa"+ , "delta"+ , "deltaBetween" , "destination"+ , "ecef"+ , "frameB"+ , "frameL"+ , "frameN" , "finalBearing"- , "geoPos"+ , "fromEcef"+ , "geo" , "greatCircle" , "initialBearing"+ , "intercept"+ , "interceptBySpeed"+ , "interceptByTime" , "interpolate" , "intersections" , "insideSurface" , "mean"+ , "ned"+ , "nedBetween"+ , "position" , "surfaceDistance"- , "toKilometres"- , "toMetres"- , "toNauticalMiles"+ , "target"+ , "targetN"+ , "track"+ , "toEcef" , "toNVector" ] @@ -162,7 +282,7 @@ expr s = do ts <- tokenise s ast <- parse ts- fmap (expectVec ts,) (transform ast)+ fmap (expectVec ts, ) (transform ast) expectVec :: [Token] -> Bool expectVec (_:Func "toNVector":_) = True@@ -171,131 +291,210 @@ evalExpr :: Expr -> Vault -> Result evalExpr (Param p) vault = case lookup p vault of- Just (Geo geo) -> Right (NVec (toNVector geo))+ Just (Gp geo) -> Right (Np (toNVector geo)) Just v -> Right v Nothing -> tryRead p evalExpr (Antipode a) vault = case evalExpr a vault of- (Right (NVec p)) -> Right (NVec (antipode p))+ (Right (Np p)) -> Right (Np (antipode p)) r -> Left ("Call error: antipode " ++ showErr [r])+evalExpr (ClosestPointOfApproach a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (Trk t1), Right (Trk t2)] ->+ maybe (Left "closest point of approach in the past") (Right . Cpa) (cpa84 t1 t2)+ r -> Left ("Call error: cpa " ++ showErr r) evalExpr (CrossTrackDistance a b) vault = case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p), Right (Gc gc)] -> Right (Len (crossTrackDistance84 p gc))+ [Right (Np p), Right (Gc gc)] -> Right (Len (crossTrackDistance84 p gc)) r -> Left ("Call error: crossTrackDistance " ++ showErr r)+evalExpr (DeltaBetween a b c d) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault, evalEarth d] of+ [Right (Np p1), Right (Np p2), Right (FrmB y p r), Right (Em m)] ->+ Right (Dlt (deltaBetween p1 p2 (frameB y p r) m))+ [Right (Np p1), Right (Np p2), Right (FrmL w), Right (Em m)] ->+ Right (Dlt (deltaBetween p1 p2 (frameL w) m))+ [Right (Np p1), Right (Np p2), Right FrmN, Right (Em m)] ->+ Right (Dlt (deltaBetween p1 p2 frameN m))+ r -> Left ("Call error: deltaBetween " ++ showErr r)+evalExpr (DeltaV a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Len x), Right (Len y), Right (Len z)] -> Right (Dlt (delta x y z))+ [Right (Double x), Right (Double y), Right (Double z)] -> Right (Dlt (deltaMetres x y z))+ r -> Left ("Call error: delta " ++ showErr r) evalExpr (Destination a b c) vault = case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of- [Right (NVec p), Right (Ang a'), Right (Len l)] -> Right (NVec (destination84 p a' l))- [Right (NVec p), Right (Double a'), Right (Len l)] ->- Right (NVec (destination84 p (decimalDegrees a') l))+ [Right (Np p), Right (Ang a'), Right (Len l)] -> Right (Np (destination84 p a' l))+ [Right (Np p), Right (Double a'), Right (Len l)] ->+ Right (Np (destination84 p (decimalDegrees a') l)) r -> Left ("Call error: destination " ++ showErr r)+evalExpr (Ecef a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Len x), Right (Len y), Right (Len z)] -> Right (Ep (ecef x y z))+ [Right (Double x), Right (Double y), Right (Double z)] -> Right (Ep (ecefMetres x y z))+ r -> Left ("Call error: ecef " ++ showErr r)+evalExpr (FrameB a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Ang a'), Right (Ang b'), Right (Ang c')] -> Right (FrmB a' b' c')+ r -> Left ("Call error: frameB " ++ showErr r)+evalExpr (FrameL a) vault =+ case evalExpr a vault of+ (Right (Ang a')) -> Right (FrmL a')+ r -> Left ("Call error: frameL " ++ showErr [r])+evalExpr FrameN _ = Right FrmN+evalExpr (FromEcef a b) vault =+ case [evalExpr a vault, evalEarth b] of+ [Right (Ep p), Right (Em m)] -> Right (Np (fromEcef p m))+ r -> Left ("Call error: fromEcef " ++ showErr r) evalExpr (FinalBearing a b) vault = case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p1), Right (NVec p2)] ->+ [Right (Np p1), Right (Np p2)] -> maybe (Left "Call error: finalBearing identical points") (Right . Ang) (finalBearing p1 p2) r -> Left ("Call error: finalBearing " ++ showErr r)-evalExpr (GeoPos as) vault =+evalExpr (Geo as) vault = case vs of- [Right p@(NVec _)] -> Right p- [Right (NVec v), Right (Len h)] -> Right (NVec (AngularPosition (pos v) h))+ [Right p@(Np _)] -> Right p+ [Right (Np v), Right (Len h)] -> Right (Np (AngularPosition (pos v) h)) [Right (Double lat), Right (Double lon)] -> bimap- (\e -> "Call error: geoPos : " ++ e)- (NVec . toNVector)+ (\e -> "Call error: geo " ++ e)+ (Np . toNVector) (decimalLatLongHeightE lat lon zero) [Right (Double lat), Right (Double lon), Right (Len h)] ->- bimap- (\e -> "Call error: geoPos : " ++ e)- (NVec . toNVector)- (decimalLatLongHeightE lat lon h)+ bimap (\e -> "Call error: geo " ++ e) (Np . toNVector) (decimalLatLongHeightE lat lon h) [Right (Double lat), Right (Double lon), Right (Double h)] -> bimap- (\e -> "Call error: geoPos : " ++ e)- (NVec . toNVector)+ (\e -> "Call error: geo " ++ e)+ (Np . toNVector) (decimalLatLongHeightE lat lon (metres h))- r -> Left ("Call error: geoPos " ++ showErr r)+ r -> Left ("Call error: geo " ++ showErr r) where vs = map (`evalExpr` vault) as-evalExpr (GreatCircleSC a b) vault =- case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p1), Right (NVec p2)] -> bimap id Gc (greatCircleE p1 p2)- [Right (NVec p), Right (Ang a')] -> Right (Gc (greatCircleBearing p a'))+evalExpr (GreatCircleE as) vault =+ case fmap (`evalExpr` vault) as of+ [Right (Np p1), Right (Np p2)] -> bimap id Gc (greatCircleE (p1, p2))+ [Right (Np p), Right (Ang a')] -> bimap id Gc (greatCircleE (p, a'))+ [Right (Trk t)] -> bimap id Gc (greatCircleE t) r -> Left ("Call error: greatCircle " ++ showErr r) evalExpr (InitialBearing a b) vault = case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p1), Right (NVec p2)] ->+ [Right (Np p1), Right (Np p2)] -> maybe (Left "Call error: initialBearing identical points") (Right . Ang) (initialBearing p1 p2) r -> Left ("Call error: initialBearing " ++ showErr r)+evalExpr (Intercept a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (Trk t), Right (Np i)] ->+ maybe (Left "intercept impossible") (Right . Intp) (intercept84 t i)+ r -> Left ("Call error: intercept " ++ showErr r)+evalExpr (InterceptBySpeed a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Trk t), Right (Np i), Right (Spd s)] ->+ maybe (Left "intercept impossible") (Right . Intp) (interceptBySpeed84 t i s)+ r -> Left ("Call error: interceptBySpeed " ++ showErr r)+evalExpr (InterceptByTime a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Trk t), Right (Np i), Right (Dur d)] ->+ maybe (Left "intercept impossible") (Right . Intp) (interceptByTime84 t i d)+ r -> Left ("Call error: interceptByTime " ++ showErr r) evalExpr (Interpolate a b c) vault = case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p1), Right (NVec p2)] -> Right (NVec (interpolate p1 p2 c))+ [Right (Np p1), Right (Np p2)] -> Right (Np (interpolate p1 p2 c)) r -> Left ("Call error: interpolate " ++ showErr r) evalExpr (Intersections a b) vault = case [evalExpr a vault, evalExpr b vault] of [Right (Gc gc1), Right (Gc gc2)] -> maybe (Right (Vals []))- (\is -> Right (Vals [NVec (fst is), NVec (snd is)]))+ (\is -> Right (Vals [Np (fst is), Np (snd is)])) (intersections gc1 gc2 :: Maybe (AngularPosition NVector, AngularPosition NVector)) r -> Left ("Call error: intersections " ++ showErr r) evalExpr (InsideSurface as) vault = let m = map (`evalExpr` vault) as- ps = [p | Right (NVec p) <- m]+ ps = [p | Right (Np p) <- m] in if length m == length ps && length ps > 3 then Right (Bool (insideSurface (head ps) (tail ps))) else Left ("Call error: insideSurface " ++ showErr m) evalExpr (Mean as) vault = let m = map (`evalExpr` vault) as- ps = [p | Right (NVec p) <- m]+ ps = [p | Right (Np p) <- m] in if length m == length ps- then maybe (Left ("Call error: mean " ++ showErr m)) (Right . NVec) (mean ps)+ then maybe (Left ("Call error: mean " ++ showErr m)) (Right . Np) (mean ps) else Left ("Call error: mean " ++ showErr m)+evalExpr (NedBetween a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalEarth c] of+ [Right (Np p1), Right (Np p2), Right (Em m)] -> Right (Ned (nedBetween p1 p2 m))+ r -> Left ("Call error: nedBetween " ++ showErr r)+evalExpr (NedV a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Len x), Right (Len y), Right (Len z)] -> Right (Ned (ned x y z))+ [Right (Double x), Right (Double y), Right (Double z)] -> Right (Ned (nedMetres x y z))+ r -> Left ("Call error: ned " ++ showErr r)+evalExpr (Position a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (Trk t), Right (Dur d)] -> Right (Np (position84 t d))+ r -> Left ("Call error: position " ++ showErr r) evalExpr (SurfaceDistance a b) vault = case [evalExpr a vault, evalExpr b vault] of- [Right (NVec p1), Right (NVec p2)] -> Right (Len (surfaceDistance84 p1 p2))+ [Right (Np p1), Right (Np p2)] -> Right (Len (surfaceDistance84 p1 p2)) r -> Left ("Call error: surfaceDistance " ++ showErr r)-evalExpr (ToKilometres e) vault =- case evalExpr e vault of- (Right (Len l)) -> Right (Double (toKilometres l))- r -> Left ("Call error: toKilometres" ++ showErr [r])-evalExpr (ToMetres e) vault =- case evalExpr e vault of- (Right (Len l)) -> Right (Double (toMetres l))- r -> Left ("Call error: toMetres" ++ showErr [r])-evalExpr (ToNauticalMiles e) vault =- case evalExpr e vault of- (Right (Len l)) -> Right (Double (toNauticalMiles l))- r -> Left ("Call error: toNauticalMiles" ++ showErr [r])+evalExpr (Target a b c d) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault, evalEarth d] of+ [Right (Np p0), Right (FrmB y p r), Right (Dlt d'), Right (Em m)] ->+ Right (Np (target p0 (frameB y p r) d' m))+ [Right (Np p0), Right (FrmL w), Right (Dlt d'), Right (Em m)] ->+ Right (Np (target p0 (frameL w) d' m))+ [Right (Np p0), Right FrmN, Right (Dlt d'), Right (Em m)] ->+ Right (Np (target p0 frameN d' m))+ r -> Left ("Call error: target " ++ showErr r)+evalExpr (TargetN a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalEarth c] of+ [Right (Np p0), Right (Ned d), Right (Em m)] -> Right (Np (targetN p0 d m))+ r -> Left ("Call error: targetN " ++ showErr r)+evalExpr (TrackE a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (Np p), Right (Ang b'), Right (Spd s)] -> Right (Trk (Track p b' s))+ r -> Left ("Call error: track " ++ showErr r)+evalExpr (ToEcef a b) vault =+ case [evalExpr a vault, evalEarth b] of+ [Right (Np p), Right (Em m)] -> Right (Ep (toEcef p m))+ r -> Left ("Call error: toEcef " ++ showErr r) evalExpr (ToNVector a) vault = case evalExpr a vault of- r@(Right (NVec _)) -> r+ r@(Right (Np _)) -> r r -> Left ("Call error: toNVector " ++ showErr [r]) +evalEarth :: String -> Result+evalEarth "wgs84" = Right (Em wgs84)+evalEarth "grs80" = Right (Em grs80)+evalEarth "wgs72" = Right (Em wgs72)+evalEarth "s84" = Right (Em s84)+evalEarth "s80" = Right (Em s80)+evalEarth "s72" = Right (Em s72)+evalEarth s = Left s+ showErr :: [Result] -> String showErr rs = " > " ++ intercalate " & " (map (either id show) rs) tryRead :: String -> Result-tryRead s =- case r of- [a@(Right (Ang _)), _, _, _] -> a- [_, l@(Right (Len _)), _, _] -> l- [_, _, Right (Geo geo), _] -> Right (NVec (toNVector geo))- [_, _, _, d@(Right (Double _))] -> d- _ -> Left ("couldn't read " ++ s)+tryRead s+ | null r = Left ("couldn't read " ++ s)+ | otherwise = Right (head r) where r =- map- ($ s)- [ readE readAngleE Ang- , readE readLengthE Len- , readE readLatLongE (\ll -> Geo (AngularPosition ll zero))- , readE readEither Double- ]+ rights+ (map ($ s)+ [ readE readAngleE Ang+ , readE readLengthE Len+ , readE readSpeedE Spd+ , readE readDurationE Dur+ , readE readLatLongE (\ll -> Np (toNVector (AngularPosition ll zero)))+ , readE readEither Double+ ]) readE :: (String -> Either String a) -> (a -> Value) -> String -> Either String Value readE p v s = bimap id v (p s)@@ -399,18 +598,44 @@ data Expr = Param String | Antipode Expr+ | ClosestPointOfApproach Expr+ Expr | CrossTrackDistance Expr Expr+ | DeltaBetween Expr+ Expr+ Expr+ String+ | DeltaV Expr+ Expr+ Expr | Destination Expr Expr Expr+ | Ecef Expr+ Expr+ Expr+ | FrameB Expr+ Expr+ Expr+ | FrameL Expr+ | FrameN | FinalBearing Expr Expr- | GeoPos [Expr]- | GreatCircleSC Expr- Expr+ | FromEcef Expr+ String+ | Geo [Expr]+ | GreatCircleE [Expr] | InitialBearing Expr Expr+ | Intercept Expr+ Expr+ | InterceptBySpeed Expr+ Expr+ Expr+ | InterceptByTime Expr+ Expr+ Expr | Interpolate Expr Expr Double@@ -418,40 +643,107 @@ Expr | InsideSurface [Expr] | Mean [Expr]+ | NedBetween Expr+ Expr+ String+ | NedV Expr+ Expr+ Expr+ | Position Expr+ Expr | SurfaceDistance Expr Expr- | ToKilometres Expr- | ToMetres Expr- | ToNauticalMiles Expr+ | Target Expr+ Expr+ Expr+ String+ | TargetN Expr+ Expr+ String+ | TrackE Expr+ Expr+ Expr+ | ToEcef Expr+ String | ToNVector Expr deriving (Show) transform :: (MonadFail m) => Ast -> m Expr transform (Call "antipode" [e]) = fmap Antipode (transform e)+transform (Call "cpa" [e1, e2]) = do+ t1 <- transform e1+ t2 <- transform e2+ return (ClosestPointOfApproach t1 t2) transform (Call "crossTrackDistance" [e1, e2]) = do p <- transform e1 gc <- transform e2 return (CrossTrackDistance p gc)+transform (Call "delta" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ p3 <- transform e3+ return (DeltaV p1 p2 p3)+transform (Call "deltaBetween" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ f <- transform e3+ return (DeltaBetween p1 p2 f "wgs84")+transform (Call "deltaBetween" [e1, e2, e3, Lit s]) = do+ p1 <- transform e1+ p2 <- transform e2+ f <- transform e3+ return (DeltaBetween p1 p2 f s) transform (Call "destination" [e1, e2, e3]) = do p1 <- transform e1 p2 <- transform e2 p3 <- transform e3 return (Destination p1 p2 p3)+transform (Call "ecef" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ p3 <- transform e3+ return (Ecef p1 p2 p3)+transform (Call "frameB" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ p3 <- transform e3+ return (FrameB p1 p2 p3)+transform (Call "frameL" [e]) = fmap FrameL (transform e)+transform (Call "frameN" []) = return FrameN+transform (Call "fromEcef" [e]) = do+ p <- transform e+ return (FromEcef p "wgs84")+transform (Call "fromEcef" [e, Lit s]) = do+ p <- transform e+ return (FromEcef p s) transform (Call "finalBearing" [e1, e2]) = do p1 <- transform e1 p2 <- transform e2 return (FinalBearing p1 p2)-transform (Call "geoPos" e) = do+transform (Call "geo" e) = do ps <- mapM transform e- return (GeoPos ps)-transform (Call "greatCircle" [e1, e2]) = do- p1 <- transform e1- p2 <- transform e2- return (GreatCircleSC p1 p2)+ return (Geo ps)+transform (Call "greatCircle" e) = do+ ps <- mapM transform e+ return (GreatCircleE ps) transform (Call "initialBearing" [e1, e2]) = do p1 <- transform e1 p2 <- transform e2 return (InitialBearing p1 p2)+transform (Call "intercept" [e1, e2]) = do+ t <- transform e1+ i <- transform e2+ return (Intercept t i)+transform (Call "interceptBySpeed" [e1, e2, e3]) = do+ t <- transform e1+ i <- transform e2+ s <- transform e3+ return (InterceptBySpeed t i s)+transform (Call "interceptTime" [e1, e2, e3]) = do+ t <- transform e1+ i <- transform e2+ d <- transform e3+ return (InterceptByTime t i d) transform (Call "interpolate" [e1, e2, Lit s]) = do p1 <- transform e1 p2 <- transform e2@@ -469,13 +761,56 @@ transform (Call "mean" e) = do ps <- mapM transform e return (Mean ps)+transform (Call "ned" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ p3 <- transform e3+ return (NedV p1 p2 p3)+transform (Call "nedBetween" [e1, e2]) = do+ p1 <- transform e1+ p2 <- transform e2+ return (NedBetween p1 p2 "wgs84")+transform (Call "nedBetween" [e1, e2, Lit s]) = do+ p1 <- transform e1+ p2 <- transform e2+ return (NedBetween p1 p2 s)+transform (Call "position" [e1, e2]) = do+ t <- transform e1+ d <- transform e2+ return (Position t d) transform (Call "surfaceDistance" [e1, e2]) = do p1 <- transform e1 p2 <- transform e2 return (SurfaceDistance p1 p2)-transform (Call "toKilometres" [e]) = fmap ToKilometres (transform e)-transform (Call "toMetres" [e]) = fmap ToMetres (transform e)-transform (Call "toNauticalMiles" [e]) = fmap ToNauticalMiles (transform e)+transform (Call "target" [e1, e2, e3]) = do+ p0 <- transform e1+ f <- transform e2+ d <- transform e3+ return (Target p0 f d "wgs84")+transform (Call "target" [e1, e2, e3, Lit s]) = do+ p0 <- transform e1+ f <- transform e2+ d <- transform e3+ return (Target p0 f d s)+transform (Call "targetN" [e1, e2]) = do+ p0 <- transform e1+ d <- transform e2+ return (TargetN p0 d "wgs84")+transform (Call "targetN" [e1, e2, Lit s]) = do+ p0 <- transform e1+ d <- transform e2+ return (TargetN p0 d s)+transform (Call "track" [e1, e2, e3]) = do+ p0 <- transform e1+ b <- transform e2+ s <- transform e3+ return (TrackE p0 b s)+transform (Call "toEcef" [e]) = do+ p <- transform e+ return (ToEcef p "wgs84")+transform (Call "toEcef" [e, Lit s]) = do+ p <- transform e+ return (ToEcef p s) transform (Call "toNVector" [e]) = fmap ToNVector (transform e) transform (Call f e) = fail ("Semantic error: " ++ f ++ " does not accept " ++ show e) transform (Lit s) = return (Param s)
app/Main.hs view
@@ -8,175 +8,188 @@ -- -- REPL around "Jord". -- -module Main where--import Data.Geo.Jord-import Data.List ((\\), dropWhileEnd, intercalate, isPrefixOf)-import Eval-import Prelude hiding (lookup)-import System.Console.Haskeline--search :: String -> [Completion]-search s = map simpleCompletion $ filterFunc s--filterFunc :: String -> [String]-filterFunc s = map (\f -> pref ++ f) filtered- where+module Main where + +import Data.Geo.Jord +import Data.List ((\\), dropWhileEnd, isPrefixOf) +import Eval +import Prelude hiding (lookup) +import System.Console.Haskeline + +search :: String -> [Completion] +search s = map simpleCompletion $ filterFunc s + +filterFunc :: String -> [String] +filterFunc s = map (\f -> pref ++ f) filtered + where pref = dropWhileEnd (/= '(') s -- everything before the last '(' inclusive func = (\\) s pref -- everything after the last '(' - filtered = filter (\f -> func `isPrefixOf` f) functions--mySettings :: Settings IO-mySettings =- Settings- { complete = completeWord Nothing " \t" $ return . search- , historyFile = Nothing- , autoAddHistory = True- }--main :: IO ()-main = do- putStrLn- ("jord interpreter, version " ++- jordVersion ++ ": https://github.com/ofmooseandmen/jord :? for help")- runInputT mySettings $ withInterrupt $ loop emptyVault- where- loop state = do- input <- handleInterrupt (return (Just "")) $ getInputLine "jord> "- case input of- Nothing -> return ()- Just ":quit" -> return ()- Just ":q" -> return ()- Just i -> do- let (result, newState) = evalS i state- printS result- loop newState--printS :: Either String String -> InputT IO ()-printS (Left err) = outputStrLn ("jord> " ++ err)-printS (Right "") = return ()-printS (Right r) = outputStrLn ("jord> " ++ r)--evalS :: String -> Vault -> (Either String String, Vault)-evalS s vault- | null s = (Right "", vault)- | head s == ':' = evalC w vault- | (v:"=":e) <- w =- let r = eval (unwords e) vault- vault' = save r v vault- in (showR r, vault')- | otherwise = (showR (eval s vault), vault)- where- w = words s--evalC :: [String] -> Vault -> (Either String String, Vault)-evalC [":show", v] vault = (evalShow v vault, vault)-evalC [":delete", v] vault = evalDel (Just v) vault-evalC [":clear"] vault = evalDel Nothing vault-evalC [":help"] vault = (Right help, vault)-evalC [":?"] vault = (Right help, vault)-evalC c vault = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), vault)--evalShow :: String -> Vault -> Either String String-evalShow n vault = maybe (Left ("Unbound variable: " ++ n)) (Right . showVar n) (lookup n vault)--evalDel :: Maybe String -> Vault -> (Either String String, Vault)-evalDel (Just n) vault = (Right ("deleted var: " ++ n), delete n vault)-evalDel Nothing _ = (Right "deleted all variable ", emptyVault)--help :: String-help =- "\njord interpreter, version " ++- jordVersion ++- "\n" ++- "\n Commands available from the prompt:\n\n" ++- " :help, :? display this list of commands\n" ++- " :quit, :q quit jord\n" ++- " :show {var} shows {var}\n" ++- " :delete {var} deletes {var}\n" ++- " :clear deletes all variable(s)\n" ++- "\n Jord expressions:\n\n" ++- " (f x y) where f is one of function described below and x and y\n" ++- " are either parameters in one of the format described below or\n" ++- " a call to another function\n" ++- "\n" ++- " (finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E)\n" ++- "\n" ++- " Top level () can be ommitted: antipode 54N028E\n" ++- "\n Position calculations (Spherical Earth):\n\n" ++- " The following calculations assume a spherical earth model with a radius\n" ++- " derived from the WGS84 ellipsoid: " ++- show r84 ++- "\n" ++- "\n antipode pos antipodal point of pos\n" ++- " crossTrackDistance pos gc signed distance from pos to great circle gc\n" ++- " destination pos ang len destination position from pos having travelled len\n" ++- " on initial bearing ang (either in text form or decimal degrees)\n" ++- " finalBearing pos1 pos2 initial bearing from pos1 to pos2\n" ++- " initialBearing pos1 pos2 bearing arriving at pos2 from pos1\n" ++- " interpolate pos1 pos2 (0..1) position at fraction between pos1 and pos2\n" ++- " intersections gc1 gc2 intersections between great circles gc1 and gc2\n" ++- " exactly 0 or 2 intersections\n" ++- " insideSurface pos [pos] is p inside surface polygon?\n" ++- " mean [pos] geographical mean surface position of [pos]\n" ++- " surfaceDistance pos1 pos2 surface distance between pos1 and pos2\n" ++- "\n Constructors and conversions:\n\n" ++- " geoPos latlong surface geographic position from latlong\n" ++- " geoPos latlong height geographic position from latlong and height\n" ++- " geoPos lat long height geographic position from decimal latitude, longitude and height\n" ++- " geoPos lat long metres geographic position from decimal latitude, longitude and metres\n" ++- " greatCircle pos1 pos2 great circle passing by pos1 and pos2\n" ++- " greatCircle pos ang great circle passing by pos and heading on bearing ang\n" ++- " toKilometres len length to kilometres\n" ++- " toMetres len length to metres\n" ++- " toNauticalMiles len length to nautical miles\n" ++- " toNVector pos n-vector corresponding to pos\n" ++- "\n Supported Lat/Long formats:\n\n" ++- " DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - 553621N0130209E\n" ++- " d°m's\"[N|S],d°m's\"[E|W] - 55°36'21\"N,13°2'9\"E\n" ++- " ^ zeroes can be ommitted and separtors can also be d, m, s\n" ++- " decimal°[N|S],decimal°[E|W] - 51.885°N,13,1°E\n" ++- "\n Supported Angle formats:\n\n" ++- " d°m's - 55°36'21.154\n" ++- " decimal° - 51.885°\n" ++- "\n Supported Length formats: {l}m, {l}km, {l}Nm\n\n" ++- "\n Every evaluated result can be saved by prefixing the expression with \"{var} = \"\n" ++- " Saved results can subsequently be used when calling a function\n" ++- " jord> a = antipode 54N028E\n" ++ " jord> antipode a\n"--save :: Result -> String -> Vault -> Vault-save (Right v) k vault = insert k v vault-save _ _ vault = vault--showR :: Result -> Either String String-showR (Left err) = Left err-showR (Right v) = Right (showV v)--showV :: Value -> String-showV (Ang a) = "angle: " ++ show a ++ " (" ++ show (toDecimalDegrees a) ++ ")"-showV (Bool b) = show b-showV (Double d) = show d-showV (Gc gc) = "great circle: " ++ show gc-showV (Len l) = "length: " ++ show l-showV (Geo g) =- "latitude, longitude: " ++- show ll ++- " (" ++- show (toDecimalDegrees (latitude ll)) ++- ", " ++ show (toDecimalDegrees (longitude ll)) ++ ") - height: " ++ show h- where- ll = pos g- h = height g-showV (NVec nv) =- "n-vector: (" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ") - height: " ++ show h- where- v = vec (pos nv)- x = vx v- y = vy v- z = vz v- h = height nv-showV (Vals []) = "empty"-showV (Vals vs) = "\n " ++ intercalate "\n " (map showV vs)--showVar :: String -> Value -> String-showVar n v = n ++ "=" ++ showV v+ filtered = filter (\f -> func `isPrefixOf` f) functions + +mySettings :: Settings IO +mySettings = + Settings + { complete = completeWord Nothing " \t" $ return . search + , historyFile = Nothing + , autoAddHistory = True + } + +main :: IO () +main = do + putStrLn + ("jord interpreter, version " ++ + jordVersion ++ ": https://github.com/ofmooseandmen/jord :? for help") + runInputT mySettings $ withInterrupt $ loop emptyVault + where + loop state = do + input <- handleInterrupt (return (Just "")) $ getInputLine "jord> " + case input of + Nothing -> return () + Just ":quit" -> return () + Just ":q" -> return () + Just i -> do + let (result, newState) = evalS i state + printS result + loop newState + +printS :: Either String String -> InputT IO () +printS (Left err) = outputStrLn ("jord> " ++ err) +printS (Right "") = return () +printS (Right r) = outputStrLn ("jord> " ++ r) + +evalS :: String -> Vault -> (Either String String, Vault) +evalS s vault + | null s = (Right "", vault) + | head s == ':' = evalC w vault + | (v:"=":e) <- w = + let r = eval (unwords e) vault + vault' = save r v vault + in (showR r, vault') + | otherwise = (showR (eval s vault), vault) + where + w = words s + +evalC :: [String] -> Vault -> (Either String String, Vault) +evalC [":show", v] vault = (evalShow v vault, vault) +evalC [":delete", v] vault = evalDel (Just v) vault +evalC [":clear"] vault = evalDel Nothing vault +evalC [":help"] vault = (Right help, vault) +evalC [":?"] vault = (Right help, vault) +evalC c vault = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), vault) + +evalShow :: String -> Vault -> Either String String +evalShow n vault = maybe (Left ("Unbound variable: " ++ n)) (Right . showVar n) (lookup n vault) + +evalDel :: Maybe String -> Vault -> (Either String String, Vault) +evalDel (Just n) vault = (Right ("deleted var: " ++ n), delete n vault) +evalDel Nothing _ = (Right "deleted all variable ", emptyVault) + +help :: String +help = + "\njord interpreter, version " ++ + jordVersion ++ + "\n" ++ + "\n Commands available from the prompt:\n\n" ++ + " :help, :? display this list of commands\n" ++ + " :quit, :q quit jord\n" ++ + " :show {var} shows {var}\n" ++ + " :delete {var} deletes {var}\n" ++ + " :clear deletes all variable(s)\n" ++ + "\n Jord expressions:\n\n" ++ + " (f x y) where f is one of function described below and x and y\n" ++ + " are either parameters in one of the format described below or\n" ++ + " a call to another function\n" ++ + "\n" ++ + " (finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E)\n" ++ + "\n" ++ + " Top level () can be ommitted: antipode 54N028E\n" ++ + "\n Position calculations (Frames):\n\n" ++ + " The following calculations work with both ellipsoidal and derived earth model\n" ++ + " WGS84 ellipsoid is used if model is omitted\n" ++ + "\n deltaBetween pos1 pos2 frame (earth) delta between pos1 and pos2 in frame originating at pos1\n" ++ + " nedBetween pos1 pos2 (earth) NED between pos1 and pos2 in frame N originating at pos1\n" ++ + " target pos frame delta (earth) target position from pos and delta in frame originating at pos\n" ++ + " targetN pos delta (earth) target position from pos and NED in frame N originating at pos\n" ++ + "\n Position calculations (Spherical Earth):\n\n" ++ + " The following calculations assume a spherical earth model with a radius\n" ++ + " derived from the WGS84 ellipsoid: " ++ + show r84 ++ + "\n" ++ + "\n antipode pos antipodal point of pos\n" ++ + " crossTrackDistance pos gc signed distance from pos to great circle gc\n" ++ + " destination pos ang len destination position from pos having travelled len\n" ++ + " on initial bearing ang (either in text form or decimal degrees)\n" ++ + " finalBearing pos1 pos2 initial bearing from pos1 to pos2\n" ++ + " initialBearing pos1 pos2 bearing arriving at pos2 from pos1\n" ++ + " interpolate pos1 pos2 (0..1) position at fraction between pos1 and pos2\n" ++ + " intersections gc1 gc2 intersections between great circles gc1 and gc2\n" ++ + " exactly 0 or 2 intersections\n" ++ + " insideSurface pos [pos] is p inside surface polygon?\n" ++ + " mean [pos] geographical mean surface position of [pos]\n" ++ + " surfaceDistance pos1 pos2 surface distance between pos1 and pos2\n" ++ + "\n Kinematics calculations (Spherical Earth):\n\n" ++ + " The following calculations assume a spherical earth model with a radius\n" ++ + " derived from the WGS84 ellipsoid: " ++ + show r84 ++ + "\n" ++ + "\n position track dur position of track after duration\n" ++ + " cpa track1 track2 closest point of approach between two tracks\n" ++ + " intercept track pos minimum speed of interceptor at pos to intercept target\n" ++ + " interceptBySpeed track pos spd time needed by interceptor at pos and travelling at spd to intercept target\n" ++ + " interceptByTime track pos dur speed needed by interceptor at pos to intercept target after duration\n" ++ + "\n Constructors and conversions:\n\n" ++ + " ecef len len len earth-centred earth-fixed position from x, y, z lengths\n" ++ + " ecef metres metres metres earth-centred earth-fixed position from x, y, z metres\n" ++ + " toEcef pos (earth) geographic position to ECEF position using earth model\n" ++ + " WGS84 ellipsoid is used if model is omitted\n" ++ + " fromEcef ecef (earth) ECEF position to geographic position using earth model\n" ++ + " WGS84 ellipsoid is used if model is omitted\n" ++ + " frameB ang ang ang body frame (vehicle) from yaw, pitch and roll angles\n" ++ + " frameL ang local frame from wander azimuth angle\n" ++ + " frameN north, east, down frame\n" ++ + " delta len len len delta from lengths\n" ++ + " delta metres metres metres delta from metres\n" ++ + " ned len len len north, east, down from lengths\n" ++ + " ned metres metres metres north, east, down from metres\n" ++ + " geo latlong surface geographic position from latlong\n" ++ + " geo latlong height geographic position from latlong and height\n" ++ + " geo lat long height geographic position from decimal latitude, longitude and height\n" ++ + " geo lat long metres geographic position from decimal latitude, longitude and metres\n" ++ + " toNVector pos n-vector corresponding to pos\n" ++ + " greatCircle pos1 pos2 great circle passing by pos1 and pos2\n" ++ + " greatCircle pos ang great circle passing by pos and heading on bearing ang\n" ++ + " greatCircle track great circle from track\n" ++ + " track pos ang spd track at pos, heading on bearing ang and travelling at speed spd\n" ++ + "\n Supported lat/long formats:\n\n" ++ + " DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - 553621N0130209E\n" ++ + " d°m's\"[N|S],d°m's\"[E|W] - 55°36'21\"N,13°2'9\"E\n" ++ + " ^ zeroes can be ommitted and separtors can also be d, m, s\n" ++ + " decimal°[N|S],decimal°[E|W] - 51.885°N,13,1°E\n" ++ + "\n Supported angle formats:\n\n" ++ + " d°m's - 55°36'21.154\n" ++ + " decimal° - 51.885°\n" ++ + "\n Supported length formats: {l}m, {l}km, {l}nm, {l}ft\n" ++ + "\n Supported speed formats: {s}m/s, {s}km/h, {s}mph, {s}kt, {s}ft/s\n" ++ + "\n Supported duration formats: (-)nHnMn.nS\n" ++ + "\n Supported earth models:\n\n" ++ + " ellipsoidal: wgs84, grs80, wgs72\n" ++ + " spherical : s84, s80, s72\n" ++ + "\n\n Every evaluated result can be saved by prefixing the expression with \"{var} = \"\n" ++ + " Saved results can subsequently be used when calling a function\n" ++ + "\n Examples:\n\n" ++ + " jord> a = antipode 54N028E\n" ++ + " jord> antipode a\n" ++ + " jord> f = frameB 10d 20d 30d\n" ++ + " jord> d = delta 3000 2000 100\n" ++ + " jord> p0 = geo 49.66618 3.45063 0\n" ++ " jord> target p0 f d wgs84\n" + +save :: Result -> String -> Vault -> Vault +save (Right v) k vault = insert k v vault +save _ _ vault = vault + +showR :: Result -> Either String String +showR (Left err) = Left err +showR (Right v) = Right (show v) + +showVar :: String -> Value -> String +showVar n v = n ++ "=" ++ show v
jord.cabal view
@@ -1,90 +1,96 @@--- This file has been generated from package.yaml by hpack version 0.28.2.------ see: https://github.com/sol/hpack------ hash: 57d929cbf50121806f320b1f6fd26ae04974ba272077bcab7af116e4adbb31ff--name: jord-version: 0.3.1.0-synopsis: Geographical Position Calculations-description: Please see the README on GitHub at <https://github.com/ofmooseandmen/jord#readme>-category: Geography-stability: experimental-homepage: https://github.com/ofmooseandmen/jord-bug-reports: https://github.com/ofmooseandmen/jord/issues-author: Cedric Liegeois-maintainer: Cedric Liegeois <ofmooseandmen@yahoo.com>-copyright: 2018 Cedric Liegeois-license: BSD3-license-file: LICENSE-build-type: Simple-cabal-version: >= 1.10-extra-source-files:- ChangeLog.md- README.md--source-repository head- type: git- location: https://github.com/ofmooseandmen/jord--library- exposed-modules:- Data.Geo.Jord- Data.Geo.Jord.Angle- Data.Geo.Jord.AngularPosition- Data.Geo.Jord.Earth- Data.Geo.Jord.EcefPosition- Data.Geo.Jord.Frames- Data.Geo.Jord.Geodetics- Data.Geo.Jord.LatLong- Data.Geo.Jord.Length- Data.Geo.Jord.NVector- Data.Geo.Jord.Quantity- Data.Geo.Jord.Rotation- Data.Geo.Jord.Transform- Data.Geo.Jord.Vector3d- other-modules:- Data.Geo.Jord.Parse- hs-source-dirs:- src- ghc-options: -Wall- build-depends:- base >=4.9 && <5- default-language: Haskell2010--executable jord-exe- main-is: Main.hs- other-modules:- Eval- Paths_jord- hs-source-dirs:- app- ghc-options: -Wall- build-depends:- base >=4.9 && <5- , haskeline >=0.7 && <0.8- , jord- default-language: Haskell2010--test-suite jord-test- type: exitcode-stdio-1.0- main-is: Spec.hs- other-modules:- Data.Geo.Jord.AngleSpec- Data.Geo.Jord.EarthSpec- Data.Geo.Jord.FramesSpec- Data.Geo.Jord.GeodeticsSpec- Data.Geo.Jord.LatLongSpec- Data.Geo.Jord.LengthSpec- Data.Geo.Jord.RotationSpec- Data.Geo.Jord.TransformSpec- Paths_jord- hs-source-dirs:- test- ghc-options: -Wall- build-depends:- HUnit ==1.6.*- , base >=4.9 && <5- , hspec ==2.*- , jord- default-language: Haskell2010+-- This file has been generated from package.yaml by hpack version 0.28.2. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 080acd4c3e6004551525d8dfe07551b2c77ba7644db59ffa50b4b70d53061a64 + +name: jord +version: 0.4.0.0 +synopsis: Geographical Position Calculations +description: Please see the README on GitHub at <https://github.com/ofmooseandmen/jord#readme> +category: Geography +stability: experimental +homepage: https://github.com/ofmooseandmen/jord +bug-reports: https://github.com/ofmooseandmen/jord/issues +author: Cedric Liegeois +maintainer: Cedric Liegeois <ofmooseandmen@yahoo.com> +copyright: 2018 Cedric Liegeois +license: BSD3 +license-file: LICENSE +build-type: Simple +cabal-version: >= 1.10 +extra-source-files: + ChangeLog.md + README.md + +source-repository head + type: git + location: https://github.com/ofmooseandmen/jord + +library + exposed-modules: + Data.Geo.Jord + Data.Geo.Jord.Angle + Data.Geo.Jord.AngularPosition + Data.Geo.Jord.Duration + Data.Geo.Jord.Earth + Data.Geo.Jord.EcefPosition + Data.Geo.Jord.Frames + Data.Geo.Jord.Geodetics + Data.Geo.Jord.Kinematics + Data.Geo.Jord.LatLong + Data.Geo.Jord.Length + Data.Geo.Jord.NVector + Data.Geo.Jord.Quantity + Data.Geo.Jord.Rotation + Data.Geo.Jord.Speed + Data.Geo.Jord.Transformation + Data.Geo.Jord.Vector3d + other-modules: + Data.Geo.Jord.Parse + hs-source-dirs: + src + ghc-options: -Wall + build-depends: + base >=4.9 && <5 + default-language: Haskell2010 + +executable jord-exe + main-is: Main.hs + other-modules: + Eval + Paths_jord + hs-source-dirs: + app + ghc-options: -Wall + build-depends: + base >=4.9 && <5 + , haskeline >=0.7 && <0.8 + , jord + default-language: Haskell2010 + +test-suite jord-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + Data.Geo.Jord.AngleSpec + Data.Geo.Jord.DurationSpec + Data.Geo.Jord.EarthSpec + Data.Geo.Jord.FramesSpec + Data.Geo.Jord.GeodeticsSpec + Data.Geo.Jord.KinematicsSpec + Data.Geo.Jord.LatLongSpec + Data.Geo.Jord.LengthSpec + Data.Geo.Jord.RotationSpec + Data.Geo.Jord.SpeedSpec + Data.Geo.Jord.TransformationSpec + Paths_jord + hs-source-dirs: + test + ghc-options: -Wall + build-depends: + HUnit ==1.6.* + , base >=4.9 && <5 + , hspec ==2.* + , jord + default-language: Haskell2010
src/Data/Geo/Jord.hs view
@@ -18,34 +18,40 @@ module Data.Geo.Jord ( module Data.Geo.Jord.Angle , module Data.Geo.Jord.AngularPosition + , module Data.Geo.Jord.Duration , module Data.Geo.Jord.Earth , module Data.Geo.Jord.EcefPosition , module Data.Geo.Jord.Frames , module Data.Geo.Jord.Geodetics + , module Data.Geo.Jord.Kinematics , module Data.Geo.Jord.LatLong , module Data.Geo.Jord.Length , module Data.Geo.Jord.NVector , module Data.Geo.Jord.Quantity , module Data.Geo.Jord.Rotation - , module Data.Geo.Jord.Transform + , module Data.Geo.Jord.Speed + , module Data.Geo.Jord.Transformation , module Data.Geo.Jord.Vector3d , jordVersion ) where import Data.Geo.Jord.Angle import Data.Geo.Jord.AngularPosition +import Data.Geo.Jord.Duration import Data.Geo.Jord.Earth import Data.Geo.Jord.EcefPosition import Data.Geo.Jord.Frames import Data.Geo.Jord.Geodetics +import Data.Geo.Jord.Kinematics import Data.Geo.Jord.LatLong import Data.Geo.Jord.Length import Data.Geo.Jord.NVector import Data.Geo.Jord.Quantity import Data.Geo.Jord.Rotation -import Data.Geo.Jord.Transform +import Data.Geo.Jord.Speed +import Data.Geo.Jord.Transformation import Data.Geo.Jord.Vector3d -- | version. jordVersion :: String -jordVersion = "0.3.1.0" +jordVersion = "0.4.0.0"
src/Data/Geo/Jord/Angle.hs view
@@ -8,214 +8,214 @@ -- -- Types and functions for working with angles representing latitudes, longitude and bearings. -- -module Data.Geo.Jord.Angle- (- -- * The 'Angle' type- Angle+module Data.Geo.Jord.Angle + ( + -- * The 'Angle' type + Angle -- * Smart constructors - , decimalDegrees- , dms- , dmsE- , dmsF+ , decimalDegrees + , dms + , dmsE + , dmsF -- * Calculations - , arcLength- , central- , isNegative- , isWithin- , negate'- , normalise+ , arcLength + , central + , isNegative + , isWithin + , negate' + , normalise -- * Trigonometric functions - , asin'- , atan2'- , cos'- , sin'+ , asin' + , atan2' + , cos' + , sin' -- * Accessors - , getDegrees- , getMinutes- , getSeconds- , getMilliseconds- , toDecimalDegrees+ , getDegrees + , getMinutes + , getSeconds + , getMilliseconds + , toDecimalDegrees -- * Read - , angle- , readAngle- , readAngleE- , readAngleF- ) where--import Control.Applicative-import Control.Monad.Fail-import Data.Fixed-import Data.Geo.Jord.Length-import Data.Geo.Jord.Parse-import Data.Geo.Jord.Quantity-import Data.Maybe-import Prelude hiding (fail, length)-import Text.ParserCombinators.ReadP-import Text.Printf-import Text.Read hiding (get, look, pfail)-+ , angle + , readAngle + , readAngleE + , readAngleF + ) where + +import Control.Applicative +import Control.Monad.Fail +import Data.Fixed +import Data.Geo.Jord.Length +import Data.Geo.Jord.Parse +import Data.Geo.Jord.Quantity +import Data.Maybe +import Prelude hiding (fail, length) +import Text.ParserCombinators.ReadP +import Text.Printf +import Text.Read hiding (pfail) + -- | An angle with a resolution of a milliseconds of a degree. -- When used as a latitude/longitude this roughly translate to a precision -- of 30 millimetres at the equator. -newtype Angle = Angle- { milliseconds :: Int- } deriving (Eq)-+newtype Angle = Angle + { milliseconds :: Int + } deriving (Eq) + -- | See 'readAngle'. -instance Read Angle where- readsPrec _ = readP_to_S angle-+instance Read Angle where + readsPrec _ = readP_to_S angle + -- | Angle is shown degrees, minutes, seconds and milliseconds - e.g. 154°25'43.5". -instance Show Angle where- show a =- s ++- show d ++- "°" ++- show (getMinutes a) ++- "'" ++ show (getSeconds a) ++ "." ++ printf "%03d" (getMilliseconds a) ++ "\""- where- d = getDegrees a- s =- if d == 0 && milliseconds a < 0- then "-"- else ""---- | Add/Subtract 'Angle'. -instance Quantity Angle where- add (Angle millis1) (Angle millis2) = Angle (millis1 + millis2)- sub (Angle millis1) (Angle millis2) = Angle (millis1 - millis2)- zero = Angle 0-+instance Show Angle where + show a = + s ++ + show d ++ + "°" ++ + show (getMinutes a) ++ + "'" ++ show (getSeconds a) ++ "." ++ printf "%03d" (getMilliseconds a) ++ "\"" + where + d = getDegrees a + s = + if d == 0 && milliseconds a < 0 + then "-" + else "" + +-- | Add/Subtract 'Angle's. +instance Quantity Angle where + add (Angle millis1) (Angle millis2) = Angle (millis1 + millis2) + sub (Angle millis1) (Angle millis2) = Angle (millis1 - millis2) + zero = Angle 0 + -- | 'Angle' from given decimal degrees. Any 'Double' is accepted: it must be -- validated by the call site when used to represent a latitude or longitude. -decimalDegrees :: Double -> Angle-decimalDegrees dec = Angle (round (dec * 3600000.0))-+decimalDegrees :: Double -> Angle +decimalDegrees dec = Angle (round (dec * 3600000.0)) + -- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. -- 'error's if given minutes, seconds and/or milliseconds are invalid. -- Degrees are not validated and can be any 'Int': they must be validated by the call site -- when used to represent a latitude or longitude. -dms :: Int -> Int -> Int -> Int -> Angle-dms degs mins secs millis =- fromMaybe- (error- ("Invalid minutes=" ++- show mins ++ " or seconds=" ++ show secs ++ " or milliseconds=" ++ show millis))- (dmsF degs mins secs millis)-+dms :: Int -> Int -> Int -> Int -> Angle +dms degs mins secs millis = + fromMaybe + (error + ("Invalid minutes=" ++ + show mins ++ " or seconds=" ++ show secs ++ " or milliseconds=" ++ show millis)) + (dmsF degs mins secs millis) + -- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. -- A 'Left' indicates that given minutes, seconds and/or milliseconds are invalid. -- Degrees are not validated and can be any 'Int': they must be validated by the call site -- when used to represent a latitude or longitude. -dmsE :: Int -> Int -> Int -> Int -> Either String Angle-dmsE degs mins secs millis- | mins < 0 || mins > 59 = Left ("Invalid minutes: " ++ show mins)- | secs < 0 || secs >= 60 = Left ("Invalid seconds: " ++ show secs)- | millis < 0 || millis >= 1000 = Left ("Invalid milliseconds: " ++ show millis)- | otherwise = Right (decimalDegrees ms)- where- ms =- signed- (fromIntegral (abs degs) + (fromIntegral mins / 60.0 :: Double) +- (fromIntegral secs / 3600.0 :: Double) +- (fromIntegral millis / 3600000.0 :: Double))- (signum degs)-+dmsE :: Int -> Int -> Int -> Int -> Either String Angle +dmsE degs mins secs millis + | mins < 0 || mins > 59 = Left ("Invalid minutes: " ++ show mins) + | secs < 0 || secs >= 60 = Left ("Invalid seconds: " ++ show secs) + | millis < 0 || millis >= 1000 = Left ("Invalid milliseconds: " ++ show millis) + | otherwise = Right (decimalDegrees ms) + where + ms = + signed + (fromIntegral (abs degs) + (fromIntegral mins / 60.0 :: Double) + + (fromIntegral secs / 3600.0 :: Double) + + (fromIntegral millis / 3600000.0 :: Double)) + (signum degs) + -- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. -- 'fail's if given minutes, seconds and/or milliseconds are invalid. -- Degrees are not validated and can be any 'Int': they must be validated by the call site -- when used to represent a latitude or longitude. -dmsF :: (MonadFail m) => Int -> Int -> Int -> Int -> m Angle-dmsF degs mins secs millis =- case e of- Left err -> fail err- Right a -> return a- where- e = dmsE degs mins secs millis-+dmsF :: (MonadFail m) => Int -> Int -> Int -> Int -> m Angle +dmsF degs mins secs millis = + case e of + Left err -> fail err + Right a -> return a + where + e = dmsE degs mins secs millis + -- | @arcLength a r@ computes the 'Length' of the arc that subtends the angle @a@ for radius @r@. -arcLength :: Angle -> Length -> Length-arcLength a r = metres (toMetres r * toRadians a)-+arcLength :: Angle -> Length -> Length +arcLength a r = metres (toMetres r * toRadians a) + -- | @central l r@ computes the central 'Angle' from the arc length @l@ and radius @r@. -central :: Length -> Length -> Angle-central s r = fromRadians (toMetres s / toMetres r)-+central :: Length -> Length -> Angle +central s r = fromRadians (toMetres s / toMetres r) + -- | Returns the given 'Angle' negated. -negate' :: Angle -> Angle-negate' (Angle millis) = Angle (-millis)-+negate' :: Angle -> Angle +negate' (Angle millis) = Angle (-millis) + -- | @normalise a n@ normalises @a@ to [0, @n@]. -normalise :: Angle -> Angle -> Angle-normalise a n = decimalDegrees dec- where- dec = mod' (toDecimalDegrees a + toDecimalDegrees n) 360.0-+normalise :: Angle -> Angle -> Angle +normalise a n = decimalDegrees dec + where + dec = mod' (toDecimalDegrees a + toDecimalDegrees n) 360.0 + -- | Is given 'Angle' < 0? -isNegative :: Angle -> Bool-isNegative (Angle millis) = millis < 0-+isNegative :: Angle -> Bool +isNegative (Angle millis) = millis < 0 + -- | Is given 'Angle' within range [@low@..@high@] inclusive? -isWithin :: Angle -> Angle -> Angle -> Bool-isWithin (Angle millis) (Angle low) (Angle high) = millis >= low && millis <= high-+isWithin :: Angle -> Angle -> Angle -> Bool +isWithin (Angle millis) (Angle low) (Angle high) = millis >= low && millis <= high + -- | @atan2' y x@ computes the 'Angle' (from the positive x-axis) of the vector from the origin to the point (x,y). -atan2' :: Double -> Double -> Angle-atan2' y x = fromRadians (atan2 y x)-+atan2' :: Double -> Double -> Angle +atan2' y x = fromRadians (atan2 y x) + -- | @asin' a@ computes arc sine of @a@. -asin' :: Double -> Angle-asin' a = fromRadians (asin a)-+asin' :: Double -> Angle +asin' a = fromRadians (asin a) + -- | @cos' a@ returns the cosinus of @a@. -cos' :: Angle -> Double-cos' a = cos (toRadians a)-+cos' :: Angle -> Double +cos' a = cos (toRadians a) + -- | @sin' a@ returns the sinus of @a@. -sin' :: Angle -> Double-sin' a = sin (toRadians a)-+sin' :: Angle -> Double +sin' a = sin (toRadians a) + -- | radians to degrees. -fromRadians :: Double -> Angle-fromRadians r = decimalDegrees (r / pi * 180.0)-+fromRadians :: Double -> Angle +fromRadians r = decimalDegrees (r / pi * 180.0) + -- | degrees to radians. -toRadians :: Angle -> Double-toRadians a = toDecimalDegrees a * pi / 180.0-+toRadians :: Angle -> Double +toRadians a = toDecimalDegrees a * pi / 180.0 + -- | Converts the given 'Angle' to decimal degrees. -toDecimalDegrees :: Angle -> Double-toDecimalDegrees (Angle millis) = fromIntegral millis / 3600000.0-+toDecimalDegrees :: Angle -> Double +toDecimalDegrees (Angle millis) = fromIntegral millis / 3600000.0 + -- | @getDegrees a@ returns the degree component of @a@. -getDegrees :: Angle -> Int-getDegrees a = signed (field a 3600000.0 360.0) (signum (milliseconds a))-+getDegrees :: Angle -> Int +getDegrees a = signed (field a 3600000.0 360.0) (signum (milliseconds a)) + -- | @getMinutes a@ returns the minute component of @a@. -getMinutes :: Angle -> Int-getMinutes a = field a 60000.0 60.0-+getMinutes :: Angle -> Int +getMinutes a = field a 60000.0 60.0 + -- | @getSeconds a@ returns the second component of @a@. -getSeconds :: Angle -> Int-getSeconds a = field a 1000.0 60.0-+getSeconds :: Angle -> Int +getSeconds a = field a 1000.0 60.0 + -- | @getMilliseconds a@ returns the milliseconds component of @a@. -getMilliseconds :: Angle -> Int-getMilliseconds (Angle millis) = mod (abs millis) 1000--field :: Angle -> Double -> Double -> Int-field (Angle millis) divisor modulo =- truncate (mod' (fromIntegral (abs millis) / divisor) modulo) :: Int--signed :: (Num a, Num b, Ord b) => a -> b -> a-signed n s- | s < 0 = -n- | otherwise = n-+getMilliseconds :: Angle -> Int +getMilliseconds (Angle millis) = mod (abs millis) 1000 + +field :: Angle -> Double -> Double -> Int +field (Angle millis) divisor modulo = + truncate (mod' (fromIntegral (abs millis) / divisor) modulo) :: Int + +signed :: (Num a, Num b, Ord b) => a -> b -> a +signed n s + | s < 0 = -n + | otherwise = n + -- | Parses and returns an 'Angle'. -angle :: ReadP Angle-angle = degsMinsSecs <|> decimal-+angle :: ReadP Angle +angle = degsMinsSecs <|> decimal + -- | Obtains a 'Angle' from the given string formatted as either: -- -- * d°m′s.ms″ - e.g. 55°36'21.3", where minutes, seconds and milliseconds are optional. @@ -232,70 +232,70 @@ -- -- This simply calls @read s :: Angle@ so 'error' should be handled at the call site. -- -readAngle :: String -> Angle-readAngle s = read s :: Angle-+readAngle :: String -> Angle +readAngle s = read s :: Angle + -- | Same as 'readAngle' but returns an 'Either'. -readAngleE :: String -> Either String Angle-readAngleE s =- case readMaybe s of- Nothing -> Left ("couldn't read angle " ++ s)- Just a -> Right a-+readAngleE :: String -> Either String Angle +readAngleE s = + case readMaybe s of + Nothing -> Left ("couldn't read angle " ++ s) + Just a -> Right a + -- | Same as 'readAngle' but returns a 'MonadFail'. -readAngleF :: (MonadFail m) => String -> m Angle-readAngleF s =- let p = readAngleE s- in case p of- Left e -> fail e- Right l -> return l-+readAngleF :: (MonadFail m) => String -> m Angle +readAngleF s = + let p = readAngleE s + in case p of + Left e -> fail e + Right l -> return l + -- | Parses DMS.MS and returns an 'Angle'. -degsMinsSecs :: ReadP Angle-degsMinsSecs = do- d' <- fmap fromIntegral integer- degSymbol- (m', s', ms') <- option (0, 0, 0) (minsSecs <|> minsOnly)- dmsF d' m' s' ms'-+degsMinsSecs :: ReadP Angle +degsMinsSecs = do + d' <- fmap fromIntegral integer + degSymbol + (m', s', ms') <- option (0, 0, 0) (minsSecs <|> minsOnly) + dmsF d' m' s' ms' + -- | Parses minutes, seconds with optionally milliseconds. -minsSecs :: ReadP (Int, Int, Int)-minsSecs = do- m' <- natural- minSymbol- s' <- natural- ms' <- option 0 (char '.' >> natural)- secSymbol- return (m', s', ms')-+minsSecs :: ReadP (Int, Int, Int) +minsSecs = do + m' <- natural + minSymbol + s' <- natural + ms' <- option 0 (char '.' >> natural) + secSymbol + return (m', s', ms') + -- | Parses minutes. -minsOnly :: ReadP (Int, Int, Int)-minsOnly = do- m' <- natural- minSymbol- return (m', 0, 0)-+minsOnly :: ReadP (Int, Int, Int) +minsOnly = do + m' <- natural + minSymbol + return (m', 0, 0) + -- | Parses decimal degrees. -decimal :: ReadP Angle-decimal = do- d <- double- degSymbol- return (decimalDegrees d)-+decimal :: ReadP Angle +decimal = do + d <- double + degSymbol + return (decimalDegrees d) + -- | skips degree symbol. -degSymbol :: ReadP ()-degSymbol = do- _ <- char '°' <|> char 'd'- return ()-+degSymbol :: ReadP () +degSymbol = do + _ <- char '°' <|> char 'd' + return () + -- | skips minute symbol. -minSymbol :: ReadP ()-minSymbol = do- _ <- char '\'' <|> char '′' <|> char 'm'- return ()-+minSymbol :: ReadP () +minSymbol = do + _ <- char '\'' <|> char '′' <|> char 'm' + return () + -- | skips second symbol. -secSymbol :: ReadP ()-secSymbol = do- _ <- string "\"" <|> string "''" <|> string "″" <|> string "s"- return ()+secSymbol :: ReadP () +secSymbol = do + _ <- string "\"" <|> string "''" <|> string "″" <|> string "s" + return ()
src/Data/Geo/Jord/AngularPosition.hs view
@@ -10,50 +10,50 @@ -- -- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> -- -module Data.Geo.Jord.AngularPosition- ( AngularPosition(..)- , latLongHeight- , decimalLatLongHeight- , decimalLatLongHeightE- , decimalLatLongHeightF- , nvectorHeight- ) where--import Control.Monad.Fail-import Data.Geo.Jord.LatLong-import Data.Geo.Jord.Length-import Data.Geo.Jord.NVector-+module Data.Geo.Jord.AngularPosition + ( AngularPosition(..) + , latLongHeight + , decimalLatLongHeight + , decimalLatLongHeightE + , decimalLatLongHeightF + , nvectorHeight + ) where + +import Control.Monad.Fail +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector + -- | An earth position defined by an horizontal position and height. -- -- horizontal position can be either a 'LatLong' or a 'NVector'. -data AngularPosition a = AngularPosition- { pos :: a- , height :: Length- } deriving (Eq, Show)-+data AngularPosition a = AngularPosition + { pos :: a + , height :: Length + } deriving (Eq, Show) + -- | 'AngularPosition' from a 'LatLong' and height. -latLongHeight :: LatLong -> Length -> AngularPosition LatLong-latLongHeight = AngularPosition-+latLongHeight :: LatLong -> Length -> AngularPosition LatLong +latLongHeight = AngularPosition + -- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. -- 'error's if given latitude is outisde [-90..90]° and/or -- given longitude is outisde [-180..180]°. -decimalLatLongHeight :: Double -> Double -> Length -> AngularPosition LatLong-decimalLatLongHeight lat lon = latLongHeight (decimalLatLong lat lon)-+decimalLatLongHeight :: Double -> Double -> Length -> AngularPosition LatLong +decimalLatLongHeight lat lon = latLongHeight (decimalLatLong lat lon) + -- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. -- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or -- given longitude is outisde [-180..180]°. -decimalLatLongHeightE :: Double -> Double -> Length -> Either String (AngularPosition LatLong)-decimalLatLongHeightE lat lon h = fmap (`latLongHeight` h) (decimalLatLongE lat lon)-+decimalLatLongHeightE :: Double -> Double -> Length -> Either String (AngularPosition LatLong) +decimalLatLongHeightE lat lon h = fmap (`latLongHeight` h) (decimalLatLongE lat lon) + -- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. -- 'fail's if given latitude is outisde [-90..90]° and/or -- given longitude is outisde [-180..180]°. -decimalLatLongHeightF :: (MonadFail m) => Double -> Double -> Length -> m (AngularPosition LatLong)-decimalLatLongHeightF lat lon h = fmap (`latLongHeight` h) (decimalLatLongF lat lon)-+decimalLatLongHeightF :: (MonadFail m) => Double -> Double -> Length -> m (AngularPosition LatLong) +decimalLatLongHeightF lat lon h = fmap (`latLongHeight` h) (decimalLatLongF lat lon) + -- | 'AngularPosition' from a 'NVector' and height. -nvectorHeight :: NVector -> Length -> AngularPosition NVector-nvectorHeight = AngularPosition+nvectorHeight :: NVector -> Length -> AngularPosition NVector +nvectorHeight = AngularPosition
+ src/Data/Geo/Jord/Duration.hs view
@@ -0,0 +1,144 @@+-- | +-- Module: Data.Geo.Jord.Duration +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with (signed) durations. +-- +module Data.Geo.Jord.Duration + ( + -- * The 'Duration' type + Duration + , toMilliseconds + -- * Smart constructors + , milliseconds + , hours + , minutes + , seconds + , hms + -- * Accessors + , toHours + , toMinutes + , toSeconds + -- * Read + , readDuration + , readDurationE + , readDurationF + ) where + +import Control.Monad.Fail +import Data.Geo.Jord.Parse +import Data.Geo.Jord.Quantity +import Prelude hiding (fail) +import Text.ParserCombinators.ReadP +import Text.Printf +import Text.Read hiding (pfail) + +-- | A durartion with a resolution of 1 millisecond. +newtype Duration = Duration + { toMilliseconds :: Int -- ^ the number of milliseconds in duration. + } deriving (Eq) + +-- | See 'readDuration'. +instance Read Duration where + readsPrec _ = readP_to_S duration + +-- | show Duration as @(-)nHnMn.nS@. +instance Show Duration where + show d@(Duration millis) = + show h ++ "H" ++ show m ++ "M" ++ show s ++ "." ++ printf "%03d" ms ++ "S" + where + h = truncate (toHours d) :: Int + m = truncate (fromIntegral (millis `mod` 3600000) / 60000.0 :: Double) :: Int + s = truncate (fromIntegral (millis `mod` 60000) / 1000.0 :: Double) :: Int + ms = mod (abs millis) 1000 + +-- | Add/Subtract Durations. +instance Quantity Duration where + add a b = Duration (toMilliseconds a + toMilliseconds b) + sub a b = Duration (toMilliseconds a - toMilliseconds b) + zero = Duration 0 + +-- | 'Duration' from hours minutes and decimal seconds. +hms :: Int -> Int -> Double -> Duration +hms h m s = milliseconds (fromIntegral h * 3600000 + fromIntegral m * 60000 + s * 1000) + +-- | 'Duration' from given amount of hours. +hours :: Double -> Duration +hours h = milliseconds (h * 3600000) + +-- | 'Duration' from given amount of minutes. +minutes :: Double -> Duration +minutes m = milliseconds (m * 60000) + +-- | 'Duration' from given amount of seconds. +seconds :: Double -> Duration +seconds s = milliseconds (s * 1000) + +-- | 'Duration' from given amount of milliseconds. +milliseconds :: Double -> Duration +milliseconds ms = Duration (round ms) + +-- | @toHours d@ gets the number of hours in @d@. +toHours :: Duration -> Double +toHours (Duration ms) = fromIntegral ms / 3600000.0 :: Double + +-- | @toMinutes d@ gets the number of minutes in @d@. +toMinutes :: Duration -> Double +toMinutes (Duration ms) = fromIntegral ms / 60000.0 :: Double + +-- | @toSeconds d@ gets the number of seconds in @d@. +toSeconds :: Duration -> Double +toSeconds (Duration ms) = fromIntegral ms / 1000.0 :: Double + +-- | Obtains a 'Duration' from the given string formatted @(-)nHnMn.nS@. +-- +-- This simply calls @read s :: Duration@ so 'error' should be handled at the call site. +-- +readDuration :: String -> Duration +readDuration s = read s :: Duration + +-- | Same as 'readDuration' but returns a 'Either'. +readDurationE :: String -> Either String Duration +readDurationE s = + case readMaybe s of + Nothing -> Left ("couldn't read duration " ++ s) + Just l -> Right l + +-- | Same as 'readDuration' but returns a 'MonadFail'. +readDurationF :: (MonadFail m) => String -> m Duration +readDurationF s = + let p = readEither s + in case p of + Left e -> fail e + Right l -> return l + +-- | Parses and returns an 'Duration'. +duration :: ReadP Duration +duration = do + h <- option 0 hoursP + m <- option 0 minutesP + s <- option 0.0 secondsP + return (milliseconds (h * 3600000.0 + m * 60000.0 + s * 1000.0)) + +hoursP :: ReadP Double +hoursP = do + h <- integer + _ <- char 'H' + return (fromIntegral h :: Double) + +minutesP :: ReadP Double +minutesP = do + m <- integer + _ <- char 'M' + return (fromIntegral m :: Double) + +secondsP :: ReadP Double +secondsP = do + s <- integer + ms <- option 0 (char '.' >> natural) + _ <- char 'S' + return (fromIntegral s + fromIntegral ms / 10.0)
src/Data/Geo/Jord/Earth.hs view
@@ -8,110 +8,110 @@ -- -- Ellipsoidal and derived spherical earth models. -- -module Data.Geo.Jord.Earth- ( Earth(..)- , Ellipsoid(..)- , eccentricity- , meanRadius- , polarRadius- , spherical- -- * Reference ellipsoids.- , wgs84- , grs80- , wgs72- -- * Spherical models dervived from reference ellipsoids.- , s84- , s80- , s72- , r84- , r80- , r72- ) where--import Data.Geo.Jord.Length---- | Earth model: ellipsoidal or spherical.-data Earth- = Ellipsoidal Ellipsoid- | Spherical Length- deriving (Eq, Show)-+module Data.Geo.Jord.Earth + ( Earth(..) + , Ellipsoid(..) + , eccentricity + , meanRadius + , polarRadius + , spherical + -- * Reference ellipsoids. + , wgs84 + , grs80 + , wgs72 + -- * Spherical models dervived from reference ellipsoids. + , s84 + , s80 + , s72 + , r84 + , r80 + , r72 + ) where + +import Data.Geo.Jord.Length + +-- | Earth model: ellipsoidal or spherical. +data Earth + = Ellipsoidal Ellipsoid + | Spherical Length + deriving (Eq, Show) + -- | Primary ellipsoid parameters. -data Ellipsoid = Ellipsoid+data Ellipsoid = Ellipsoid { equatorialRadius :: Length -- ^ equatorial radius or semi-major axis (a). , inverseFlattening :: Double -- ^ inverse flattening. - } deriving (Eq, Show)-+ } deriving (Eq, Show) + -- | Computes the eccentricity of the given 'Earth' model. -eccentricity :: Earth -> Double-eccentricity (Ellipsoidal e) = sqrt (1.0 - (b * b) / (a * a))- where- a = semiMajorAxis e- b = semiMinorAxis a (inverseFlattening e)-eccentricity (Spherical _) = 0-+eccentricity :: Earth -> Double +eccentricity (Ellipsoidal e) = sqrt (1.0 - (b * b) / (a * a)) + where + a = semiMajorAxis e + b = semiMinorAxis a (inverseFlattening e) +eccentricity (Spherical _) = 0 + -- | Computes the mean radius of the given 'Earth' model. -- -- This radius can be used for geodetic calculations assuming a spherical earth model. -meanRadius :: Earth -> Length-meanRadius (Ellipsoidal e) = metres ((2.0 * a + b) / 3.0)- where- a = semiMajorAxis e- b = semiMinorAxis a (inverseFlattening e)-meanRadius (Spherical r) = r-+meanRadius :: Earth -> Length +meanRadius (Ellipsoidal e) = metres ((2.0 * a + b) / 3.0) + where + a = semiMajorAxis e + b = semiMinorAxis a (inverseFlattening e) +meanRadius (Spherical r) = r + -- | Computes the polar radius or semi-minor axis (b) of the given 'Earth' model. -polarRadius :: Earth -> Length-polarRadius (Ellipsoidal e) = metres (semiMinorAxis a f)- where- a = semiMajorAxis e- f = inverseFlattening e-polarRadius (Spherical r) = r---- | Spherical model derived from given model.-spherical :: Earth -> Earth-spherical e = Spherical (meanRadius e)-+polarRadius :: Earth -> Length +polarRadius (Ellipsoidal e) = metres (semiMinorAxis a f) + where + a = semiMajorAxis e + f = inverseFlattening e +polarRadius (Spherical r) = r + +-- | Spherical model derived from given model. +spherical :: Earth -> Earth +spherical e = Spherical (meanRadius e) + -- | World Geodetic System WGS84 ellipsoid. -wgs84 :: Earth-wgs84 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257223563))-+wgs84 :: Earth +wgs84 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257223563)) + -- | Geodetic Reference System 1980 ellipsoid. -grs80 :: Earth-grs80 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257222101))-+grs80 :: Earth +grs80 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257222101)) + -- | World Geodetic System WGS72 ellipsoid. -wgs72 :: Earth-wgs72 = Ellipsoidal (Ellipsoid (metres 6378135.0) (1.0 / 298.26))---- | Spherical earth model derived from 'wgs84'.-s84 :: Earth-s84 = spherical wgs84-+wgs72 :: Earth +wgs72 = Ellipsoidal (Ellipsoid (metres 6378135.0) (1.0 / 298.26)) + +-- | Spherical earth model derived from 'wgs84'. +s84 :: Earth +s84 = spherical wgs84 + -- | Spherical earth model derived from 'grs80'. -s80 :: Earth-s80 = spherical grs80-+s80 :: Earth +s80 = spherical grs80 + -- | Spherical earth model derived from 'wgs72'. -s72 :: Earth-s72 = spherical wgs72-+s72 :: Earth +s72 = spherical wgs72 + -- | Mean earth radius derived from the 'wgs84' ellipsoid. -r84 :: Length-r84 = meanRadius s84-+r84 :: Length +r84 = meanRadius s84 + -- | Mean earth radius derived from the 'grs80' ellipsoid. -r80 :: Length-r80 = meanRadius s80-+r80 :: Length +r80 = meanRadius s80 + -- | Mean earth radius derived from the 'wgs72' ellipsoid. -r72 :: Length-r72 = meanRadius s72-+r72 :: Length +r72 = meanRadius s72 + -- | semi-major axis (a) in metres. -semiMajorAxis :: Ellipsoid -> Double-semiMajorAxis = toMetres . equatorialRadius-+semiMajorAxis :: Ellipsoid -> Double +semiMajorAxis = toMetres . equatorialRadius + -- | Computes the polar semi-minor axis (b) from @a@ anf @f@. -semiMinorAxis :: Double -> Double -> Double-semiMinorAxis a f = a * (1.0 - f)+semiMinorAxis :: Double -> Double -> Double +semiMinorAxis a f = a * (1.0 - f)
src/Data/Geo/Jord/EcefPosition.hs view
@@ -10,50 +10,50 @@ -- -- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> -- -module Data.Geo.Jord.EcefPosition- ( EcefPosition- , ecef- , ecefMetres- , ex- , ey- , ez- ) where--import Data.Geo.Jord.Length-import Data.Geo.Jord.Vector3d-+module Data.Geo.Jord.EcefPosition + ( EcefPosition + , ecef + , ecefMetres + , ex + , ey + , ez + ) where + +import Data.Geo.Jord.Length +import Data.Geo.Jord.Vector3d + -- | An earth position expressed in the Earth Centred, Earth Fixed (ECEF) coordinates system. -- -- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. -- -- Note: on a spherical model earth, an n-vector is equivalent to a normalised version of an (ECEF) cartesian coordinate. -newtype EcefPosition =- EcefPosition Vector3d- deriving (Eq, Show)--instance IsVector3d EcefPosition where- vec (EcefPosition v) = v-+newtype EcefPosition = + EcefPosition Vector3d + deriving (Eq, Show) + +instance IsVector3d EcefPosition where + vec (EcefPosition v) = v + -- | 'EcefPosition' from given x, y and z length. -- -- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. -ecef :: Length -> Length -> Length -> EcefPosition-ecef x y z = EcefPosition (Vector3d (toMetres x) (toMetres y) (toMetres z))---- | 'EcefPosition' from given x, y and z length in _metres_. +ecef :: Length -> Length -> Length -> EcefPosition +ecef x y z = EcefPosition (Vector3d (toMetres x) (toMetres y) (toMetres z)) + +-- | 'EcefPosition' from given x, y and z length in __metres__. -- -- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. -ecefMetres :: Double -> Double -> Double -> EcefPosition-ecefMetres x y z = ecef (metres x) (metres y) (metres z)---- | x coordinate of the given '$tc'EcefPosition'. -ex :: EcefPosition -> Length-ex (EcefPosition v) = metres (vx v)---- | y coordinate of the given '$tc'EcefPosition'. -ey :: EcefPosition -> Length-ey (EcefPosition v) = metres (vy v)---- | z coordinate of the given '$tc'EcefPosition'. -ez :: EcefPosition -> Length-ez (EcefPosition v) = metres (vz v)+ecefMetres :: Double -> Double -> Double -> EcefPosition +ecefMetres x y z = ecef (metres x) (metres y) (metres z) + +-- | x coordinate of the given 'EcefPosition'. +ex :: EcefPosition -> Length +ex (EcefPosition v) = metres (vx v) + +-- | y coordinate of the given 'EcefPosition'. +ey :: EcefPosition -> Length +ey (EcefPosition v) = metres (vy v) + +-- | z coordinate of the given 'EcefPosition'. +ez :: EcefPosition -> Length +ez (EcefPosition v) = metres (vz v)
src/Data/Geo/Jord/Frames.hs view
@@ -34,6 +34,9 @@ , Delta , delta , deltaMetres + , dx + , dy + , dz -- * Delta in the north, east, down frame , Ned , ned @@ -43,7 +46,7 @@ , down , bearing , elevation - , norm + , slantRange -- * Calculations , deltaBetween , nedBetween @@ -59,7 +62,7 @@ import Data.Geo.Jord.Length import Data.Geo.Jord.NVector import Data.Geo.Jord.Rotation -import Data.Geo.Jord.Transform +import Data.Geo.Jord.Transformation import Data.Geo.Jord.Vector3d -- | class for reference frames. @@ -198,10 +201,22 @@ delta :: Length -> Length -> Length -> Delta delta x y z = Delta (Vector3d (toMetres x) (toMetres y) (toMetres z)) --- | 'Delta' from given x, y and z length in _metres_. +-- | 'Delta' from given x, y and z length in __metres__. deltaMetres :: Double -> Double -> Double -> Delta deltaMetres x y z = delta (metres x) (metres y) (metres z) +-- | x component of given 'Delta'. +dx :: Delta -> Length +dx (Delta v) = metres (vx v) + +-- | y component of given 'Delta'. +dy :: Delta -> Length +dy (Delta v) = metres (vy v) + +-- | z component of given 'Delta'. +dz :: Delta -> Length +dz (Delta v) = metres (vz v) + -- | North, east and down delta (thus in frame 'FrameN'). newtype Ned = Ned Vector3d @@ -211,7 +226,7 @@ ned :: Length -> Length -> Length -> Ned ned n e d = Ned (Vector3d (toMetres n) (toMetres e) (toMetres d)) --- | 'Ned' from given north, east and down in _metres_. +-- | 'Ned' from given north, east and down in __metres__. nedMetres :: Double -> Double -> Double -> Ned nedMetres n e d = ned (metres n) (metres e) (metres d) @@ -240,12 +255,20 @@ elevation :: Ned -> Angle elevation (Ned v) = negate' (asin' (vz v / vnorm v)) --- | @norm v@ computes the norm of the NED vector @v@. -norm :: Ned -> Length -norm (Ned v) = metres (vnorm v) +-- | @slantRange v@ computes the distance from origin in the local system of the NED vector @v@. +slantRange :: Ned -> Length +slantRange (Ned v) = metres (vnorm v) -- | @deltaBetween p1 p2 f e@ computes the exact 'Delta' between the two positions @p1@ and @p2@ in frame @f@ -- using earth model @e@. +-- +-- @ +-- let p1 = decimalLatLongHeight 1 2 (metres (-3)) +-- let p2 = decimalLatLongHeight 4 5 (metres (-6)) +-- let w = decimalDegrees 5 -- wander azimuth +-- let d = deltaBetween p1 p2 (frameL w) wgs84 +-- d = deltaMetres 359490.579 302818.523 17404.272 +-- @ deltaBetween :: (ETransform a, Frame c) => a -> a -> (a -> Earth -> c) -> Earth -> Delta deltaBetween p1 p2 f e = deltaMetres (vx d) (vy d) (vz d) where @@ -263,6 +286,16 @@ -- the north, east, and down directions will change (relative to Earth) for different places. -- -- Position @p1@ must be outside the poles for the north and east directions to be defined. +-- +-- @ +-- let p1 = decimalLatLongHeight 1 2 (metres (-3)) +-- let p2 = decimalLatLongHeight 4 5 (metres (-6)) +-- let d1 = nedBetween p1 p2 wgs84 +-- let d2 = deltaBetween p1 p2 frameN wgs84 +-- north d1 = dx d2 +-- east d1 = dy d2 +-- down d1 = dz d2 +-- @ nedBetween :: (ETransform a) => a -> a -> Earth -> Ned nedBetween p1 p2 e = nedMetres (vx d) (vy d) (vz d) where @@ -270,6 +303,15 @@ -- | @target p0 f d e@ computes the target position from position @p0@ and delta @d@ using in frame @f@ -- and using earth model @e@. +-- +-- @ +-- let p0 = decimalLatLongHeight 49.66618 3.45063 zero +-- let y = decimalDegrees 10 -- yaw +-- let r = decimalDegrees 20 -- roll +-- let p = decimalDegrees 30 -- pitch +-- let d = deltaMetres 3000 2000 100 +-- target p0 (frameB y r p) d wgs84 = decimalLatLongHeight 49.6918016 3.4812669 (metres 6.007) +-- @ target :: (ETransform a, Frame c) => a -> (a -> Earth -> c) -> Delta -> Earth -> a target p0 f (Delta d) e = fromEcef (ecefMetres (vx e0 + vx c) (vy e0 + vy c) (vz e0 + vz c)) e where @@ -278,6 +320,11 @@ c = vrotate d rm -- | @targetN p0 d e@ computes the target position from position @p0@ and north, east, down @d@ using earth model @e@. +-- +-- @ +-- let p0 = decimalLatLongHeight 49.66618 3.45063 zero +-- targetN p0 (nedMeters 100 200 300) wgs84 = target p0 frameN (deltaMetres 100 200 300) wgs84 +-- @ targetN :: (ETransform a) => a -> Ned -> Earth -> a targetN p0 (Ned d) = target p0 frameN (Delta d)
src/Data/Geo/Jord/Geodetics.hs view
@@ -1,3 +1,5 @@+{-# LANGUAGE FlexibleInstances #-}+ -- | -- Module: Data.Geo.Jord.Geodetics -- Copyright: (c) 2018 Cedric Liegeois @@ -11,47 +13,43 @@ -- All functions are implemented using the vector-based approached described in -- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation> -- -module Data.Geo.Jord.Geodetics - ( - -- * The 'GreatCircle' type - GreatCircle - -- * Smart constructors - , greatCircle - , greatCircleE - , greatCircleF - , greatCircleBearing +module Data.Geo.Jord.Geodetics+ (+ -- * The 'GreatCircle' type+ GreatCircle+ , IsGreatCircle(..) -- * Calculations - , angularDistance - , antipode - , crossTrackDistance - , crossTrackDistance84 - , destination - , destination84 - , finalBearing - , initialBearing - , interpolate - , intersections - , insideSurface - , mean - , surfaceDistance - , surfaceDistance84 - ) where - -import Control.Monad.Fail -import Data.Fixed -import Data.Geo.Jord.Angle -import Data.Geo.Jord.AngularPosition -import Data.Geo.Jord.Earth (r84) -import Data.Geo.Jord.LatLong -import Data.Geo.Jord.Length -import Data.Geo.Jord.NVector -import Data.Geo.Jord.Quantity -import Data.Geo.Jord.Transform -import Data.Geo.Jord.Vector3d -import Data.List (subsequences) -import Data.Maybe (fromMaybe) -import Prelude hiding (fail) - + , angularDistance+ , antipode+ , crossTrackDistance+ , crossTrackDistance84+ , destination+ , destination84+ , finalBearing+ , initialBearing+ , interpolate+ , intersections+ , insideSurface+ , mean+ , surfaceDistance+ , surfaceDistance84+ ) where++import Control.Monad.Fail+import Data.Fixed+import Data.Geo.Jord.Angle+import Data.Geo.Jord.AngularPosition+import Data.Geo.Jord.Earth (r84)+import Data.Geo.Jord.LatLong+import Data.Geo.Jord.Length+import Data.Geo.Jord.NVector+import Data.Geo.Jord.Quantity+import Data.Geo.Jord.Transformation+import Data.Geo.Jord.Vector3d+import Data.List (subsequences)+import Data.Maybe (fromMaybe)+import Prelude hiding (fail)+ -- | A circle on the __surface__ of the Earth which lies in a plane passing through -- the Earth's centre. Every two distinct and non-antipodal points on the surface -- of the Earth define a Great Circle. @@ -61,74 +59,82 @@ -- -- See 'greatCircle', 'greatCircleE', 'greatCircleF' or 'greatCircleBearing' constructors. -- -data GreatCircle = GreatCircle - { normal :: Vector3d - , dscr :: String - } deriving (Eq) - -instance Show GreatCircle where - show = dscr - --- | 'GreatCircle' passing by both given positions. 'error's if given positions are --- equal or antipodal. -greatCircle :: (NTransform a, Show a) => a -> a -> GreatCircle -greatCircle p1 p2 = - fromMaybe - (error (show p1 ++ " and " ++ show p2 ++ " do not define a unique Great Circle")) - (greatCircleF p1 p2) - --- | 'GreatCircle' passing by both given positions. A 'Left' indicates that given positions are --- equal or antipodal. -greatCircleE :: (NTransform a) => a -> a -> Either String GreatCircle -greatCircleE p1 p2 - | v1 == v2 = Left "Invalid Great Circle: positions are equal" - | (realToFrac (vnorm (vadd v1 v2)) :: Nano) == 0 = - Left "Invalid Great Circle: positions are antipodal" - | otherwise = - Right (GreatCircle (vcross v1 v2) ("passing by " ++ show (ll p1) ++ " & " ++ show (ll p2))) - where - v1 = vector3d p1 - v2 = vector3d p2 - --- | 'GreatCircle' passing by both given positions. 'fail's if given positions are +data GreatCircle = GreatCircle+ { normal :: Vector3d+ , dscr :: String+ } deriving (Eq)++instance Show GreatCircle where+ show = dscr++-- | Class for data from which a 'GreatCircle' can be computed. +class (Show a) =>+ IsGreatCircle a+ where+ greatCircle :: a -> GreatCircle -- ^ 'GreatCircle' from @a@, if 'greateCircleE' returns a 'Left', this function 'error's. + greatCircle a = fromMaybe (error (show a ++ " do not define a Great Circle")) (greatCircleF a)+ greatCircleE :: a -> Either String GreatCircle -- ^ 'GreatCircle' from @a@, A 'Left' indicates an error. + greatCircleF :: (MonadFail m) => a -> m GreatCircle -- ^ 'GreatCircle' from @a@, if 'greateCircleE' returns a 'Left', this function 'fail's. + greatCircleF a =+ case e of+ Left err -> fail err+ Right gc -> return gc+ where+ e = greatCircleE a++-- | 'GreatCircle' passing by both given positions'. A 'Left' indicates that given positions are -- equal or antipodal. -greatCircleF :: (NTransform a, MonadFail m) => a -> a -> m GreatCircle -greatCircleF p1 p2 = - case e of - Left err -> fail err - Right gc -> return gc - where - e = greatCircleE p1 p2 - +-- +-- @ +-- let p1 = decimalLatLongHeight 45.0 (-143.5) (metres 1500) +-- let p2 = decimalLatLongHeight 46.0 14.5 (metres 3000) +-- greatCircle (p1, p2) -- heights are ignored, great circle are always at earth surface. +-- @ +instance (NTransform a, Show a) => IsGreatCircle (a, a) where+ greatCircleE (p1, p2)+ | v1 == v2 = Left "Invalid Great Circle: positions are equal"+ | (realToFrac (vnorm (vadd v1 v2)) :: Nano) == 0 =+ Left "Invalid Great Circle: positions are antipodal"+ | otherwise =+ Right+ (GreatCircle (vcross v1 v2) ("passing by " ++ show (ll p1) ++ " & " ++ show (ll p2)))+ where+ v1 = vector3d p1+ v2 = vector3d p2+ -- | 'GreatCircle' passing by the given position and heading on given bearing. -greatCircleBearing :: (NTransform a) => a -> Angle -> GreatCircle -greatCircleBearing p b = - GreatCircle (vsub n' e') ("passing by " ++ show (ll p) ++ " heading on " ++ show b) - where - v = vector3d p - e = vcross (vec northPole) v -- easting - n = vcross v e -- northing - e' = vscale e (cos' b / vnorm e) - n' = vscale n (sin' b / vnorm n) - +-- +-- @ +-- greatCircle (readLatLong "283321N0290700W", decimalDegrees 33.0) +-- @ +instance (NTransform a, Show a) => IsGreatCircle (a, Angle) where+ greatCircleE (p, b) =+ Right (GreatCircle (vsub n' e') ("passing by " ++ show (ll p) ++ " heading on " ++ show b))+ where+ v = vector3d p+ e = vcross (vec northPole) v -- easting + n = vcross v e -- northing + e' = vscale e (cos' b / vnorm e)+ n' = vscale n (sin' b / vnorm n)+ -- | @angularDistance p1 p2 n@ computes the angle between the horizontal positions @p1@ and @p2@. -- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], -- signed + if @p1@ is clockwise looking along @n@, - in opposite direction. -angularDistance :: (NTransform a) => a -> a -> Maybe a -> Angle -angularDistance p1 p2 n = angularDistance' v1 v2 vn - where - v1 = vector3d p1 - v2 = vector3d p2 - vn = fmap vector3d n - +angularDistance :: (NTransform a) => a -> a -> Maybe a -> Angle+angularDistance p1 p2 n = angularDistance' v1 v2 vn+ where+ v1 = vector3d p1+ v2 = vector3d p2+ vn = fmap vector3d n+ -- | @antipode p@ computes the antipodal horizontal position of @p@: -- the horizontal position on the surface of the Earth which is diametrically opposite to @p@. -antipode :: (NTransform a) => a -> a -antipode p = fromNVector (angular (vscale (vector3d nv) (-1.0)) h) - where - (AngularPosition nv h) = toNVector p - --- | @crossTrackDistance p gc@ computes the signed distance horizontal position @p@ to great circle @gc@. +antipode :: (NTransform a) => a -> a+antipode p = fromNVector (angular (vscale (vector3d nv) (-1.0)) h)+ where+ (AngularPosition nv h) = toNVector p++-- | @crossTrackDistance p gc@ computes the signed distance from horizontal position @p@ to great circle @gc@. -- Returns a negative 'Length' if position if left of great circle, -- positive 'Length' if position if right of great circle; the orientation of the -- great circle is therefore important: @@ -136,97 +142,102 @@ -- @ -- let gc1 = greatCircle (decimalLatLong 51 0) (decimalLatLong 52 1) -- let gc2 = greatCircle (decimalLatLong 52 1) (decimalLatLong 51 0) --- crossTrackDistance p gc1 == (- crossTrackDistance p gc2) +-- crossTrackDistance p gc1 = (- crossTrackDistance p gc2) +-- +-- let p = decimalLatLong 53.2611 (-0.7972) +-- let gc = greatCircleBearing (decimalLatLong 53.3206 (-1.7297)) (decimalDegrees 96.0) +-- crossTrackDistance p gc r84 -- -305.663 metres -- @ -crossTrackDistance :: (NTransform a) => a -> GreatCircle -> Length -> Length -crossTrackDistance p gc = - arcLength (sub (angularDistance' (normal gc) (vector3d p) Nothing) (decimalDegrees 90)) - +crossTrackDistance :: (NTransform a) => a -> GreatCircle -> Length -> Length+crossTrackDistance p gc =+ arcLength (sub (angularDistance' (normal gc) (vector3d p) Nothing) (decimalDegrees 90))+ -- | 'crossTrackDistance' using the mean radius of the WGS84 reference ellipsoid. -crossTrackDistance84 :: (NTransform a) => a -> GreatCircle -> Length -crossTrackDistance84 p gc = crossTrackDistance p gc r84 - +crossTrackDistance84 :: (NTransform a) => a -> GreatCircle -> Length+crossTrackDistance84 p gc = crossTrackDistance p gc r84+ -- | @destination p b d r@ computes the destination position from position @p@ having -- travelled the distance @d@ on the initial bearing (compass angle) @b@ (bearing will normally vary -- before destination is reached) and using the earth radius @r@. -destination :: (NTransform a) => a -> Angle -> Length -> Length -> a -destination p b d r - | toMetres d == 0.0 = p - | otherwise = fromNVector (angular vd h) - where - (AngularPosition nv h) = toNVector p - v = vec nv +-- +-- @ +-- let p0 = ecefToNVector (ecefMetres 3812864.094 (-115142.863) 5121515.161) s84 +-- let p1 = ecefMetres 3826406.4710518294 8900.536398998282 5112694.233184049 +-- let p = destination p0 (decimalDegrees 96.0217) (metres 124800) r84 +-- nvectorToEcef p s84 = p1 +-- @ +destination :: (NTransform a) => a -> Angle -> Length -> Length -> a+destination p b d r+ | toMetres d == 0.0 = p+ | otherwise = fromNVector (angular vd h)+ where+ (AngularPosition nv h) = toNVector p+ v = vec nv ed = vunit (vcross (vec northPole) v) -- east direction vector at v nd = vcross v ed -- north direction vector at v ta = central d r -- central angle de = vadd (vscale nd (cos' b)) (vscale ed (sin' b)) -- vunit vector in the direction of the azimuth - vd = vadd (vscale v (cos' ta)) (vscale de (sin' ta)) - + vd = vadd (vscale v (cos' ta)) (vscale de (sin' ta))+ -- | 'destination' using the mean radius of the WGS84 reference ellipsoid. -destination84 :: (NTransform a) => a -> Angle -> Length -> a -destination84 p b d = destination p b d r84 - +destination84 :: (NTransform a) => a -> Angle -> Length -> a+destination84 p b d = destination p b d r84+ -- | @finalBearing p1 p2@ computes the final bearing arriving at @p2@ from @p1@ in compass angle. -- -- Compass angles are clockwise angles from true north: 0 = north, 90 = east, 180 = south, 270 = west. -- --- The final bearing will differ from the 'initialBearing' by varying degrees according to distance and latitude. +-- The final bearing will differ from the 'initialBearing' by varying degrees according to distance and latitude. -- -- Returns 'Nothing' if both horizontal positions are equals. -finalBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle -finalBearing p1 p2 = fmap (\b -> normalise b (decimalDegrees 180)) (initialBearing p2 p1) - +finalBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle+finalBearing p1 p2 = fmap (\b -> normalise b (decimalDegrees 180)) (initialBearing p2 p1)+ -- | @initialBearing p1 p2@ computes the initial bearing from @p1@ to @p2@ in compass angle. -- -- Compass angles are clockwise angles from true north: 0 = north, 90 = east, 180 = south, 270 = west. -- -- Returns 'Nothing' if both horizontal positions are equals. -initialBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle -initialBearing p1 p2 - | p1 == p2 = Nothing - | otherwise = Just (normalise (angularDistance' gc1 gc2 (Just v1)) (decimalDegrees 360)) - where - v1 = vector3d p1 - v2 = vector3d p2 +initialBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle+initialBearing p1 p2+ | p1 == p2 = Nothing+ | otherwise = Just (normalise (angularDistance' gc1 gc2 (Just v1)) (decimalDegrees 360))+ where+ v1 = vector3d p1+ v2 = vector3d p2 gc1 = vcross v1 v2 -- great circle through p1 & p2 gc2 = vcross v1 (vec northPole) -- great circle through p1 & north pole - + -- | @interpolate p0 p1 f# computes the horizontal position at fraction @f@ between the @p0@ and @p1@. -- -- Special conditions: -- -- @ --- interpolate p0 p1 0.0 == p0 --- interpolate p0 p1 1.0 == p1 +-- interpolate p0 p1 0.0 = p0 +-- interpolate p0 p1 1.0 = p1 -- @ -- --- 'error's if @f < 0 || f > 1.0@ +-- 'error's if @f < 0 || f > 1@ -- -interpolate :: (NTransform a) => a -> a -> Double -> a -interpolate p0 p1 f - | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f) - | f == 0 = p0 - | f == 1 = p1 - | otherwise = fromNVector (angular iv ih) - where - (AngularPosition nv0 h0) = toNVector p0 - (AngularPosition nv1 h1) = toNVector p1 - v0 = vec nv0 - v1 = vec nv1 - iv = vunit (vadd v0 (vscale (vsub v1 v0) f)) - ih = lrph h0 h1 f - --- | Computes the intersections between the two given 'GreatCircle's. --- Two 'GreatCircle's intersect exactly twice unless there are equal (regardless of orientation), --- in which case 'Nothing' is returned. -intersections :: (NTransform a) => GreatCircle -> GreatCircle -> Maybe (a, a) -intersections gc1 gc2 - | (vnorm i :: Double) == 0.0 = Nothing - | otherwise - , let ni = fromNVector (angular (vunit i) zero) = Just (ni, antipode ni) - where - i = vcross (normal gc1) (normal gc2) - +-- @ +-- let p1 = latLongHeight (readLatLong "53°28'46''N 2°14'43''W") (metres 10000) +-- let p2 = latLongHeight (readLatLong "55°36'21''N 13°02'09''E") (metres 20000) +-- interpolate p1 p2 0.5 = decimalLatLongHeight 54.7835574 5.1949856 (metres 15000) +-- @ +interpolate :: (NTransform a) => a -> a -> Double -> a+interpolate p0 p1 f+ | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f)+ | f == 0 = p0+ | f == 1 = p1+ | otherwise = fromNVector (angular iv ih)+ where+ (AngularPosition nv0 h0) = toNVector p0+ (AngularPosition nv1 h1) = toNVector p1+ v0 = vec nv0+ v1 = vec nv1+ iv = vunit (vadd v0 (vscale (vsub v1 v0) f))+ ih = lrph h0 h1 f+ -- | @insideSurface p ps@ determines whether the @p@ is inside the polygon defined by the list of positions @ps@. -- The polygon is closed if needed (i.e. if @head ps /= last ps@). -- @@ -235,22 +246,53 @@ -- -- Always returns 'False' if @ps@ does not at least defines a triangle. -- -insideSurface :: (Eq a, NTransform a) => a -> [a] -> Bool -insideSurface p ps - | null ps = False - | head ps == last ps = insideSurface p (init ps) - | length ps < 3 = False - | otherwise = - let aSum = - foldl - (\a v' -> add a (uncurry angularDistance' v' (Just v))) - (decimalDegrees 0) - (egdes (map (vsub v) vs)) - in abs (toDecimalDegrees aSum) > 180.0 - where - v = vector3d p - vs = fmap vector3d ps - +-- @ +-- let malmo = decimalLatLong 55.6050 13.0038 +-- let ystad = decimalLatLong 55.4295 13.82 +-- let lund = decimalLatLong 55.7047 13.1910 +-- let helsingborg = decimalLatLong 56.0465 12.6945 +-- let kristianstad = decimalLatLong 56.0294 14.1567 +-- let polygon = [malmo, ystad, kristianstad, helsingborg, lund] +-- let hoor = decimalLatLong 55.9295 13.5297 +-- let hassleholm = decimalLatLong 56.1589 13.7668 +-- insideSurface hoor polygon = True +-- insideSurface hassleholm polygon = False +-- @ +insideSurface :: (Eq a, NTransform a) => a -> [a] -> Bool+insideSurface p ps+ | null ps = False+ | head ps == last ps = insideSurface p (init ps)+ | length ps < 3 = False+ | otherwise =+ let aSum =+ foldl+ (\a v' -> add a (uncurry angularDistance' v' (Just v)))+ (decimalDegrees 0)+ (egdes (map (vsub v) vs))+ in abs (toDecimalDegrees aSum) > 180.0+ where+ v = vector3d p+ vs = fmap vector3d ps++-- | Computes the intersections between the two given 'GreatCircle's. +-- Two 'GreatCircle's intersect exactly twice unless there are equal (regardless of orientation), +-- in which case 'Nothing' is returned. +-- +-- @ +-- let gc1 = greatCircleBearing (decimalLatLong 51.885 0.235) (decimalDegrees 108.63) +-- let gc2 = greatCircleBearing (decimalLatLong 49.008 2.549) (decimalDegrees 32.72) +-- let (i1, i2) = fromJust (intersections gc1 gc2) +-- i1 = decimalLatLong 50.9017226 4.4942782 +-- i2 = antipode i1 +-- @ +intersections :: (NTransform a) => GreatCircle -> GreatCircle -> Maybe (a, a)+intersections gc1 gc2+ | (vnorm i :: Double) == 0.0 = Nothing+ | otherwise+ , let ni = fromNVector (angular (vunit i) zero) = Just (ni, antipode ni)+ where+ i = vcross (normal gc1) (normal gc2)+ -- | @mean ps@ computes the mean geographic horitzontal position of @ps@, if it is defined. -- -- The geographic mean is not defined for antipodals position (since they @@ -259,59 +301,58 @@ -- Special conditions: -- -- @ --- mean [] == Nothing --- mean [p] == Just p --- mean [p1, p2, p3] == Just circumcentre --- mean [p1, .., antipode p1] == Nothing +-- mean [] = Nothing +-- mean [p] = Just p +-- mean [p1, p2, p3] = Just circumcentre +-- mean [p1, .., antipode p1] = Nothing -- @ --- -mean :: (NTransform a) => [a] -> Maybe a -mean [] = Nothing -mean [p] = Just p -mean ps = - if null antipodals - then Just (fromNVector (angular (vunit (foldl vadd vzero vs)) zero)) - else Nothing - where - vs = fmap vector3d ps - ts = filter (\l -> length l == 2) (subsequences vs) - antipodals = - filter (\t -> (realToFrac (vnorm (vadd (head t) (last t)) :: Double) :: Nano) == 0) ts - +mean :: (NTransform a) => [a] -> Maybe a+mean [] = Nothing+mean [p] = Just p+mean ps =+ if null antipodals+ then Just (fromNVector (angular (vunit (foldl vadd vzero vs)) zero))+ else Nothing+ where+ vs = fmap vector3d ps+ ts = filter (\l -> length l == 2) (subsequences vs)+ antipodals =+ filter (\t -> (realToFrac (vnorm (vadd (head t) (last t)) :: Double) :: Nano) == 0) ts+ -- | @surfaceDistance p1 p2@ computes the surface distance (length of geodesic) between the positions @p1@ and @p2@. -surfaceDistance :: (NTransform a) => a -> a -> Length -> Length -surfaceDistance p1 p2 = arcLength (angularDistance p1 p2 Nothing) - +surfaceDistance :: (NTransform a) => a -> a -> Length -> Length+surfaceDistance p1 p2 = arcLength (angularDistance p1 p2 Nothing)+ -- | 'surfaceDistance' using the mean radius of the WGS84 reference ellipsoid. -surfaceDistance84 :: (NTransform a) => a -> a -> Length -surfaceDistance84 p1 p2 = surfaceDistance p1 p2 r84 - +surfaceDistance84 :: (NTransform a) => a -> a -> Length+surfaceDistance84 p1 p2 = surfaceDistance p1 p2 r84+ -- | Angle between the two given n-vectors. -- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], -- signed + if @v1@ is clockwise looking along @n@, - in opposite direction. -angularDistance' :: Vector3d -> Vector3d -> Maybe Vector3d -> Angle -angularDistance' v1 v2 n = atan2' sinO cosO - where - sign = maybe 1 (signum . vdot (vcross v1 v2)) n - sinO = sign * vnorm (vcross v1 v2) - cosO = vdot v1 v2 - +angularDistance' :: Vector3d -> Vector3d -> Maybe Vector3d -> Angle+angularDistance' v1 v2 n = atan2' sinO cosO+ where+ sign = maybe 1 (signum . vdot (vcross v1 v2)) n+ sinO = sign * vnorm (vcross v1 v2)+ cosO = vdot v1 v2+ -- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] -egdes :: [Vector3d] -> [(Vector3d, Vector3d)] -egdes ps = zip ps (tail ps ++ [head ps]) - -lrph :: Length -> Length -> Double -> Length -lrph h0 h1 f = metres h - where - h0' = toMetres h0 - h1' = toMetres h1 - h = h0' + (h1' - h0') * f - -vector3d :: (NTransform a) => a -> Vector3d -vector3d = vec . pos . toNVector - -angular :: Vector3d -> Length -> AngularPosition NVector -angular v = nvectorHeight (nvector (vx v) (vy v) (vz v)) - -ll :: (NTransform a) => a -> LatLong -ll = nvectorToLatLong . pos . toNVector +egdes :: [Vector3d] -> [(Vector3d, Vector3d)]+egdes ps = zip ps (tail ps ++ [head ps])++lrph :: Length -> Length -> Double -> Length+lrph h0 h1 f = metres h+ where+ h0' = toMetres h0+ h1' = toMetres h1+ h = h0' + (h1' - h0') * f++vector3d :: (NTransform a) => a -> Vector3d+vector3d = vec . pos . toNVector++angular :: Vector3d -> Length -> AngularPosition NVector+angular v = nvectorHeight (nvector (vx v) (vy v) (vz v))++ll :: (NTransform a) => a -> LatLong+ll = nvectorToLatLong . pos . toNVector
+ src/Data/Geo/Jord/Kinematics.hs view
@@ -0,0 +1,444 @@+-- | +-- Module: Data.Geo.Jord.Kinematics +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with kinematics calculations assuming a __spherical__ earth model. +-- +-- All functions are implemented using the vector-based approached described in +-- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation> +-- and in <https://calhoun.nps.edu/bitstream/handle/10945/29516/sometacticalalgo00shud.pdf Shudde, Rex H. (1986). Some tactical algorithms for spherical geometry> +-- +module Data.Geo.Jord.Kinematics+ (+ -- * The 'Track' type.+ Track(..)+ -- * The 'Course' type. + , Course+ -- * The 'Cpa' type. + , Cpa+ , cpaTime+ , cpaDistance+ , cpaPosition1+ , cpaPosition2+ -- * The 'Intercept' type. + , Intercept+ , interceptTime+ , interceptDistance+ , interceptPosition+ , interceptorBearing+ , interceptorSpeed+ -- * Calculations + , course+ , position+ , position84+ , cpa+ , cpa84+ , intercept+ , intercept84+ , interceptBySpeed+ , interceptBySpeed84+ , interceptByTime+ , interceptByTime84+ ) where++import Control.Applicative+import Data.Geo.Jord.Angle+import Data.Geo.Jord.AngularPosition+import Data.Geo.Jord.Duration+import Data.Geo.Jord.Earth+import Data.Geo.Jord.Geodetics+import Data.Geo.Jord.LatLong+import Data.Geo.Jord.Length+import Data.Geo.Jord.NVector+import Data.Geo.Jord.Quantity+import Data.Geo.Jord.Speed+import Data.Geo.Jord.Transformation+import Data.Geo.Jord.Vector3d++-- | 'Track' represents the state of a vehicle by its current position, bearing and speed. +data Track a = Track+ { trackPos :: a -- ^ position of the track. + , trackBearing :: Angle -- ^ bearing of the track. + , trackSpeed :: Speed -- ^ speed of the track. + } deriving (Eq, Show)++-- | 'GreatCircle' from track. +instance (NTransform a, Show a) => IsGreatCircle (Track a) where+ greatCircleE t = greatCircleE (trackPos t, trackBearing t)++-- | 'Course' represents the cardinal direction in which the vehicle is to be steered. +newtype Course =+ Course Vector3d+ deriving (Eq, Show)++instance IsVector3d Course where+ vec (Course v) = v++-- | Time to, and distance at, closest point of approach (CPA) as well as position of both tracks at CPA. +data Cpa a = Cpa+ { cpaTime :: Duration -- ^ time to CPA. + , cpaDistance :: Length -- ^ distance at CPA. + , cpaPosition1 :: a -- ^ position of track 1 at CPA. + , cpaPosition2 :: a -- ^ position of track 2 at CPA. + } deriving (Eq, Show)++-- | Time, distance and position of intercept as well as speed and initial bearing of interceptor. +data Intercept a = Intercept+ { interceptTime :: Duration -- ^ time to intercept. + , interceptDistance :: Length -- ^ distance at intercept. + , interceptPosition :: a -- ^ position of intercept. + , interceptorBearing :: Angle -- ^ initial bearing of interceptor. + , interceptorSpeed :: Speed -- ^ speed of interceptor. + } deriving (Eq, Show)++-- | @course p b@ computes the course of a vehicle currently at position @p@ and following bearing @b@. +course :: (NTransform a) => a -> Angle -> Course+course p b = Course (Vector3d (vz (head r)) (vz (r !! 1)) (vz (r !! 2)))+ where+ ll = nvectorToLatLong . pos . toNVector $ p+ lat = latitude ll+ lon = longitude ll+ r = mdot (mdot (rz (negate' lon)) (ry lat)) (rx b)++-- | @position t d r@ computes the position of a track @t@ after duration @d@ has elapsed and using the earth radius @r@. +-- +-- @ +-- let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000) +-- let b = decimalDegrees 96.0217 +-- let s = kilometresPerHour 124.8 +-- let p1 = decimalLatLongHeight 53.1882691 0.1332741 (metres 15000) +-- position (Track p0 b s) (hours 1) r84 = p1 +-- @ +position :: (NTransform a) => Track a -> Duration -> Length -> a+position (Track p0 b s) d = position' p0 s (course p0 b) (toSeconds d)++-- | 'position' using the mean radius of the WGS84 reference ellipsoid. +position84 :: (NTransform a) => Track a -> Duration -> a+position84 t d = position t d r84++-- | @cpa t1 t2 r@ computes the closest point of approach between tracks @t1@ and @t2@ and using the earth radius @r@. +-- +-- @ +-- let p1 = decimalLatLong 20 (-60) +-- let b1 = decimalDegrees 10 +-- let s1 = knots 15 +-- let p2 = decimalLatLong 34 (-50) +-- let b2 = decimalDegrees 220 +-- let s2 = knots 300 +-- let t1 = Track p1 b1 s1 +-- let t2 = Track p2 b2 s2 +-- let c = cpa t1 t2 r84 +-- fmap cpaTime c = Just (milliseconds 11396155) +-- fmap cpaDistance c = Just (kilometres 124.2317453) +-- @ +cpa :: (Eq a, NTransform a) => Track a -> Track a -> Length -> Maybe (Cpa a)+cpa (Track p1 b1 s1) (Track p2 b2 s2) r+ | p1 == p2 = Just (Cpa zero zero p1 p2)+ | t < 0 = Nothing+ | otherwise = Just (Cpa (seconds t) d cp1 cp2)+ where+ c1 = course p1 b1+ c2 = course p2 b2+ t = timeToCpa p1 c1 s1 p2 c2 s2 r+ cp1 = position' p1 s1 c1 t r+ cp2 = position' p2 s2 c2 t r+ d = surfaceDistance cp1 cp2 r++-- | 'cpa' using the mean radius of the WGS84 reference ellipsoid. +cpa84 :: (Eq a, NTransform a) => Track a -> Track a -> Maybe (Cpa a)+cpa84 t1 t2 = cpa t1 t2 r84++-- | @intercept t p r@ computes the __minimum__ speed of interceptor at +-- position @p@ needed for an intercept with target track @t@ to take place +-- using the earth radius @r@. Intercept time, position, distance and interceptor +-- bearing are derived from this minimum speed. Returns 'Nothing' if intercept +-- cannot be achieved e.g.: +-- +-- * interceptor and target are at the same position +-- +-- * interceptor is on the great circle of target and behind as the minimum speed would be target speed + epsillon +-- +-- @ +-- let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600) +-- let ip = (decimalLatLong 20 (-60)) +-- let i = intercept t ip r84 +-- fmap interceptorSpeed i = Just (knots 52.837096) +-- fmap interceptTime i = Just (seconds 5947.698) +-- @ +intercept :: (Eq a, NTransform a) => Track a -> a -> Length -> Maybe (Intercept a)+intercept t@(Track tp tb ts) p r = interceptByTime t p (seconds d) r+ where+ ct0 = course tp tb+ d = timeToIntercept tp ts ct0 p r++-- | 'intercept' using the mean radius of the WGS84 reference ellipsoid. +intercept84 :: (Eq a, NTransform a) => Track a -> a -> Maybe (Intercept a)+intercept84 t p = intercept t p r84++-- | @interceptBySpeed t p s r@ computes the time needed by interceptor at+-- position @p@ and travelling at speed @s@ to intercept target track @t@+-- using the earth radius @r@. Returns 'Nothing' if intercept+-- cannot be achieved e.g.:+--+-- * interceptor and target are at the same position+--+-- * interceptor speed is below minimum speed+interceptBySpeed :: (Eq a, NTransform a) => Track a -> a -> Speed -> Length -> Maybe (Intercept a)+interceptBySpeed t@(Track tp tb ts) p s r = interceptByTime t p (seconds d) r+ where+ ct0 = course tp tb+ d = timeToInterceptSpeed tp ts ct0 p s r++-- | 'interceptBySpeed' using the mean radius of the WGS84 reference ellipsoid.+interceptBySpeed84 :: (Eq a, NTransform a) => Track a -> a -> Speed -> Maybe (Intercept a)+interceptBySpeed84 t p s = interceptBySpeed t p s r84++-- | @interceptByTime t p d r@ computes the speed of interceptor at +-- position @p@ needed for an intercept with target track @t@ to take place +-- after duration @d@ and using the earth radius @r@. Returns 'Nothing' if +-- given duration is <= 0 or interceptor and target are at the same position. +-- +-- @ +-- let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600) +-- let ip = (decimalLatLong 20 (-60)) +-- let d = seconds 2700 +-- let i = interceptByTime t ip d r84 +-- fmap interceptorSpeed i = Just (knots 730.959238) +-- fmap interceptorBearing i = Just (decimalDegrees 26.1199030) +-- fmap interceptPosition i = Just (decimalLatLong 28.1366797 (-55.4559475)) +-- fmap interceptDistance i = Just (metres 1015302.3815) +-- fmap interceptTime i = Just (seconds 2700) +-- @ +interceptByTime :: (Eq a, NTransform a) => Track a -> a -> Duration -> Length -> Maybe (Intercept a)+interceptByTime t p d r+ | toMilliseconds d <= 0 = Nothing+ | trackPos t == p = Nothing+ | otherwise = fmap (\b -> Intercept d idist ipos b is) ib+ where+ ipos = position t d r+ idist = surfaceDistance p ipos r+ ib = initialBearing p ipos <|> initialBearing p (trackPos t)+ is = metresPerSecond (toMetres idist / toSeconds d)++-- | 'interceptByTime' using the mean radius of the WGS84 reference ellipsoid. +interceptByTime84 :: (Eq a, NTransform a) => Track a -> a -> Duration -> Maybe (Intercept a)+interceptByTime84 t p d = interceptByTime t p d r84++-- | position from speed course and seconds. +position' :: (NTransform a) => a -> Speed -> Course -> Double -> Length -> a+position' p0 s c sec r = fromNVector (nvectorHeight (nvector (vx v1) (vy v1) (vz v1)) h0)+ where+ nv0 = toNVector p0+ v0 = vec . pos $nv0+ h0 = height nv0+ v1 = position'' v0 s (vec c) sec r++-- | position from speed course and seconds. +position'' :: Vector3d -> Speed -> Vector3d -> Double -> Length -> Vector3d+position'' v0 s c sec r = v1+ where+ w = toMetresPerSecond s / toMetres r+ v1 = vadd (vscale v0 (cos (w * sec))) (vscale c (sin (w * sec)))++-- | time to CPA. +timeToCpa :: (NTransform a) => a -> Course -> Speed -> a -> Course -> Speed -> Length -> Double+timeToCpa p1 c1 s1 p2 c2 s2 r = cpaNrRec v10 c10 w1 v20 c20 w2 0 0+ where+ v10 = vec . pos . toNVector $ p1+ c10 = vec c1+ rm = toMetres r+ w1 = toMetresPerSecond s1 / rm+ v20 = vec . pos . toNVector $ p2+ c20 = vec c2+ w2 = toMetresPerSecond s2 / rm++-- | time to intercept with minimum speed. +timeToIntercept :: (NTransform a) => a -> Speed -> Course -> a -> Length -> Double+timeToIntercept p2 s2 c20 p1 r = intMinNrRec v10 v20 (vec c20) s2 w2 r s0 t0 0+ where+ v10 = vec . pos . toNVector $ p1+ v20 = vec . pos . toNVector $ p2+ s2mps = toMetresPerSecond s2+ rm = toMetres r+ w2 = s2mps / rm+ s0 = ad v10 v20+ t0 = rm * s0 / s2mps++-- | time to intercept with speed.+timeToInterceptSpeed :: (NTransform a) => a -> Speed -> Course -> a -> Speed -> Length -> Double+timeToInterceptSpeed p2 s2 c20 p1 s1 r = intSpdNrRec v10 w1 v20 (vec c20) s2 w2 r s0 t0 0+ where+ v10 = vec . pos . toNVector $ p1+ v20 = vec . pos . toNVector $ p2+ s1mps = toMetresPerSecond s1+ s2mps = toMetresPerSecond s2+ rm = toMetres r+ w2 = s2mps / rm+ w1 = s1mps / rm+ s0 = ad v10 v20+ t0 = rm * s0 / s2mps++rx :: Angle -> [Vector3d]+rx a = [Vector3d 1 0 0, Vector3d 0 c s, Vector3d 0 (-s) c]+ where+ c = cos' a+ s = sin' a++ry :: Angle -> [Vector3d]+ry a = [Vector3d c 0 (-s), Vector3d 0 1 0, Vector3d s 0 c]+ where+ c = cos' a+ s = sin' a++rz :: Angle -> [Vector3d]+rz a = [Vector3d c s 0, Vector3d (-s) c 0, Vector3d 0 0 1]+ where+ c = cos' a+ s = sin' a++cpaA :: Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double+cpaA v10 c10 w1 v20 c20 w2 = negate (vdot (vscale v10 w1) c20 + vdot (vscale v20 w2) c10)++cpaB :: Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double+cpaB v10 c10 w1 v20 c20 w2 = vdot (vscale c10 w1) v20 + vdot (vscale c20 w2) v10++cpaC :: Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double+cpaC v10 c10 w1 v20 c20 w2 = negate (vdot (vscale v10 w1) v20 - vdot (vscale c20 w2) c10)++cpaD :: Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double+cpaD v10 c10 w1 v20 c20 w2 = vdot (vscale c10 w1) c20 - vdot (vscale v20 w2) v10++cpaFt :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double+cpaFt cw1t cw2t sw1t sw2t a b c d =+ a * sw1t * sw2t + b * cw1t * cw2t + c * sw1t * cw2t + d * cw1t * sw2t++cpaDft ::+ Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+ -> Double+cpaDft w1 w2 cw1t cw2t sw1t sw2t a b c d =+ negate ((c * w2 + d * w1) * sw1t * sw2t) + (d * w2 + c * w1) * cw1t * cw2t ++ (a * w2 - b * w1) * sw1t * cw2t -+ (b * w2 - a * w1) * cw1t * sw2t++cpaStep :: Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double -> Double+cpaStep v10 c10 w1 v20 c20 w2 t =+ cpaFt cw1t cw2t sw1t sw2t a b c d / cpaDft w1 w2 cw1t cw2t sw1t sw2t a b c d+ where+ cw1t = cos (w1 * t)+ cw2t = cos (w2 * t)+ sw1t = sin (w1 * t)+ sw2t = sin (w2 * t)+ a = cpaA v10 c10 w1 v20 c20 w2+ b = cpaB v10 c10 w1 v20 c20 w2+ c = cpaC v10 c10 w1 v20 c20 w2+ d = cpaD v10 c10 w1 v20 c20 w2++-- | Newton-Raphson for CPA time. +-- note: this should always converge to the minimum time given +-- that the assumptions made in the proof of quadratic convergence are met +cpaNrRec ::+ Vector3d -> Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double -> Int -> Double+cpaNrRec v10 c10 w1 v20 c20 w2 ti i+ | i == 50 = -1.0 -- no convergence + | abs fi < 1e-12 = ti1+ | otherwise = cpaNrRec v10 c10 w1 v20 c20 w2 ti1 (i + 1)+ where+ fi = cpaStep v10 c10 w1 v20 c20 w2 ti+ ti1 = ti - fi++-- | Newton-Raphson for min speed intercept. +-- note: this should always converge to the minimum time given +-- that the assumptions made in the proof of quadratic convergence are met +intMinNrRec ::+ Vector3d+ -> Vector3d+ -> Vector3d+ -> Speed+ -> Double+ -> Length+ -> Double+ -> Double+ -> Int+ -> Double+intMinNrRec v10 v20 c20 s2 w2 r si ti i+ | i == 50 = -1.0 -- no convergence + | abs fi < 1e-12 = ti1+ | otherwise = intMinNrRec v10 v20 c20 s2 w2 r si1 ti1 (i + 1)+ where+ fi = intMinStep v10 v20 c20 w2 si ti+ ti1 = ti - fi+ v2t = position'' v20 s2 c20 ti1 r+ si1 = ad v10 v2t++intMinStep :: Vector3d -> Vector3d -> Vector3d -> Double -> Double -> Double -> Double+intMinStep v10 v20 c20 w2 s t =+ dsdt s w2 v10v20 v10c20 sinw2t cosw2t / d2sdt2 s w2 v10v20 v10c20 sinw2t cosw2t+ where+ cosw2t = cos (w2 * t)+ sinw2t = sin (w2 * t)+ v10v20 = vdot v10 v20+ v10c20 = vdot v10 c20++-- | Newton-Raphson for speed intercept.+-- note: this should always converge to the minimum time given+-- that the assumptions made in the proof of quadratic convergence are met+intSpdNrRec ::+ Vector3d+ -> Double+ -> Vector3d+ -> Vector3d+ -> Speed+ -> Double+ -> Length+ -> Double+ -> Double+ -> Int+ -> Double+intSpdNrRec v10 w1 v20 c20 s2 w2 r si ti i+ | i == 50 = -1.0 -- no convergence+ | abs fi < 1e-12 = ti1+ | otherwise = intSpdNrRec v10 w1 v20 c20 s2 w2 r si1 ti1 (i + 1)+ where+ fi = intSpdStep v10 w1 v20 c20 w2 si ti+ ti1 = ti - fi+ v2t = position'' v20 s2 c20 ti1 r+ si1 = ad v10 v2t++intSpdStep :: Vector3d -> Double -> Vector3d -> Vector3d -> Double -> Double -> Double -> Double+intSpdStep v10 w1 v20 c20 w2 s t = f / df+ where+ cosw2t = cos (w2 * t)+ sinw2t = sin (w2 * t)+ v10v20 = vdot v10 v20+ v10c20 = vdot v10 c20+ f = s / t - w1+ df = (1.0 / t) * (dsdt s w2 v10v20 v10c20 sinw2t cosw2t - s / t)++dsdt :: Double -> Double -> Double -> Double -> Double -> Double -> Double+dsdt s w2 v10v20 v10c20 sinw2t cosw2t =+ ((-1.0) / sin s) * ((-w2) * (v10v20 * sinw2t - v10c20 * cosw2t))++d2sdt2 :: Double -> Double -> Double -> Double -> Double -> Double -> Double+d2sdt2 s w2 v10v20 v10c20 sinw2t cosw2t =+ ((-1.0) / sin s) * (cos s / (sins * sins) * x10d2x2dt2 * x10d2x2dt2 + x10d2x2dt2)+ where+ sins = sin s+ x10d2x2dt2 = negate (w2 * w2) * (v10v20 * cosw2t + v10c20 * sinw2t)++-- | angle in radians between 2 n-vectors (as vector3d), copied from Geodetics +-- without the sign and returing radians. +ad :: Vector3d -> Vector3d -> Double+ad v1 v2 = atan2 (vnorm (vcross v1 v2)) (vdot v1 v2)
src/Data/Geo/Jord/Length.hs view
@@ -6,122 +6,122 @@ -- Stability: experimental -- Portability: portable -- --- Types and functions for working with (signed) lengths in metres, kilometres or nautical miles. +-- Types and functions for working with (signed) lengths in metres, kilometres, nautical miles or feet. -- -module Data.Geo.Jord.Length- (- -- * The 'Length' type- Length(millimetres)+module Data.Geo.Jord.Length + ( + -- * The 'Length' type + Length -- * Smart constructors - , feet- , kilometres- , metres- , nauticalMiles+ , feet + , kilometres + , metres + , nauticalMiles -- * Read - , readLength- , readLengthE- , readLengthF+ , readLength + , readLengthE + , readLengthF -- * Conversions - , toFeet- , toKilometres- , toMetres- , toNauticalMiles- ) where--import Control.Applicative-import Control.Monad.Fail-import Data.Geo.Jord.Parse-import Data.Geo.Jord.Quantity-import Prelude hiding (fail, length)-import Text.ParserCombinators.ReadP-import Text.Read hiding (pfail)---- | A length with a resolution of 1 millimetre. -newtype Length = Length- { millimetres :: Int- } deriving (Eq)-+ , toFeet + , toKilometres + , toMetres + , toNauticalMiles + ) where + +import Control.Applicative +import Control.Monad.Fail +import Data.Geo.Jord.Parse +import Data.Geo.Jord.Quantity +import Prelude hiding (fail, length) +import Text.ParserCombinators.ReadP +import Text.Read hiding (pfail) + +-- | A length with a resolution of 0.1 millimetre. +newtype Length = Length + { tenthOfMm :: Int + } deriving (Eq) + -- | See 'readLength'. -instance Read Length where- readsPrec _ = readP_to_S length---- | Length is shown in metres when <= 10,000 m and in kilometres otherwise. -instance Show Length where- show l- | m <= 10000.0 = show m ++ "m"- | otherwise = show (m / 1000.0) ++ "km"- where- m = toMetres l---- | Add/Subtract Length. -instance Quantity Length where- add a b = Length (millimetres a + millimetres b)- sub a b = Length (millimetres a - millimetres b)- zero = Length 0---- | 'Length' from given amount of feet.-feet :: Double -> Length-feet ft = metres (ft * 0.3048)---- | 'Length' from given amount of kilometres.-kilometres :: Double -> Length-kilometres km = metres (km * 1000.0)-+instance Read Length where + readsPrec _ = readP_to_S length + +-- | Length is shown in metres when absolute value is <= 10,000 m and in kilometres otherwise. +instance Show Length where + show l + | abs m <= 10000.0 = show m ++ "m" + | otherwise = show (m / 1000.0) ++ "km" + where + m = toMetres l + +-- | Add/Subtract 'Length's. +instance Quantity Length where + add a b = Length (tenthOfMm a + tenthOfMm b) + sub a b = Length (tenthOfMm a - tenthOfMm b) + zero = Length 0 + +-- | 'Length' from given amount of feet. +feet :: Double -> Length +feet ft = Length (round (ft * 3048.0)) + +-- | 'Length' from given amount of kilometres. +kilometres :: Double -> Length +kilometres km = Length (round (km * 10000000.0)) + -- | 'Length' from given amount of metres. -metres :: Double -> Length-metres m = Length (round (m * 1000.0))---- | 'Length' from given amount of nautical miles.-nauticalMiles :: Double -> Length-nauticalMiles nm = metres (nm * 1852.0)---- | Obtains a 'Length' from the given string formatted as (-)float[m|km|nm] - e.g. 3000m, 2.5km or -154nm. +metres :: Double -> Length +metres m = Length (round (m * 10000.0)) + +-- | 'Length' from given amount of nautical miles. +nauticalMiles :: Double -> Length +nauticalMiles nm = Length (round (nm * 18520000.0)) + +-- | Obtains a 'Length' from the given string formatted as (-)float[m|km|nm|ft] - e.g. 3000m, 2.5km, -154nm or 10000ft. -- -- This simply calls @read s :: Length@ so 'error' should be handled at the call site. -- -readLength :: String -> Length-readLength s = read s :: Length-+readLength :: String -> Length +readLength s = read s :: Length + -- | Same as 'readLength' but returns a 'Either'. -readLengthE :: String -> Either String Length-readLengthE s =- case readMaybe s of- Nothing -> Left ("couldn't read length " ++ s)- Just l -> Right l-+readLengthE :: String -> Either String Length +readLengthE s = + case readMaybe s of + Nothing -> Left ("couldn't read length " ++ s) + Just l -> Right l + -- | Same as 'readLength' but returns a 'MonadFail'. -readLengthF :: (MonadFail m) => String -> m Length-readLengthF s =- let p = readEither s- in case p of- Left e -> fail e- Right l -> return l---- | @toFeet l@ converts @l@ to feet.-toFeet :: Length -> Double-toFeet l = toMetres l / 0.3048-+readLengthF :: (MonadFail m) => String -> m Length +readLengthF s = + let p = readEither s + in case p of + Left e -> fail e + Right l -> return l + +-- | @toFeet l@ converts @l@ to feet. +toFeet :: Length -> Double +toFeet (Length l) = fromIntegral l / 3048.0 + -- | @toKilometres l@ converts @l@ to kilometres. -toKilometres :: Length -> Double-toKilometres l = toMetres l / 1000.0-+toKilometres :: Length -> Double +toKilometres (Length l) = fromIntegral l / 10000000.0 + -- | @toMetres l@ converts @l@ to metres. -toMetres :: Length -> Double-toMetres (Length mm) = fromIntegral mm / 1000.0-+toMetres :: Length -> Double +toMetres (Length l) = fromIntegral l / 10000.0 + -- | @toNauticalMiles l@ converts @l@ to nautical miles. -toNauticalMiles :: Length -> Double-toNauticalMiles l = toMetres l / 1852.0-+toNauticalMiles :: Length -> Double +toNauticalMiles (Length l) = fromIntegral l / 18520000.0 + -- | Parses and returns a 'Length'. -length :: ReadP Length-length = do- v <- number- skipSpaces- u <- string "m" <|> string "km" <|> string "Nm" <|> string "ft"- case u of- "m" -> return (metres v)- "km" -> return (kilometres v)- "Nm" -> return (nauticalMiles v)- "ft" -> return (feet v)- _ -> pfail+length :: ReadP Length +length = do + v <- number + skipSpaces + u <- string "m" <|> string "km" <|> string "Nm" <|> string "ft" + case u of + "m" -> return (metres v) + "km" -> return (kilometres v) + "Nm" -> return (nauticalMiles v) + "ft" -> return (feet v) + _ -> pfail
src/Data/Geo/Jord/NVector.hs view
@@ -8,34 +8,34 @@ -- -- Types and functions for working with n-vectors. -- -module Data.Geo.Jord.NVector- ( NVector- , nvector- , northPole- , southPole- ) where--import Data.Geo.Jord.Vector3d-+module Data.Geo.Jord.NVector + ( NVector + , nvector + , northPole + , southPole + ) where + +import Data.Geo.Jord.Vector3d + -- | Represents a position as the normal vector to the sphere. -- -- Orientation: z-axis points to the North Pole along the Earth's rotation axis, -- x-axis points towards the point where latitude = longitude = 0. -newtype NVector =- NVector Vector3d- deriving (Eq, Show)--instance IsVector3d NVector where- vec (NVector v) = v-+newtype NVector = + NVector Vector3d + deriving (Eq, Show) + +instance IsVector3d NVector where + vec (NVector v) = v + -- | Unit 'NVector' from given x, y and z. -nvector :: Double -> Double -> Double -> NVector-nvector x y z = NVector (vunit (Vector3d x y z))---- | Horizontal position of the North Pole. -northPole :: NVector-northPole = NVector (Vector3d 0.0 0.0 1.0)-+nvector :: Double -> Double -> Double -> NVector +nvector x y z = NVector (vunit (Vector3d x y z)) + -- | Horizontal position of the North Pole. -southPole :: NVector-southPole = NVector (Vector3d 0.0 0.0 (-1.0))+northPole :: NVector +northPole = NVector (Vector3d 0.0 0.0 1.0) + +-- | Horizontal position of the South Pole. +southPole :: NVector +southPole = NVector (Vector3d 0.0 0.0 (-1.0))
src/Data/Geo/Jord/Parse.hs view
@@ -1,52 +1,52 @@--- |--- Module: Data.Geo.Jord.Parse--- Copyright: (c) 2018 Cedric Liegeois--- License: BSD3--- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr>--- Stability: experimental--- Portability: portable------ internal 'ReadP' parsers used by "Jord".----module Data.Geo.Jord.Parse- ( digits- , double- , integer- , natural- , number- ) where--import Control.Applicative-import Data.Char-import Text.ParserCombinators.ReadP---- | Parses the given number of digits and returns the read 'Int'.-digits :: Int -> ReadP Int-digits n = fmap read (count n digit)---- | Parses optionally a @-@ followed by a 'positive'.'positive' and returns the read 'Double'.-double :: ReadP Double-double = do- s <- option 1.0 (fmap (\_ -> -1.0) (char '-'))- i <- natural- f <- char '.' >> natural- return (s * (read (show i ++ "." ++ show f) :: Double))---- | Parses optionally a @-@ followed by a 'positive' and returns the read 'Int'.-integer :: ReadP Int-integer = do- s <- option 1 (fmap (\_ -> -1) (char '-'))- p <- natural- return (s * p)---- | Parses 1 or more 'digit's and returns the read 'Int'.-natural :: ReadP Int-natural = fmap read (munch1 isDigit)---- | Parses an 'integer' or 'double' and returns the read 'Double'.-number :: ReadP Double-number = double <|> fmap fromIntegral integer---- | Parses and returns a digit.-digit :: ReadP Char-digit = satisfy isDigit+-- | +-- Module: Data.Geo.Jord.Parse +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- internal 'ReadP' parsers used by "Jord". +-- +module Data.Geo.Jord.Parse + ( digits + , double + , integer + , natural + , number + ) where + +import Control.Applicative +import Data.Char +import Text.ParserCombinators.ReadP + +-- | Parses the given number of digits and returns the read 'Int'. +digits :: Int -> ReadP Int +digits n = fmap read (count n digit) + +-- | Parses optionally a @-@ followed by a 'positive'.'positive' and returns the read 'Double'. +double :: ReadP Double +double = do + s <- option 1.0 (fmap (\_ -> -1.0) (char '-')) + i <- natural + f <- char '.' >> natural + return (s * (read (show i ++ "." ++ show f) :: Double)) + +-- | Parses optionally a @-@ followed by a 'positive' and returns the read 'Int'. +integer :: ReadP Int +integer = do + s <- option 1 (fmap (\_ -> -1) (char '-')) + p <- natural + return (s * p) + +-- | Parses 1 or more 'digit's and returns the read 'Int'. +natural :: ReadP Int +natural = fmap read (munch1 isDigit) + +-- | Parses an 'integer' or 'double' and returns the read 'Double'. +number :: ReadP Double +number = double <|> fmap fromIntegral integer + +-- | Parses and returns a digit. +digit :: ReadP Char +digit = satisfy isDigit
src/Data/Geo/Jord/Quantity.hs view
@@ -8,11 +8,11 @@ -- -- Classes for working with quantities. -- -module Data.Geo.Jord.Quantity- ( Quantity(..)- ) where-+module Data.Geo.Jord.Quantity + ( Quantity(..) + ) where + -- | Something that can be added or subtracted. -class (Eq a) => Quantity a where- add, sub :: a -> a -> a- zero :: a+class (Eq a) => Quantity a where + add, sub :: a -> a -> a + zero :: a
src/Data/Geo/Jord/Rotation.hs view
@@ -78,7 +78,7 @@ -- that the relation between a vector v decomposed in A and B is given by: -- @v_A = mdot R_AB v_B@ -- --- The rotation matrix R_AB is created based on 3 angles x,y,z about new axes +-- The rotation matrix R_AB is created based on 3 'Angle's x,y,z about new axes -- (intrinsic) in the order x-y-z. The angles are called Euler angles or -- Tait-Bryan angles and are defined by the following procedure of successive -- rotations: @@ -109,7 +109,7 @@ -- that the relation between a vector v decomposed in A and B is given by: -- @v_A = mdot R_AB v_B@ -- --- The rotation matrix R_AB is created based on 3 angles +-- The rotation matrix R_AB is created based on 3 'Angle's -- z,y,x about new axes (intrinsic) in the order z-y-x. The angles are called -- Euler angles or Tait-Bryan angles and are defined by the following -- procedure of successive rotations:
+ src/Data/Geo/Jord/Speed.hs view
@@ -0,0 +1,134 @@+-- | +-- Module: Data.Geo.Jord.Speed +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with speed in metres per second, kilometres per hour, miles per hour, knots or feet per second. +-- +module Data.Geo.Jord.Speed + ( + -- * The 'Speed' type + Speed + -- * Smart constructors + , metresPerSecond + , kilometresPerHour + , milesPerHour + , knots + , feetPerSecond + -- * Read + , readSpeed + , readSpeedE + , readSpeedF + -- * Conversions + , toMetresPerSecond + , toKilometresPerHour + , toMilesPerHour + , toKnots + , toFeetPerSecond + ) where + +import Control.Applicative +import Control.Monad.Fail +import Data.Geo.Jord.Parse +import Data.Geo.Jord.Quantity +import Prelude hiding (fail) +import Text.ParserCombinators.ReadP +import Text.Read hiding (pfail) + +-- | A speed with a resolution of 1 millimetre per hour. +newtype Speed = Speed + { mmPerHour :: Int + } deriving (Eq) + +-- | See 'readSpeed'. +instance Read Speed where + readsPrec _ = readP_to_S speed + +-- | Speed is shown in kilometres per hour. +instance Show Speed where + show s = show (toKilometresPerHour s) ++ "km/h" + +-- | Add/Subtract Speed. +instance Quantity Speed where + add a b = Speed (mmPerHour a + mmPerHour b) + sub a b = Speed (mmPerHour a - mmPerHour b) + zero = Speed 0 + +-- | 'Speed' from given amount of metres per second. +metresPerSecond :: Double -> Speed +metresPerSecond mps = Speed (round (mps * 3600000.0)) + +-- | 'Speed' from given amount of kilometres per hour. +kilometresPerHour :: Double -> Speed +kilometresPerHour kph = Speed (round (kph * 1e+6)) + +-- | 'Speed' from given amount of miles per hour. +milesPerHour :: Double -> Speed +milesPerHour mph = Speed (round (mph * 1609344.0)) + +-- | 'Speed' from given amount of knots. +knots :: Double -> Speed +knots kt = Speed (round (kt * 1852000.0)) + +-- | 'Speed' from given amount of feet per second. +feetPerSecond :: Double -> Speed +feetPerSecond fps = Speed (round (fps * 1097280.0)) + +-- | Obtains a 'Speed' from the given string formatted as (-)float[m/s|km/h|mph|kt] - e.g. 300m/s, 250km/h, -154mph, 400kt or 100ft/s. +-- +-- This simply calls @read s :: Speed@ so 'error' should be handled at the call site. +-- +readSpeed :: String -> Speed +readSpeed s = read s :: Speed + +-- | Same as 'readSpeed' but returns a 'Either'. +readSpeedE :: String -> Either String Speed +readSpeedE s = + case readMaybe s of + Nothing -> Left ("couldn't read speed " ++ s) + Just l -> Right l + +-- | Same as 'readSpeed' but returns a 'MonadFail'. +readSpeedF :: (MonadFail m) => String -> m Speed +readSpeedF s = + let p = readEither s + in case p of + Left e -> fail e + Right l -> return l + +-- | @toMetresPerSecond s@ converts @s@ to metres per second. +toMetresPerSecond :: Speed -> Double +toMetresPerSecond (Speed s) = fromIntegral s / 3600000.0 + +-- | @toKilometresPerHour s@ converts @s@ to kilometres per hour. +toKilometresPerHour :: Speed -> Double +toKilometresPerHour (Speed s) = fromIntegral s / 1e+6 + +-- | @toMilesPerHour s@ converts @s@ to miles per hour. +toMilesPerHour :: Speed -> Double +toMilesPerHour (Speed s) = fromIntegral s / 1609344.0 + +-- | @toKnots s@ converts @s@ to knots. +toKnots :: Speed -> Double +toKnots (Speed s) = fromIntegral s / 1852000.0 + +-- | @toFeetPerSecond s@ converts @s@ to feet per second. +toFeetPerSecond :: Speed -> Double +toFeetPerSecond (Speed s) = fromIntegral s / 1097280.0 + +-- | Parses and returns a 'Speed'. +speed :: ReadP Speed +speed = do + s <- number + skipSpaces + u <- string "m/s" <|> string "km/h" <|> string "mph" <|> string "kt" <|> string "ft/s" + case u of + "m/s" -> return (metresPerSecond s) + "km/h" -> return (kilometresPerHour s) + "mph" -> return (milesPerHour s) + "kt" -> return (knots s) + "ft/s" -> return (feetPerSecond s) + _ -> pfail
− src/Data/Geo/Jord/Transform.hs
@@ -1,192 +0,0 @@-{-# LANGUAGE FlexibleInstances #-}---- | --- Module: Data.Geo.Jord.Transform --- Copyright: (c) 2018 Cedric Liegeois --- License: BSD3 --- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> --- Stability: experimental --- Portability: portable --- --- Transformations between coordinates systems both in spherical and ellipsoidal form. --- --- All functions are implemented using the vector-based approached described in --- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation> --- --- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> --- -module Data.Geo.Jord.Transform- ( NTransform(..)- , ETransform(..)- , nvectorToLatLong- , latLongToNVector- , ecefToNVector- , nvectorToEcef- , geodeticHeight- ) where--import Data.Geo.Jord.Angle-import Data.Geo.Jord.AngularPosition-import Data.Geo.Jord.Earth-import Data.Geo.Jord.EcefPosition-import Data.Geo.Jord.LatLong-import Data.Geo.Jord.Length-import Data.Geo.Jord.NVector-import Data.Geo.Jord.Quantity-import Data.Geo.Jord.Vector3d---- | Transformation between positions and 'AngularPosition' of 'NVector'. -class NTransform a where- toNVector :: a -> AngularPosition NVector -- ^ position to 'AngularPosition' of 'NVector'. - fromNVector :: AngularPosition NVector -> a -- ^ 'AngularPosition' of 'NVector' and height to position. ---- | 'NVector' <-> 'AngularPosition' of 'NVector'.-instance NTransform NVector where- toNVector nv = AngularPosition nv zero- fromNVector = pos---- | 'LatLong' <-> 'AngularPosition' of 'NVector'. -instance NTransform LatLong where- toNVector ll = AngularPosition (latLongToNVector ll) zero- fromNVector = nvectorToLatLong . pos---- | 'NTransform' identity. -instance NTransform (AngularPosition NVector) where- toNVector = id- fromNVector = id---- | 'AngularPosition' of 'LatLong' <-> 'AngularPosition' of 'NVector'. -instance NTransform (AngularPosition LatLong) where- toNVector (AngularPosition ll h) = AngularPosition (latLongToNVector ll) h- fromNVector (AngularPosition nv h) = AngularPosition (nvectorToLatLong nv) h---- | Transformation between 'EcefPosition' and angular or n-vector positions. -class ETransform a where- toEcef :: a -> Earth -> EcefPosition -- ^ position and earth model to to 'EcefPosition'. - fromEcef :: EcefPosition -> Earth -> a -- ^ 'EcefPosition' and earth model to position. ---- | 'NVector' <-> 'EcefPosition'. -instance ETransform NVector where- fromEcef p e = pos (ecefToNVector p e)- toEcef v = nvectorToEcef (nvectorHeight v zero)---- | 'LatLong' <-> 'EcefPosition'. -instance ETransform LatLong where- fromEcef p e = fromNVector (nvectorHeight (fromEcef p e :: NVector) zero)- toEcef = toEcef . toNVector---- | 'AngularPosition' of 'NVector' <-> 'EcefPosition'. -instance ETransform (AngularPosition NVector) where- fromEcef = ecefToNVector- toEcef = nvectorToEcef---- | 'AngularPosition' of 'LatLong' <-> 'EcefPosition'. -instance ETransform (AngularPosition LatLong) where- fromEcef p e = fromNVector (ecefToNVector p e)- toEcef = nvectorToEcef . toNVector---- | 'ETransform' identity.-instance ETransform EcefPosition where- fromEcef p _ = p- toEcef p _ = p---- | @nvectorToLatLong v@ transforms 'NVector' @v@ to an equivalent 'LatLong'. --- --- Same as 'toNVector'. -nvectorToLatLong :: NVector -> LatLong-nvectorToLatLong nv = latLong lat lon- where- v = vec nv- lat = atan2' (vz v) (sqrt (vx v * vx v + vy v * vy v))- lon = atan2' (vy v) (vx v)---- | @latLongToNVector ll@ transforms 'LatLong' @ll@ to an equivalent 'NVector'. --- --- See also 'fromNVector'. -latLongToNVector :: LatLong -> NVector-latLongToNVector ll = nvector x' y' z'- where- lat = latitude ll- lon = longitude ll- cl = cos' lat- x' = cl * cos' lon- y' = cl * sin' lon- z' = sin' lat---- | @ecefToNVector p e@ transforms 'EcefPosition' @p@ to an equivalent 'NVector' and geodetic height --- using earth model @e@. --- --- See also 'fromEcef' -ecefToNVector :: EcefPosition -> Earth -> AngularPosition NVector--- Ellipsoidal-ecefToNVector ep e@(Ellipsoidal el) = nvectorHeight (nvecEllipsoidal d e2 k px py pz) (metres h)- where- ev = vec ep- e' = eccentricity e- e2 = e' * e'- e4 = e2 * e2- a = toMetres (equatorialRadius el)- a2 = a * a- px = vx ev- py = vy ev- pz = vz ev- p = (px * px + py * py) / a2- q = ((1 - e2) / a2) * (pz * pz)- r = (p + q - e4) / 6.0- s = (e4 * p * q) / (4.0 * r * r * r)- t = (1.0 + s + sqrt (s * (2.0 + s))) ** (1 / 3)- u = r * (1.0 + t + 1.0 / t)- v = sqrt (u * u + q * e4)- w = e2 * (u + v - q) / (2.0 * v)- k = sqrt (u + v + w * w) - w- d = k * sqrt (px * px + py * py) / (k + e2)- h = ((k + e2 - 1.0) / k) * sqrt (d * d + pz * pz)--- Spherical-ecefToNVector p (Spherical r) = nvectorHeight (nvector (vx nv) (vy nv) (vz nv)) h- where- ev = vec p- nv = vunit ev- h = sub (metres (vnorm ev)) r--nvecEllipsoidal :: Double -> Double -> Double -> Double -> Double -> Double -> NVector-nvecEllipsoidal d e2 k px py pz = nvector nx' ny' nz'- where- s = 1.0 / sqrt (d * d + pz * pz)- a = k / (k + e2)- nx' = s * a * px- ny' = s * a * py- nz' = s * pz---- | @nvectorToEcef (n, h) e@ transforms 'NVector' @n@ and geodetic height @h@ --- to an equivalent 'EcefPosition' using earth model @e@. --- --- See also 'toEcef' -nvectorToEcef :: AngularPosition NVector -> Earth -> EcefPosition--- Ellipsoidal-nvectorToEcef (AngularPosition nv h) e@(Ellipsoidal el) = ecef ex' ey' ez'- where- v = vec nv- uv = vunit v- a = toMetres (equatorialRadius el)- b = toMetres (polarRadius e)- nx' = vx uv- ny' = vy uv- nz' = vz uv- m = (a * a) / (b * b)- n = b / sqrt ((nx' * nx' * m) + (ny' * ny' * m) + (nz' * nz'))- h' = toMetres h- ex' = metres (n * m * nx' + h' * nx')- ey' = metres (n * m * ny' + h' * ny')- ez' = metres (n * nz' + h' * nz')--- Spherical-nvectorToEcef (AngularPosition nv h) (Spherical r) = ecefMetres (vx ev) (vy ev) (vz ev)- where- unv = vunit . vec $ nv- n = add h r- ev = vscale unv (toMetres n)---- | @geodeticHeight p e@ computes the geodetic height of 'EcefPosition' @p@ using earth model @e@. --- --- The geodetic height (or ellipsoidal height) is __not__ the mean sea level (MSL) height. -geodeticHeight :: EcefPosition -> Earth -> Length-geodeticHeight p e = height (ecefToNVector p e)
+ src/Data/Geo/Jord/Transformation.hs view
@@ -0,0 +1,192 @@+{-# LANGUAGE FlexibleInstances #-} + +-- | +-- Module: Data.Geo.Jord.Transformation +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Transformations between coordinates systems both in spherical and ellipsoidal form. +-- +-- All functions are implemented using the vector-based approached described in +-- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation> +-- +-- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> +-- +module Data.Geo.Jord.Transformation + ( NTransform(..) + , ETransform(..) + , nvectorToLatLong + , latLongToNVector + , ecefToNVector + , nvectorToEcef + , geodeticHeight + ) where + +import Data.Geo.Jord.Angle +import Data.Geo.Jord.AngularPosition +import Data.Geo.Jord.Earth +import Data.Geo.Jord.EcefPosition +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Quantity +import Data.Geo.Jord.Vector3d + +-- | Transformation between positions and 'AngularPosition' of 'NVector'. +class NTransform a where + toNVector :: a -> AngularPosition NVector -- ^ position to 'AngularPosition' of 'NVector'. + fromNVector :: AngularPosition NVector -> a -- ^ 'AngularPosition' of 'NVector' and height to position. + +-- | 'NVector' <-> 'AngularPosition' of 'NVector'. +instance NTransform NVector where + toNVector nv = AngularPosition nv zero + fromNVector = pos + +-- | 'LatLong' <-> 'AngularPosition' of 'NVector'. +instance NTransform LatLong where + toNVector ll = AngularPosition (latLongToNVector ll) zero + fromNVector = nvectorToLatLong . pos + +-- | 'NTransform' identity. +instance NTransform (AngularPosition NVector) where + toNVector = id + fromNVector = id + +-- | 'AngularPosition' of 'LatLong' <-> 'AngularPosition' of 'NVector'. +instance NTransform (AngularPosition LatLong) where + toNVector (AngularPosition ll h) = AngularPosition (latLongToNVector ll) h + fromNVector (AngularPosition nv h) = AngularPosition (nvectorToLatLong nv) h + +-- | Transformation between 'EcefPosition' and angular or n-vector positions. +class ETransform a where + toEcef :: a -> Earth -> EcefPosition -- ^ position and earth model to to 'EcefPosition'. + fromEcef :: EcefPosition -> Earth -> a -- ^ 'EcefPosition' and earth model to position. + +-- | 'NVector' <-> 'EcefPosition'. +instance ETransform NVector where + fromEcef p e = pos (ecefToNVector p e) + toEcef v = nvectorToEcef (nvectorHeight v zero) + +-- | 'LatLong' <-> 'EcefPosition'. +instance ETransform LatLong where + fromEcef p e = fromNVector (nvectorHeight (fromEcef p e :: NVector) zero) + toEcef = toEcef . toNVector + +-- | 'AngularPosition' of 'NVector' <-> 'EcefPosition'. +instance ETransform (AngularPosition NVector) where + fromEcef = ecefToNVector + toEcef = nvectorToEcef + +-- | 'AngularPosition' of 'LatLong' <-> 'EcefPosition'. +instance ETransform (AngularPosition LatLong) where + fromEcef p e = fromNVector (ecefToNVector p e) + toEcef = nvectorToEcef . toNVector + +-- | 'ETransform' identity. +instance ETransform EcefPosition where + fromEcef p _ = p + toEcef p _ = p + +-- | @nvectorToLatLong v@ transforms 'NVector' @v@ to an equivalent 'LatLong'. +-- +-- See also 'toNVector'. +nvectorToLatLong :: NVector -> LatLong +nvectorToLatLong nv = latLong lat lon + where + v = vec nv + lat = atan2' (vz v) (sqrt (vx v * vx v + vy v * vy v)) + lon = atan2' (vy v) (vx v) + +-- | @latLongToNVector ll@ transforms 'LatLong' @ll@ to an equivalent 'NVector'. +-- +-- See also 'fromNVector'. +latLongToNVector :: LatLong -> NVector +latLongToNVector ll = nvector x' y' z' + where + lat = latitude ll + lon = longitude ll + cl = cos' lat + x' = cl * cos' lon + y' = cl * sin' lon + z' = sin' lat + +-- | @ecefToNVector p e@ transforms 'EcefPosition' @p@ to an equivalent 'NVector' and geodetic height +-- using earth model @e@. +-- +-- See also 'fromEcef' +ecefToNVector :: EcefPosition -> Earth -> AngularPosition NVector +-- Ellipsoidal +ecefToNVector ep e@(Ellipsoidal el) = nvectorHeight (nvecEllipsoidal d e2 k px py pz) (metres h) + where + ev = vec ep + e' = eccentricity e + e2 = e' * e' + e4 = e2 * e2 + a = toMetres (equatorialRadius el) + a2 = a * a + px = vx ev + py = vy ev + pz = vz ev + p = (px * px + py * py) / a2 + q = ((1 - e2) / a2) * (pz * pz) + r = (p + q - e4) / 6.0 + s = (e4 * p * q) / (4.0 * r * r * r) + t = (1.0 + s + sqrt (s * (2.0 + s))) ** (1 / 3) + u = r * (1.0 + t + 1.0 / t) + v = sqrt (u * u + q * e4) + w = e2 * (u + v - q) / (2.0 * v) + k = sqrt (u + v + w * w) - w + d = k * sqrt (px * px + py * py) / (k + e2) + h = ((k + e2 - 1.0) / k) * sqrt (d * d + pz * pz) +-- Spherical +ecefToNVector p (Spherical r) = nvectorHeight (nvector (vx nv) (vy nv) (vz nv)) h + where + ev = vec p + nv = vunit ev + h = sub (metres (vnorm ev)) r + +nvecEllipsoidal :: Double -> Double -> Double -> Double -> Double -> Double -> NVector +nvecEllipsoidal d e2 k px py pz = nvector nx' ny' nz' + where + s = 1.0 / sqrt (d * d + pz * pz) + a = k / (k + e2) + nx' = s * a * px + ny' = s * a * py + nz' = s * pz + +-- | @nvectorToEcef (n, h) e@ transforms 'NVector' @n@ and geodetic height @h@ +-- to an equivalent 'EcefPosition' using earth model @e@. +-- +-- See also 'toEcef' +nvectorToEcef :: AngularPosition NVector -> Earth -> EcefPosition +-- Ellipsoidal +nvectorToEcef (AngularPosition nv h) e@(Ellipsoidal el) = ecef ex' ey' ez' + where + v = vec nv + uv = vunit v + a = toMetres (equatorialRadius el) + b = toMetres (polarRadius e) + nx' = vx uv + ny' = vy uv + nz' = vz uv + m = (a * a) / (b * b) + n = b / sqrt ((nx' * nx' * m) + (ny' * ny' * m) + (nz' * nz')) + h' = toMetres h + ex' = metres (n * m * nx' + h' * nx') + ey' = metres (n * m * ny' + h' * ny') + ez' = metres (n * nz' + h' * nz') +-- Spherical +nvectorToEcef (AngularPosition nv h) (Spherical r) = ecefMetres (vx ev) (vy ev) (vz ev) + where + unv = vunit . vec $ nv + n = add h r + ev = vscale unv (toMetres n) + +-- | @geodeticHeight p e@ computes the geodetic height of 'EcefPosition' @p@ using earth model @e@. +-- +-- The geodetic height (or ellipsoidal height) is __not__ the mean sea level (MSL) height. +geodeticHeight :: EcefPosition -> Earth -> Length +geodeticHeight p e = height (ecefToNVector p e)
src/Data/Geo/Jord/Vector3d.hs view
@@ -99,7 +99,7 @@ vzero :: Vector3d vzero = Vector3d 0 0 0 --- | transpose __square__ matrix made of 'Vector3d'. +-- | transpose __square (3x3)__ matrix of 'Vector3d'. transpose :: [Vector3d] -> [Vector3d] transpose m = fmap ds2v (transpose' xs) where @@ -110,7 +110,7 @@ transpose' ([]:_) = [] transpose' x = map head x : transpose' (map tail x) --- | multiplies 2 matrices of 'Vector3d'. +-- | multiplies 2 __square (3x3)__ matrices of 'Vector3d'. mdot :: [Vector3d] -> [Vector3d] -> [Vector3d] mdot a b = fmap ds2v [[vdot ar bc | bc <- transpose b] | ar <- a]
test/Data/Geo/Jord/AngleSpec.hs view
@@ -1,92 +1,92 @@-module Data.Geo.Jord.AngleSpec- ( spec- ) where--import Data.Geo.Jord-import System.IO-import Test.Hspec--spec :: Spec-spec = do- describe "Reading valid angles" $ do- it "reads 55°36'21\"" $ readAngle "55°36'21\"" `shouldBe` decimalDegrees 55.6058333- it "reads 55°36'21''" $ readAngle "55°36'21''" `shouldBe` decimalDegrees 55.6058333- it "reads 55d36m21.0s" $ readAngle "55d36m21.0s" `shouldBe` decimalDegrees 55.6058333- it "reads 55.6058333°" $ readAngle "55.6058333°" `shouldBe` decimalDegrees 55.6058333- it "reads -55.6058333°" $ readAngle "-55.6058333°" `shouldBe` decimalDegrees (-55.6058333)- it "reads 96°01′18″" $ do- hSetEncoding stdin utf8- hSetEncoding stdout utf8- hSetEncoding stderr utf8- readAngle "96°01′18″" `shouldBe` decimalDegrees 96.02166666- describe "Adding/Subtracting angles" $ do- it "adds angles" $- add (decimalDegrees 55.6058333) (decimalDegrees 5.0) `shouldBe`- decimalDegrees 60.6058333- it "subtracts angles" $- sub (decimalDegrees 5.0) (decimalDegrees 55.6058333) `shouldBe`- decimalDegrees (-50.6058333)- describe "Angle normalisation" $ do- it "370 degrees normalised to [0..360] = 10" $- normalise (decimalDegrees 370) (decimalDegrees 360) `shouldBe` decimalDegrees 10- it "350 degrees normalised to [0..360] = 350" $- normalise (decimalDegrees 350) (decimalDegrees 360) `shouldBe` decimalDegrees 350- describe "Angle equality" $ do- it "considers 59.9999999° == 60.0°" $ decimalDegrees 59.9999999 `shouldBe` decimalDegrees 60- it "considers 59.9999998° /= 60.0°" $- decimalDegrees 59.9999998 `shouldNotBe` decimalDegrees 60- describe "Showing angles" $ do- it "shows 59.99999999999999 as 60°0'0.000\"" $- show (decimalDegrees 59.99999999999999) `shouldBe` "60°0'0.000\""- it "shows 154.915 as 154°54'54.000\"" $- show (decimalDegrees 154.915) `shouldBe` "154°54'54.000\""- it "shows -154.915 as -154°54'54.000\"" $- show (decimalDegrees (-154.915)) `shouldBe` "-154°54'54.000\""- it "show 0.5245 as 0°31'28.800\"" $ show (decimalDegrees 0.5245) `shouldBe` "0°31'28.200\""- it "show -0.5245 as -0°31'28.800\"" $- show (decimalDegrees (-0.5245)) `shouldBe` "-0°31'28.200\""- describe "Angle from decimal degrees" $ do- it "returns 1 millisecond when called with 1 / 3600000" $ do- let actual = decimalDegrees (1 / 3600000)- getDegrees actual `shouldBe` 0- getMinutes actual `shouldBe` 0- getSeconds actual `shouldBe` 0- getMilliseconds actual `shouldBe` 1- it "returns 1 second when called with 1000 / 3600000" $ do- let actual = decimalDegrees (1000 / 3600000)- getDegrees actual `shouldBe` 0- getMinutes actual `shouldBe` 0- getSeconds actual `shouldBe` 1- getMilliseconds actual `shouldBe` 0- it "returns 1 minute when called with 60000 / 3600000" $ do- let actual = decimalDegrees (60000 / 3600000)- getDegrees actual `shouldBe` 0- getMinutes actual `shouldBe` 1- getSeconds actual `shouldBe` 0- getMilliseconds actual `shouldBe` 0- it "returns 1 degree when called with 1" $ do- let actual = decimalDegrees 1- getDegrees actual `shouldBe` 1- getMinutes actual `shouldBe` 0- getSeconds actual `shouldBe` 0- getMilliseconds actual `shouldBe` 0- it "accepts positve values" $ do- let actual = decimalDegrees 154.9150300- getDegrees actual `shouldBe` 154- getMinutes actual `shouldBe` 54- getSeconds actual `shouldBe` 54- getMilliseconds actual `shouldBe` 108- it "accepts negative values" $ do- let actual = decimalDegrees (-154.915)- getDegrees actual `shouldBe` (-154)- getMinutes actual `shouldBe` 54- getSeconds actual `shouldBe` 54- getMilliseconds actual `shouldBe` 0- describe "Arc length" $ do- it "computes the length of an arc with a central angle of 1 milliseconds" $- arcLength (decimalDegrees (1.0 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0.031- it- "arc length with central angle of 0.6 milliseconds == arc length with central angle of 1 milliseconds" $- arcLength (decimalDegrees (0.6 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0.031- it "arc length with central angle of 0.5 milliseconds == 0" $- arcLength (decimalDegrees (0.4 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0+module Data.Geo.Jord.AngleSpec + ( spec + ) where + +import Data.Geo.Jord +import System.IO +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid angles" $ do + it "reads 55°36'21\"" $ readAngle "55°36'21\"" `shouldBe` decimalDegrees 55.6058333 + it "reads 55°36'21''" $ readAngle "55°36'21''" `shouldBe` decimalDegrees 55.6058333 + it "reads 55d36m21.0s" $ readAngle "55d36m21.0s" `shouldBe` decimalDegrees 55.6058333 + it "reads 55.6058333°" $ readAngle "55.6058333°" `shouldBe` decimalDegrees 55.6058333 + it "reads -55.6058333°" $ readAngle "-55.6058333°" `shouldBe` decimalDegrees (-55.6058333) + it "reads 96°01′18″" $ do + hSetEncoding stdin utf8 + hSetEncoding stdout utf8 + hSetEncoding stderr utf8 + readAngle "96°01′18″" `shouldBe` decimalDegrees 96.02166666 + describe "Adding/Subtracting angles" $ do + it "adds angles" $ + add (decimalDegrees 55.6058333) (decimalDegrees 5.0) `shouldBe` + decimalDegrees 60.6058333 + it "subtracts angles" $ + sub (decimalDegrees 5.0) (decimalDegrees 55.6058333) `shouldBe` + decimalDegrees (-50.6058333) + describe "Angle normalisation" $ do + it "370 degrees normalised to [0..360] = 10" $ + normalise (decimalDegrees 370) (decimalDegrees 360) `shouldBe` decimalDegrees 10 + it "350 degrees normalised to [0..360] = 350" $ + normalise (decimalDegrees 350) (decimalDegrees 360) `shouldBe` decimalDegrees 350 + describe "Angle equality" $ do + it "considers 59.9999999° == 60.0°" $ decimalDegrees 59.9999999 `shouldBe` decimalDegrees 60 + it "considers 59.9999998° /= 60.0°" $ + decimalDegrees 59.9999998 `shouldNotBe` decimalDegrees 60 + describe "Showing angles" $ do + it "shows 59.99999999999999 as 60°0'0.000\"" $ + show (decimalDegrees 59.99999999999999) `shouldBe` "60°0'0.000\"" + it "shows 154.915 as 154°54'54.000\"" $ + show (decimalDegrees 154.915) `shouldBe` "154°54'54.000\"" + it "shows -154.915 as -154°54'54.000\"" $ + show (decimalDegrees (-154.915)) `shouldBe` "-154°54'54.000\"" + it "show 0.5245 as 0°31'28.800\"" $ show (decimalDegrees 0.5245) `shouldBe` "0°31'28.200\"" + it "show -0.5245 as -0°31'28.800\"" $ + show (decimalDegrees (-0.5245)) `shouldBe` "-0°31'28.200\"" + describe "Angle from decimal degrees" $ do + it "returns 1 millisecond when called with 1 / 3600000" $ do + let actual = decimalDegrees (1 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 1 + it "returns 1 second when called with 1000 / 3600000" $ do + let actual = decimalDegrees (1000 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 1 + getMilliseconds actual `shouldBe` 0 + it "returns 1 minute when called with 60000 / 3600000" $ do + let actual = decimalDegrees (60000 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 1 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 0 + it "returns 1 degree when called with 1" $ do + let actual = decimalDegrees 1 + getDegrees actual `shouldBe` 1 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 0 + it "accepts positve values" $ do + let actual = decimalDegrees 154.9150300 + getDegrees actual `shouldBe` 154 + getMinutes actual `shouldBe` 54 + getSeconds actual `shouldBe` 54 + getMilliseconds actual `shouldBe` 108 + it "accepts negative values" $ do + let actual = decimalDegrees (-154.915) + getDegrees actual `shouldBe` (-154) + getMinutes actual `shouldBe` 54 + getSeconds actual `shouldBe` 54 + getMilliseconds actual `shouldBe` 0 + describe "Arc length" $ do + it "computes the length of an arc with a central angle of 1 milliseconds" $ + arcLength (decimalDegrees (1.0 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0.0309 + it + "arc length with central angle of 0.6 milliseconds == arc length with central angle of 1 milliseconds" $ + arcLength (decimalDegrees (0.6 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0.0309 + it "arc length with central angle of 0.5 milliseconds == 0" $ + arcLength (decimalDegrees (0.4 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0
+ test/Data/Geo/Jord/DurationSpec.hs view
@@ -0,0 +1,26 @@+module Data.Geo.Jord.DurationSpec+ ( spec+ ) where++import Data.Geo.Jord+import Test.Hspec++spec :: Spec+spec = do+ describe "Reading valid durations" $ do+ it "reads 1H45M36.5S" $ readDuration "1H45M36.5S" `shouldBe` hms 1 45 36.5+ it "reads 45M" $ readDuration "45M" `shouldBe` minutes 45+ it "reads 36S" $ readDuration "36S" `shouldBe` seconds 36+ it "reads 36.6S" $ readDuration "36.6S" `shouldBe` milliseconds 36600+ it "read 1H-30M" $ readDuration "1H-30M" `shouldBe` hours 0.5+ describe "Reading invalid duration" $+ it "fails to read 5" $ readDurationE "5" `shouldBe` Left "couldn't read duration 5"+ describe "Showing duration" $+ it "shows duration" $ show (hms 1 45 36.5) `shouldBe` "1H45M36.500S"+ describe "Converting duration" $ do+ it "converts hours to seconds" $ toSeconds (hours 1) `shouldBe` 3600.0+ it "converts minutes to hours" $ toHours (minutes 30) `shouldBe` 0.5+ it "converts duration to milliseconds" $ toMilliseconds (hms 1 54 3.154) `shouldBe` 6843154+ describe "Adding/Subtracting duration" $ do+ it "adds duration" $ add (minutes 45) (seconds 36) `shouldBe` hms 0 45 36+ it "subtracts duration" $ sub (hours 1) (minutes 60) `shouldBe` zero
test/Data/Geo/Jord/EarthSpec.hs view
@@ -1,27 +1,27 @@-module Data.Geo.Jord.EarthSpec- ( spec- ) where--import Data.Geo.Jord-import Test.Hspec--spec :: Spec-spec = do- describe "Eccentricity" $ do- it "returns 0.08181919084262157 for the WGS84 ellipsoid" $- eccentricity wgs84 `shouldBe` 0.08181919084262157- it "returns 0.08181919104281514 for the GRS80 ellipsoid" $- eccentricity grs80 `shouldBe` 0.08181919104281514- it "returns 0.08181881066274845 for the WG72 ellipsoid" $- eccentricity wgs72 `shouldBe` 0.08181881066274845- describe "Polar radius" $ do- it "returns 6356752.314 m for the WGS84 ellipsoid" $- polarRadius wgs84 `shouldBe` metres 6356752.314- it "returns 6356752.314 m for the GRS80 ellipsoid" $- polarRadius grs80 `shouldBe` metres 6356752.314- it "returns 6356750.52 m for the WG72 ellipsoid" $- polarRadius wgs72 `shouldBe` metres 6356750.52- describe "Mean radius" $ do- it "returns 6371008.771 m for the WGS84 ellipsoid" $ r84 `shouldBe` metres 6371008.771- it "returns 6371008.771 m for the GRS80 ellipsoid" $ r80 `shouldBe` metres 6371008.771- it "returns 6371006.84 m for the WG72 ellipsoid" $ r72 `shouldBe` metres 6371006.84+module Data.Geo.Jord.EarthSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Eccentricity" $ do + it "returns 0.08181919084262157 for the WGS84 ellipsoid" $ + eccentricity wgs84 `shouldBe` 0.08181919084262157 + it "returns 0.08181919104281514 for the GRS80 ellipsoid" $ + eccentricity grs80 `shouldBe` 0.08181919104281514 + it "returns 0.08181881066274845 for the WG72 ellipsoid" $ + eccentricity wgs72 `shouldBe` 0.08181881066274845 + describe "Polar radius" $ do + it "returns 6356752.3142 m for the WGS84 ellipsoid" $ + polarRadius wgs84 `shouldBe` metres 6356752.3142 + it "returns 6356752.3141 m for the GRS80 ellipsoid" $ + polarRadius grs80 `shouldBe` metres 6356752.3141 + it "returns 6356750.52 m for the WG72 ellipsoid" $ + polarRadius wgs72 `shouldBe` metres 6356750.52 + describe "Mean radius" $ do + it "returns 6371008.7714 m for the WGS84 ellipsoid" $ r84 `shouldBe` metres 6371008.7714 + it "returns 6371008.7714 m for the GRS80 ellipsoid" $ r80 `shouldBe` metres 6371008.7714 + it "returns 6371006.84 m for the WG72 ellipsoid" $ r72 `shouldBe` metres 6371006.84
test/Data/Geo/Jord/FramesSpec.hs view
@@ -16,7 +16,7 @@ it "computes the target point from p0 and NED" $ do let p0 = decimalLatLong 49.66618 3.45063 let d = nedMetres (-86126) (-78900) 1069 - targetN p0 d wgs84 `shouldBe` decimalLatLong 48.8866688 2.374721111 + targetN p0 d wgs84 `shouldBe` decimalLatLong 48.8866688 2.374721388 it "computes the target point from p0 and vector in Frame B" $ do let p0 = decimalLatLongHeight 49.66618 3.45063 zero let y = decimalDegrees 10 -- yaw @@ -24,25 +24,25 @@ let p = decimalDegrees 30 -- pitch let d = deltaMetres 3000 2000 100 target p0 (frameB y r p) d wgs84 `shouldBe` - decimalLatLongHeight 49.6918016 3.4812669 (metres 6.007) + decimalLatLongHeight 49.6918016 3.4812669 (metres 6.0077) describe "nedBetween" $ do it "computes NED between LatLong positions" $ do let p1 = decimalLatLong 49.66618 3.45063 let p2 = decimalLatLong 48.88667 2.37472 let d = nedBetween p1 p2 wgs84 - d `shouldBe` nedMetres (-86125.88049540376) (-78900.08718759022) 1069.1981930266265 + d `shouldBe` nedMetres (-86125.8805) (-78900.0878) 1069.1984 it "computes NED between angular positions" $ do let p1 = decimalLatLongHeight 49.66618 3.45063 zero let p2 = decimalLatLongHeight 48.88667 2.37472 zero let d = nedBetween p1 p2 wgs84 - d `shouldBe` nedMetres (-86125.88049540376) (-78900.08718759022) 1069.1981930266265 + d `shouldBe` nedMetres (-86125.8805) (-78900.0878) 1069.1984 describe "deltaBetween" $ it "computes delta between angular positions in frame L" $ do let p1 = decimalLatLongHeight 1 2 (metres (-3)) let p2 = decimalLatLongHeight 4 5 (metres (-6)) let w = decimalDegrees 5 -- wander azimuth let d = deltaBetween p1 p2 (frameL w) wgs84 - d `shouldBe` deltaMetres 359490.579 302818.523 17404.272 + d `shouldBe` deltaMetres 359490.5782 302818.5226 17404.2713 describe "deltaBetween and target consistency" $ it "computes targetN p1 (nedBetween p1 p2) = p2" $ do let p1 = decimalLatLongHeight 49.66618 3.45063 zero @@ -53,22 +53,22 @@ let p = decimalLatLong 49.66618 3.45063 let f = frameN p wgs84 rEF f `shouldBe` - [ Vector3d (-0.7609044147490918) (-6.018845508229421e-2) (-0.6460664218872152) - , Vector3d (-4.588084170935741e-2) 0.9981870315100305 (-3.8956366478847045e-2) - , Vector3d 0.6472398473358879 0.0 (-0.7622864160016345) + [ Vector3d (-0.7609044147650337) (-6.0188455103478165e-2) (-0.6460664218664659) + , Vector3d (-4.588084172652564e-2) 0.9981870315087531 (-3.8956366491356864e-2) + , Vector3d 0.6472398473159291 0.0 (-0.7622864160185809) ] it "computes the rotation matrix to go from Frame B to earth-fixed frame" $ do let p = decimalLatLong 49.66618 3.45063 let f = frameB (decimalDegrees 10) (decimalDegrees 20) (decimalDegrees 30) p wgs84 rEF f `shouldBe` - [ Vector3d (-0.4930071357732754) (-0.3703899170519431) (-0.7872453705312508) - , Vector3d 0.13374504887910177 0.8618333991702691 (-0.48923966925725315) - , Vector3d 0.8596837941807199 (-0.34648881860873165) (-0.3753521980782417) + [ Vector3d (-0.4930071357985816) (-0.3703899170611777) (-0.787245370511058) + , Vector3d 0.13374504886728417 0.8618333991629544 (-0.4892396692733688) + , Vector3d 0.8596837941680457 (-0.3464888186170537) (-0.37535219809958753) ] describe "North, East, Down delta" $ do - describe "norm" $ - it "computes the norm of a NED vector" $ - norm (nedMetres (-86126) (-78900) 1069) `shouldBe` metres 116807.708 + describe "slantRange" $ + it "computes the slant range of a NED vector" $ + slantRange (nedMetres (-86126) (-78900) 1069) `shouldBe` metres 116807.708 describe "bearing" $ it "computes the bearing of a NED vector" $ bearing (nedMetres (-86126) (-78900) 1069) `shouldBe` decimalDegrees 222.4927888
test/Data/Geo/Jord/GeodeticsSpec.hs view
@@ -1,213 +1,213 @@-module Data.Geo.Jord.GeodeticsSpec- ( spec- ) where--import Control.Exception.Base-import Data.Geo.Jord-import Data.Maybe (fromJust)-import Test.Hspec--spec :: Spec-spec = do- describe "antipode" $ do- it "returns the antipodal point" $ do- let p = latLongHeight (readLatLong "484137N0061105E") (metres 15000)- let e = decimalLatLongHeight (-48.6936111) (-173.8152777) (metres 15000)- antipode p `shouldBe` e- it "returns the south pole when called with the north pole" $- antipode northPole `shouldBe` southPole- it "returns the north pole when called with the south pole" $- antipode southPole `shouldBe` northPole- describe "crossTrackDistance" $ do- it "returns a negative length when position is left of great circle (bearing)" $ do- let p = decimalLatLong 53.2611 (-0.7972)- let gc = greatCircleBearing (decimalLatLong 53.3206 (-1.7297)) (decimalDegrees 96.0)- crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres (-305.663)- it "returns a negative length when position is left of great circle" $ do- let p = decimalLatLong 53.2611 (-0.7972)- let gc = greatCircle (decimalLatLong 53.3206 (-1.7297)) (decimalLatLong 53.1887 0.1334)- crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres (-307.547)- it "returns a positve length when position is right of great circle (bearing)" $ do- let p = readLatLong "531540N0014750W"- let gc = greatCircleBearing (readLatLong "531914N0014347W") (readAngle "96d01m18s")- crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres 7042.324- it "returns a positive length when position is left of great circle" $ do- let p = antipode (decimalLatLong 53.2611 (-0.7972))- let gc = greatCircle (decimalLatLong 53.3206 (-1.7297)) (decimalLatLong 53.1887 0.1334)- crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres 307.547- describe "destination" $ do- it "return the given point if distance is 0 meter" $ do- let p0 = readLatLong "531914N0014347W"- destination p0 (decimalDegrees 96.0217) zero r84 `shouldBe` p0- it "return the angular position along great-circle at distance and bearing" $ do- let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000.0)- let p1 = decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0)- destination p0 (decimalDegrees 96.0217) (metres 124800) r84 `shouldBe` p1- it "return the ECEF position along great-circle at distance and bearing" $ do- let p0 = ecefToNVector (ecefMetres 3812864.094 (-115142.863) 5121515.161) s84- let p1 = ecefMetres 3826406.4710518294 8900.536398998282 5112694.233184049- let p = destination84 p0 (decimalDegrees 96.0217) (metres 124800)- nvectorToEcef p s84 `shouldBe` p1- describe "finalBearing" $ do- it "returns the Nothing if both point are the same" $ do- let p = readLatLong "500359N0054253W"- finalBearing p p `shouldBe` Nothing- it "returns 0° if both point have the same longitude (going north)" $ do- let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)- let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000)- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 0)- it "returns 180° if both point have the same longitude (going south)" $ do- let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000)- let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 180)- it "returns 90° at the equator going east" $ do- let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000)- let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000)- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 90)- it "returns 270° at the equator going west" $ do- let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000)- let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000)- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 270)- it "returns the final bearing in compass angle" $ do- let p1 = readLatLong "500359N0054253W"- let p2 = readLatLong "583838N0030412W"- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 11.2752013)- it "returns the final bearing in compass angle" $ do- let p1 = readLatLong "583838N0030412W"- let p2 = readLatLong "500359N0054253W"- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 189.1198181)- it "returns the final bearing in compass angle" $ do- let p1 = readLatLong "535941S0255915W"- let p2 = readLatLong "54N154E"- finalBearing p1 p2 `shouldBe` Just (decimalDegrees 125.6839436)- describe "greatCircle smart constructors" $ do- it "fails if both positions are equal" $- greatCircleE (decimalLatLong 3 154) (decimalLatLong 3 154) `shouldBe`- Left "Invalid Great Circle: positions are equal"- it "fails if both positions are antipodal" $- greatCircleE (decimalLatLong 3 154) (antipode (decimalLatLong 3 154)) `shouldBe`- Left "Invalid Great Circle: positions are antipodal"- describe "initialBearing" $ do- it "returns Nothing if both point are the same" $ do- let p = readLatLong "500359N1795959W"- initialBearing p p `shouldBe` Nothing- it "returns 0° if both point have the same longitude (going north)" $ do- let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)- let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 0)- it "returns 180° if both point have the same longitude (going south)" $ do- let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000)- let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 180)- it "returns 90° at the equator going east" $ do- let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000)- let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 90)- it "returns 270° at the equator going west" $ do- let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000)- let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 270)- it "returns the initial bearing in compass angle" $ do- let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)- let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 9.1198181)- it "returns the initial bearing in compass angle" $ do- let p1 = latLongHeight (readLatLong "583838N0030412W") (metres 12000)- let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)- initialBearing p1 p2 `shouldBe` Just (decimalDegrees 191.2752013)- describe "interpolate" $ do- let p1 = readLatLong "44N044E"- let p2 = readLatLong "46N046E"- it "fails if f < 0.0" $- evaluate (interpolate p1 p2 (-0.5)) `shouldThrow`- errorCall "fraction must be in range [0..1], was -0.5"- it "fails if f > 1.0" $- evaluate (interpolate p1 p2 1.1) `shouldThrow`- errorCall "fraction must be in range [0..1], was 1.1"- it "returns p0 if f == 0" $ interpolate p1 p2 0.0 `shouldBe` p1- it "returns p1 if f == 1" $ interpolate p1 p2 1.0 `shouldBe` p2- it "returns the interpolated position" $ do- let p3 = latLongHeight (readLatLong "53°28'46''N 2°14'43''W") (metres 10000)- let p4 = latLongHeight (readLatLong "55°36'21''N 13°02'09''E") (metres 20000)- interpolate p3 p4 0.5 `shouldBe`- decimalLatLongHeight 54.7835574 5.1949856 (metres 15000)- describe "insideSurface" $ do- let p1 = decimalLatLong 45 1- let p2 = decimalLatLong 45 2- let p3 = decimalLatLong 46 1- let p4 = decimalLatLong 46 2- let p5 = decimalLatLong 45.1 1.1- it "return False if polygon is empty" $ insideSurface p1 [] `shouldBe` False- it "return False if polygon does not define at least a triangle" $- insideSurface p1 [p1, p2] `shouldBe` False- it "returns True if point is inside polygon" $ do- let polygon = [p1, p2, p4, p3]- insideSurface p5 polygon `shouldBe` True- it "returns False if point is inside polygon" $ do- let polygon = [p1, p2, p4, p3]- let p = antipode p5- insideSurface p polygon `shouldBe` False- it "returns False if point is a vertex of the polygon" $ do- let polygon = [p1, p2, p4, p3]- insideSurface p1 polygon `shouldBe` False- it "handles closed polygons" $ do- let polygon = [p1, p2, p4, p3, p1]- insideSurface p5 polygon `shouldBe` True- it "handles concave polygons" $ do- let malmo = decimalLatLong 55.6050 13.0038- let ystad = decimalLatLong 55.4295 13.82- let lund = decimalLatLong 55.7047 13.1910- let helsingborg = decimalLatLong 56.0465 12.6945- let kristianstad = decimalLatLong 56.0294 14.1567- let polygon = [malmo, ystad, kristianstad, helsingborg, lund]- let hoor = decimalLatLong 55.9295 13.5297- let hassleholm = decimalLatLong 56.1589 13.7668- insideSurface hoor polygon `shouldBe` True- insideSurface hassleholm polygon `shouldBe` False- describe "intersections" $ do- it "returns nothing if both great circle are equals" $ do- let gc = greatCircleBearing (decimalLatLong 51.885 0.235) (decimalDegrees 108.63)- (intersections gc gc :: Maybe (LatLong, LatLong)) `shouldBe` Nothing- it "returns nothing if both great circle are equals (opposite orientation)" $ do- let gc1 = greatCircle (decimalLatLong 51.885 0.235) (decimalLatLong 52.885 1.235)- let gc2 = greatCircle (decimalLatLong 52.885 1.235) (decimalLatLong 51.885 0.235)- (intersections gc1 gc2 :: Maybe (LatLong, LatLong)) `shouldBe` Nothing- it "returns the two points where the two great circles intersects" $ do- let gc1 = greatCircleBearing (decimalLatLong 51.885 0.235) (decimalDegrees 108.63)- let gc2 = greatCircleBearing (decimalLatLong 49.008 2.549) (decimalDegrees 32.72)- let (i1, i2) = fromJust (intersections gc1 gc2)- i1 `shouldBe` decimalLatLong 50.9017226 4.4942782- i2 `shouldBe` antipode i1- describe "mean" $ do- it "returns Nothing if no point is given" $ (mean [] :: Maybe NVector) `shouldBe` Nothing- it "returns the unique given point" $ do- let p = readLatLong "500359N0054253W"- mean [p] `shouldBe` Just p- it "returns the geographical mean" $ do- let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 15000.0)- let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 25000.0)- let e = decimalLatLongHeight 54.3622869 (-4.5306725) zero- mean [p1, p2] `shouldBe` Just e- it "returns Nothing if list contains antipodal points" $ do- let points =- [ decimalLatLong 45 1- , decimalLatLong 45 2- , decimalLatLong 46 2- , decimalLatLong 46 1- , antipode (decimalLatLong 45 2)- ]- mean points `shouldBe` Nothing- describe "surfaceDistance" $ do- it "returns 0 if both points are equal" $ do- let p = readLatLong "500359N1795959W"- surfaceDistance p p r84 `shouldBe` zero- it "returns the distance between 2 points" $ do- let p1 = readLatLong "500359N0054253W"- let p2 = readLatLong "583838N0030412W"- surfaceDistance84 p1 p2 `shouldBe` metres 968854.868- it "handles singularity at the pole" $- surfaceDistance northPole southPole r84 `shouldBe` kilometres 20015.114351- it "handles the discontinuity at the Date Line" $ do- let p1 = readLatLong "500359N1795959W"- let p2 = readLatLong "500359N1795959E"- surfaceDistance p1 p2 (meanRadius wgs84) `shouldBe` metres 39.66+module Data.Geo.Jord.GeodeticsSpec + ( spec + ) where + +import Control.Exception.Base +import Data.Geo.Jord +import Data.Maybe (fromJust) +import Test.Hspec + +spec :: Spec +spec = do + describe "antipode" $ do + it "returns the antipodal point" $ do + let p = latLongHeight (readLatLong "484137N0061105E") (metres 15000) + let e = decimalLatLongHeight (-48.6936111) (-173.8152777) (metres 15000) + antipode p `shouldBe` e + it "returns the south pole when called with the north pole" $ + antipode northPole `shouldBe` southPole + it "returns the north pole when called with the south pole" $ + antipode southPole `shouldBe` northPole + describe "crossTrackDistance" $ do + it "returns a negative length when position is left of great circle (bearing)" $ do + let p = decimalLatLong 53.2611 (-0.7972) + let gc = greatCircle (decimalLatLong 53.3206 (-1.7297), decimalDegrees 96.0) + crossTrackDistance p gc r84 `shouldBe` metres (-305.6629) + it "returns a negative length when position is left of great circle" $ do + let p = decimalLatLong 53.2611 (-0.7972) + let gc = greatCircle (decimalLatLong 53.3206 (-1.7297), decimalLatLong 53.1887 0.1334) + crossTrackDistance p gc r84 `shouldBe` metres (-307.5471) + it "returns a positve length when position is right of great circle (bearing)" $ do + let p = readLatLong "531540N0014750W" + let gc = greatCircle (readLatLong "531914N0014347W", readAngle "96d01m18s") + crossTrackDistance p gc r84 `shouldBe` metres 7042.3242 + it "returns a positive length when position is left of great circle" $ do + let p = antipode (decimalLatLong 53.2611 (-0.7972)) + let gc = greatCircle (decimalLatLong 53.3206 (-1.7297), decimalLatLong 53.1887 0.1334) + crossTrackDistance p gc r84 `shouldBe` metres 307.5471 + describe "destination" $ do + it "return the given point if distance is 0 meter" $ do + let p0 = readLatLong "531914N0014347W" + destination p0 (decimalDegrees 96.0217) zero r84 `shouldBe` p0 + it "return the angular position along great-circle at distance and bearing" $ do + let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000.0) + let p1 = decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0) + destination p0 (decimalDegrees 96.0217) (metres 124800) r84 `shouldBe` p1 + it "return the ECEF position along great-circle at distance and bearing" $ do + let p0 = ecefToNVector (ecefMetres 3812864.094 (-115142.863) 5121515.161) s84 + let p1 = ecefMetres 3826406.471 8900.5364 5112694.2331 + let p = destination84 p0 (decimalDegrees 96.0217) (metres 124800) + nvectorToEcef p s84 `shouldBe` p1 + describe "finalBearing" $ do + it "returns the Nothing if both point are the same" $ do + let p = readLatLong "500359N0054253W" + finalBearing p p `shouldBe` Nothing + it "returns 0° if both point have the same longitude (going north)" $ do + let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000) + let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000) + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 0) + it "returns 180° if both point have the same longitude (going south)" $ do + let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000) + let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000) + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 180) + it "returns 90° at the equator going east" $ do + let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000) + let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000) + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 90) + it "returns 270° at the equator going west" $ do + let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000) + let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000) + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 270) + it "returns the final bearing in compass angle" $ do + let p1 = readLatLong "500359N0054253W" + let p2 = readLatLong "583838N0030412W" + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 11.2752013) + it "returns the final bearing in compass angle" $ do + let p1 = readLatLong "583838N0030412W" + let p2 = readLatLong "500359N0054253W" + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 189.1198181) + it "returns the final bearing in compass angle" $ do + let p1 = readLatLong "535941S0255915W" + let p2 = readLatLong "54N154E" + finalBearing p1 p2 `shouldBe` Just (decimalDegrees 125.6839436) + describe "greatCircle smart constructors" $ do + it "fails if both positions are equal" $ + greatCircleE (decimalLatLong 3 154, decimalLatLong 3 154) `shouldBe` + Left "Invalid Great Circle: positions are equal" + it "fails if both positions are antipodal" $ + greatCircleE (decimalLatLong 3 154, antipode (decimalLatLong 3 154)) `shouldBe` + Left "Invalid Great Circle: positions are antipodal" + describe "initialBearing" $ do + it "returns Nothing if both point are the same" $ do + let p = readLatLong "500359N1795959W" + initialBearing p p `shouldBe` Nothing + it "returns 0° if both point have the same longitude (going north)" $ do + let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000) + let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 0) + it "returns 180° if both point have the same longitude (going south)" $ do + let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000) + let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 180) + it "returns 90° at the equator going east" $ do + let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000) + let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 90) + it "returns 270° at the equator going west" $ do + let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000) + let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 270) + it "returns the initial bearing in compass angle" $ do + let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000) + let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 9.1198181) + it "returns the initial bearing in compass angle" $ do + let p1 = latLongHeight (readLatLong "583838N0030412W") (metres 12000) + let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000) + initialBearing p1 p2 `shouldBe` Just (decimalDegrees 191.2752013) + describe "interpolate" $ do + let p1 = readLatLong "44N044E" + let p2 = readLatLong "46N046E" + it "fails if f < 0.0" $ + evaluate (interpolate p1 p2 (-0.5)) `shouldThrow` + errorCall "fraction must be in range [0..1], was -0.5" + it "fails if f > 1.0" $ + evaluate (interpolate p1 p2 1.1) `shouldThrow` + errorCall "fraction must be in range [0..1], was 1.1" + it "returns p0 if f == 0" $ interpolate p1 p2 0.0 `shouldBe` p1 + it "returns p1 if f == 1" $ interpolate p1 p2 1.0 `shouldBe` p2 + it "returns the interpolated position" $ do + let p3 = latLongHeight (readLatLong "53°28'46''N 2°14'43''W") (metres 10000) + let p4 = latLongHeight (readLatLong "55°36'21''N 13°02'09''E") (metres 20000) + interpolate p3 p4 0.5 `shouldBe` + decimalLatLongHeight 54.7835574 5.1949856 (metres 15000) + describe "insideSurface" $ do + let p1 = decimalLatLong 45 1 + let p2 = decimalLatLong 45 2 + let p3 = decimalLatLong 46 1 + let p4 = decimalLatLong 46 2 + let p5 = decimalLatLong 45.1 1.1 + it "return False if polygon is empty" $ insideSurface p1 [] `shouldBe` False + it "return False if polygon does not define at least a triangle" $ + insideSurface p1 [p1, p2] `shouldBe` False + it "returns True if point is inside polygon" $ do + let polygon = [p1, p2, p4, p3] + insideSurface p5 polygon `shouldBe` True + it "returns False if point is inside polygon" $ do + let polygon = [p1, p2, p4, p3] + let p = antipode p5 + insideSurface p polygon `shouldBe` False + it "returns False if point is a vertex of the polygon" $ do + let polygon = [p1, p2, p4, p3] + insideSurface p1 polygon `shouldBe` False + it "handles closed polygons" $ do + let polygon = [p1, p2, p4, p3, p1] + insideSurface p5 polygon `shouldBe` True + it "handles concave polygons" $ do + let malmo = decimalLatLong 55.6050 13.0038 + let ystad = decimalLatLong 55.4295 13.82 + let lund = decimalLatLong 55.7047 13.1910 + let helsingborg = decimalLatLong 56.0465 12.6945 + let kristianstad = decimalLatLong 56.0294 14.1567 + let polygon = [malmo, ystad, kristianstad, helsingborg, lund] + let hoor = decimalLatLong 55.9295 13.5297 + let hassleholm = decimalLatLong 56.1589 13.7668 + insideSurface hoor polygon `shouldBe` True + insideSurface hassleholm polygon `shouldBe` False + describe "intersections" $ do + it "returns nothing if both great circle are equals" $ do + let gc = greatCircle (decimalLatLong 51.885 0.235, decimalDegrees 108.63) + (intersections gc gc :: Maybe (LatLong, LatLong)) `shouldBe` Nothing + it "returns nothing if both great circle are equals (opposite orientation)" $ do + let gc1 = greatCircle (decimalLatLong 51.885 0.235, decimalLatLong 52.885 1.235) + let gc2 = greatCircle (decimalLatLong 52.885 1.235, decimalLatLong 51.885 0.235) + (intersections gc1 gc2 :: Maybe (LatLong, LatLong)) `shouldBe` Nothing + it "returns the two points where the two great circles intersects" $ do + let gc1 = greatCircle (decimalLatLong 51.885 0.235, decimalDegrees 108.63) + let gc2 = greatCircle (decimalLatLong 49.008 2.549, decimalDegrees 32.72) + let (i1, i2) = fromJust (intersections gc1 gc2) + i1 `shouldBe` decimalLatLong 50.9017226 4.4942782 + i2 `shouldBe` antipode i1 + describe "mean" $ do + it "returns Nothing if no point is given" $ (mean [] :: Maybe NVector) `shouldBe` Nothing + it "returns the unique given point" $ do + let p = readLatLong "500359N0054253W" + mean [p] `shouldBe` Just p + it "returns the geographical mean" $ do + let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 15000.0) + let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 25000.0) + let e = decimalLatLongHeight 54.3622869 (-4.5306725) zero + mean [p1, p2] `shouldBe` Just e + it "returns Nothing if list contains antipodal points" $ do + let points = + [ decimalLatLong 45 1 + , decimalLatLong 45 2 + , decimalLatLong 46 2 + , decimalLatLong 46 1 + , antipode (decimalLatLong 45 2) + ] + mean points `shouldBe` Nothing + describe "surfaceDistance" $ do + it "returns 0 if both points are equal" $ do + let p = readLatLong "500359N1795959W" + surfaceDistance p p r84 `shouldBe` zero + it "returns the distance between 2 points" $ do + let p1 = readLatLong "500359N0054253W" + let p2 = readLatLong "583838N0030412W" + surfaceDistance84 p1 p2 `shouldBe` metres 968854.8685 + it "handles singularity at the pole" $ + surfaceDistance northPole southPole r84 `shouldBe` kilometres 20015.114352200002 + it "handles the discontinuity at the Date Line" $ do + let p1 = readLatLong "500359N1795959W" + let p2 = readLatLong "500359N1795959E" + surfaceDistance p1 p2 (meanRadius wgs84) `shouldBe` metres 39.6596
+ test/Data/Geo/Jord/KinematicsSpec.hs view
@@ -0,0 +1,191 @@+module Data.Geo.Jord.KinematicsSpec+ ( spec+ ) where++import Data.Geo.Jord+import Data.Maybe (fromJust)+import Test.Hspec++spec :: Spec+spec =+ describe "kinematics" $ do+ describe "position" $ do+ it "computes position at t from p0, bearing and speed" $ do+ let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000)+ let p1 = decimalLatLongHeight 53.1882691 0.1332741 (metres 15000)+ let t = Track p0 (decimalDegrees 96.0217) (kilometresPerHour 124.8)+ position84 t (hours 1) `shouldBe` p1+ it "handles poles" $+ -- distance between poles assuming a spherical earth (WGS84) = 20015.114352200002km + -- track at north pole travelling at 20015.114352200002km/h and true north reaches the + -- south pole after 1 hour. + do+ let t = Track (decimalLatLong 90 0) zero (kilometresPerHour 20015.114352200002)+ position84 t (hours 1) `shouldBe` decimalLatLong (-90) 180.0+ it "return p0 if speed is 0" $ do+ let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000)+ let t = Track p0 (decimalDegrees 96.0217) zero+ position84 t (hours 1) `shouldBe` p0+ it "return p0 if duration is 0" $ do+ let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000)+ let t = Track p0 (decimalDegrees 96.0217) (kilometresPerHour 124.8)+ position84 t zero `shouldBe` p0+ describe "cpa" $ do+ it "handles trailing tracks" $ do+ let p1 = decimalLatLong 20 30+ let px = destination84 p1 (decimalDegrees 20) (kilometres 1)+ let p2 = interpolate p1 px 0.25+ let b1 = fromJust (initialBearing p1 px)+ let b2 = fromJust (initialBearing p2 px)+ let t1 = Track p1 b1 (knots 400)+ let t2 = Track p2 b2 (knots 400)+ let c = cpa84 t1 t2+ -- any time is correct but it should be close to zero since that's + -- our initial value + fmap (\r -> toMilliseconds (cpaTime r) < 5000) c `shouldBe` Just True+ fmap cpaDistance c `shouldBe` Just (metres 250.0036)+ it "handles heading tracks" $ do+ let p1 = decimalLatLong 20 30+ let p2 = decimalLatLong 21 31+ let b1 = fromJust (initialBearing p1 p2)+ let b2 = fromJust (initialBearing p2 p1)+ let t1 = Track p1 b1 (knots 400)+ let t2 = Track p2 b2 (knots 400)+ let c = cpa84 t1 t2+ -- distance between p1 and p2 = 152.354309 km + -- speed = 740.8 km/h + -- time = 152.354309 / 740.8 / 2 + fmap cpaTime c `shouldBe` Just (milliseconds 370191)+ fmap cpaDistance c `shouldBe` Just zero+ it "handles tracks at the same position" $ do+ let p = decimalLatLong 20 30+ let t1 = Track p (decimalDegrees 45) (knots 300)+ let t2 = Track p (decimalDegrees 135) (knots 500)+ let c = cpa84 t1 t2+ fmap cpaTime c `shouldBe` Just zero+ fmap cpaDistance c `shouldBe` Just zero+ it "computes time to CPA, positions and distance at CPA" $ do+ let p1 = decimalLatLong 20 (-60)+ let b1 = decimalDegrees 10+ let s1 = knots 15+ let p2 = decimalLatLong 34 (-50)+ let b2 = decimalDegrees 220+ let s2 = knots 300+ let t1 = Track p1 b1 s1+ let t2 = Track p2 b2 s2+ let c = cpa84 t1 t2+ fmap cpaTime c `shouldBe` Just (milliseconds 11396155)+ fmap cpaDistance c `shouldBe` Just (kilometres 124.2317453)+ it "returns Nothing if time to CPA is in the past" $ do+ let t1 = Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400)+ let t2 = Track (decimalLatLong 30.01 30) (decimalDegrees 315) (knots 400)+ cpa84 t1 t2 `shouldBe` Nothing+ describe "intercept" $ do+ it "returns Nothing if target and interceptor are at the same position" $+ intercept84+ (Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400))+ (decimalLatLong 30 30) `shouldBe`+ Nothing+ it "returns Nothing if interceptor is on the great circle of target and behind" $ do+ -- minimum speed would be ideally target speed + epsillon.+ let ip = decimalLatLong 20 30+ let px = destination84 ip (decimalDegrees 20) (kilometres 1)+ let tp = interpolate ip px 0.25+ let b = fromJust (initialBearing tp px)+ let t = Track tp b (knots 400)+ intercept84 t ip `shouldBe` Nothing+ it "handles interceptor on the great circle of target and in front" $ do+ let tp = decimalLatLong 20 30+ let px = destination84 tp (decimalDegrees 20) (kilometres 1)+ let ip = interpolate tp px 0.25+ let b = fromJust (initialBearing tp px)+ let t = Track tp b (knots 400)+ let i = intercept84 t ip+ fmap interceptorSpeed i `shouldBe` Just zero+ fmap interceptPosition i `shouldBe` Just ip+ fmap interceptTime i `shouldBe` Just (seconds 1.215)+ it "returns the minimum speed required for intercept to take place" $ do+ let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600)+ let ip = decimalLatLong 20 (-60)+ let i = intercept84 t ip+ fmap interceptorSpeed i `shouldBe` Just (knots 52.837096)+ fmap interceptTime i `shouldBe` Just (seconds 5947.698)+ let interceptor =+ Track+ ip+ (fromJust (fmap interceptorBearing i))+ (fromJust (fmap interceptorSpeed i))+ fmap interceptPosition i `shouldBe`+ Just (position84 interceptor (fromJust (fmap interceptTime i)))+ describe "interceptBySpeed" $ do+ it "returns Nothing if target and interceptor are at the same position" $+ interceptBySpeed84+ (Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400))+ (decimalLatLong 30 30)+ (knots 400) `shouldBe`+ Nothing+ it "returns Nothing if interceptor speed is below minimum speed" $ do+ let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600)+ let ip = decimalLatLong 20 (-60)+ interceptBySpeed84 t ip (knots 50) `shouldBe` Nothing+ it "handles interceptor on the great circle of target and behind" $ do+ let ip = decimalLatLong 20 30+ let px = destination84 ip (decimalDegrees 20) (kilometres 1)+ let tp = interpolate ip px 0.25+ let b = fromJust (initialBearing tp px)+ let t = Track tp b (metresPerSecond 400)+ let i = interceptBySpeed84 t ip (metresPerSecond 500)+ fmap interceptTime i `shouldBe` Just (seconds 2.5)+ it "returns the time needed for intercept to take place" $ do+ let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600)+ let ip = decimalLatLong 20 (-60)+ let i = interceptBySpeed84 t ip (knots 700)+ fmap interceptTime i `shouldBe` Just (seconds 2764.692)+ fmap interceptorBearing i `shouldBe` Just (decimalDegrees 25.93541277)+ fmap interceptDistance i `shouldBe` Just (kilometres 995.5960805999999)+ describe "interceptByTime" $ do+ it "returns Nothing if duration is zero" $+ interceptByTime84+ (Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400))+ (decimalLatLong 34 (-50))+ zero `shouldBe`+ Nothing+ it "returns Nothing if duration is negative" $+ interceptByTime84+ (Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400))+ (decimalLatLong 34 (-50))+ (seconds (-1)) `shouldBe`+ Nothing+ it "returns Nothing if target and interceptor are at the same position" $+ interceptByTime84+ (Track (decimalLatLong 30 30) (decimalDegrees 45) (knots 400))+ (decimalLatLong 30 30)+ (seconds 10) `shouldBe`+ Nothing+ it "returns the speed needed for intercept to take place" $ do+ let t = Track (decimalLatLong 34 (-50)) (decimalDegrees 220) (knots 600)+ let ip = decimalLatLong 20 (-60)+ let d = seconds 2700+ let i = interceptByTime84 t ip d+ fmap interceptorSpeed i `shouldBe` Just (knots 730.959238)+ fmap interceptorBearing i `shouldBe` Just (decimalDegrees 26.1199030)+ fmap interceptPosition i `shouldBe` Just (decimalLatLong 28.1366797 (-55.4559475))+ fmap interceptDistance i `shouldBe` Just (metres 1015302.3815)+ fmap interceptTime i `shouldBe` Just (seconds 2700)+ it "handles the poles" $ do+ -- distance between poles assuming a spherical earth (WGS84) = 20015.114352200002km+ -- target at north pole travelling at 500km/h and true north can be intercepted from+ -- the south pole by an interceptor travelling at ~ 19515.114352200002km/h and 180 degrees.+ let t = Track (decimalLatLong 90 0) zero (kilometresPerHour 500)+ let ip = decimalLatLong (-90) 0+ let i = interceptByTime84 t ip (seconds 3600)+ fmap interceptorSpeed i `shouldBe` Just (kilometresPerHour 19515.11434)+ fmap interceptorBearing i `shouldBe` Just (decimalDegrees 180)+ it "handles the interceptor being at the intercept position at t" $ do+ let tp = decimalLatLong 34 (-50)+ let t = Track tp (decimalDegrees 220) (knots 600)+ let d = seconds 3600+ let ip = position84 t d+ let i = interceptByTime84 t ip d+ fmap interceptorSpeed i `shouldBe` Just zero+ fmap interceptorBearing i `shouldBe` initialBearing ip tp
test/Data/Geo/Jord/LatLongSpec.hs view
@@ -1,63 +1,63 @@-module Data.Geo.Jord.LatLongSpec- ( spec- ) where--import Data.Geo.Jord-import Test.Hspec--spec :: Spec-spec = do- describe "Reading valid DMS text" $ do- it "reads 553621N0130002E" $- readLatLong "553621N0130002E" `shouldBe` decimalLatLong 55.6058333 13.0005555- it "reads 55°36'21''N 013°00'02''E" $- readLatLong "55°36'21''N 013°00'02''E" `shouldBe` decimalLatLong 55.6058333 13.0005555- it "reads 5536N01300E" $ readLatLong "5536N01300E" `shouldBe` decimalLatLong 55.6 13.0- it "reads 55N013E" $ readLatLong "55N013E" `shouldBe` decimalLatLong 55.0 13.0- it "reads 011659S0364900E" $- readLatLong "011659S0364900E" `shouldBe` decimalLatLong (-1.2830555) 36.8166666- it "reads 0116S03649E" $- readLatLong "0116S03649E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666- it "reads 1°16'S,36°49'E" $- readLatLong "1°16'S,36°49'E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666- it "reads 01S036E" $ readLatLong "01S036E" `shouldBe` decimalLatLong (-1.0) 36.0- it "reads 473622N1221955W" $- readLatLong "473622N1221955W" `shouldBe` decimalLatLong 47.6061111 (-122.3319444)- it "reads 4736N12219W" $- readLatLong "4736N12219W" `shouldBe` decimalLatLong 47.6 (-122.3166666)- it "reads 47N122W" $ readLatLong "47N122W" `shouldBe` decimalLatLong 47.0 (-122.0)- it "reads 47°N 122°W" $ readLatLong "47°N 122°W" `shouldBe` decimalLatLong 47.0 (-122.0)- it "reads 544807S0681811W" $- readLatLong "544807S0681811W" `shouldBe` decimalLatLong (-54.8019444) (-68.3030555)- it "reads 5448S06818W" $ readLatLong "5448S06818W" `shouldBe` decimalLatLong (-54.8) (-68.3)- it "reads 54S068W" $ readLatLong "54S068W" `shouldBe` decimalLatLong (-54.0) (-68.0)- describe "Reading invalid DMS text" $ do- it "fails to read 553621K0130002E" $- readLatLongE "553621K0130002E" `shouldBe` Left "couldn't read geo pos 553621K0130002E"- it "fails to read 011659S0364900Z" $- readLatLongE "011659S0364900Z" `shouldBe` Left "couldn't read geo pos 011659S0364900Z"- it "fails to read 4736221221955W" $- readLatLongE "4736221221955W" `shouldBe` Left "couldn't read geo pos 4736221221955W"- it "fails to read 54480S0681811W" $- readLatLongE "54480S0681811W" `shouldBe` Left "couldn't read geo pos 54480S0681811W"- it "fails to read 553621N013000E" $- readLatLongE "553621N013000E" `shouldBe` Left "couldn't read geo pos 553621N013000E"- it "fails to read 914807S0681811W" $- readLatLongE "914807S0681811W" `shouldBe` Left "couldn't read geo pos 914807S0681811W"- it "fails to read 544807S1811811W" $- readLatLongE "544807S1811811W" `shouldBe` Left "couldn't read geo pos 544807S1811811W"- it "fails to read 546007S1801811W" $- readLatLongE "546007S1801811W" `shouldBe` Left "couldn't read geo pos 546007S1801811W"- it "fails to read 545907S1801860W" $- readLatLongE "545907S1801860W" `shouldBe` Left "couldn't read geo pos 545907S1801860W"- describe "Showing geographic positions" $ do- it "shows the N/E position formatted in DMS with symbols" $- show (decimalLatLong 55.60583333 13.00055556) `shouldBe` "55°36'21.000\"N,13°0'2.000\"E"- it "shows the S/E position formatted in DMS with symbols" $- show (decimalLatLong (-1.28305556) 36.81666) `shouldBe` "1°16'59.000\"S,36°48'59.976\"E"- it "shows the N/W position formatted in DMS with symbols" $- show (decimalLatLong 47.60611 (-122.33194)) `shouldBe`- "47°36'21.996\"N,122°19'54.984\"W"- it "shows the S/W position formatted in DMS with symbols" $- show (decimalLatLong (-54.80194) (-68.30305)) `shouldBe`- "54°48'6.984\"S,68°18'10.980\"W"+module Data.Geo.Jord.LatLongSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid DMS text" $ do + it "reads 553621N0130002E" $ + readLatLong "553621N0130002E" `shouldBe` decimalLatLong 55.6058333 13.0005555 + it "reads 55°36'21''N 013°00'02''E" $ + readLatLong "55°36'21''N 013°00'02''E" `shouldBe` decimalLatLong 55.6058333 13.0005555 + it "reads 5536N01300E" $ readLatLong "5536N01300E" `shouldBe` decimalLatLong 55.6 13.0 + it "reads 55N013E" $ readLatLong "55N013E" `shouldBe` decimalLatLong 55.0 13.0 + it "reads 011659S0364900E" $ + readLatLong "011659S0364900E" `shouldBe` decimalLatLong (-1.2830555) 36.8166666 + it "reads 0116S03649E" $ + readLatLong "0116S03649E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666 + it "reads 1°16'S,36°49'E" $ + readLatLong "1°16'S,36°49'E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666 + it "reads 01S036E" $ readLatLong "01S036E" `shouldBe` decimalLatLong (-1.0) 36.0 + it "reads 473622N1221955W" $ + readLatLong "473622N1221955W" `shouldBe` decimalLatLong 47.6061111 (-122.3319444) + it "reads 4736N12219W" $ + readLatLong "4736N12219W" `shouldBe` decimalLatLong 47.6 (-122.3166666) + it "reads 47N122W" $ readLatLong "47N122W" `shouldBe` decimalLatLong 47.0 (-122.0) + it "reads 47°N 122°W" $ readLatLong "47°N 122°W" `shouldBe` decimalLatLong 47.0 (-122.0) + it "reads 544807S0681811W" $ + readLatLong "544807S0681811W" `shouldBe` decimalLatLong (-54.8019444) (-68.3030555) + it "reads 5448S06818W" $ readLatLong "5448S06818W" `shouldBe` decimalLatLong (-54.8) (-68.3) + it "reads 54S068W" $ readLatLong "54S068W" `shouldBe` decimalLatLong (-54.0) (-68.0) + describe "Reading invalid DMS text" $ do + it "fails to read 553621K0130002E" $ + readLatLongE "553621K0130002E" `shouldBe` Left "couldn't read geo pos 553621K0130002E" + it "fails to read 011659S0364900Z" $ + readLatLongE "011659S0364900Z" `shouldBe` Left "couldn't read geo pos 011659S0364900Z" + it "fails to read 4736221221955W" $ + readLatLongE "4736221221955W" `shouldBe` Left "couldn't read geo pos 4736221221955W" + it "fails to read 54480S0681811W" $ + readLatLongE "54480S0681811W" `shouldBe` Left "couldn't read geo pos 54480S0681811W" + it "fails to read 553621N013000E" $ + readLatLongE "553621N013000E" `shouldBe` Left "couldn't read geo pos 553621N013000E" + it "fails to read 914807S0681811W" $ + readLatLongE "914807S0681811W" `shouldBe` Left "couldn't read geo pos 914807S0681811W" + it "fails to read 544807S1811811W" $ + readLatLongE "544807S1811811W" `shouldBe` Left "couldn't read geo pos 544807S1811811W" + it "fails to read 546007S1801811W" $ + readLatLongE "546007S1801811W" `shouldBe` Left "couldn't read geo pos 546007S1801811W" + it "fails to read 545907S1801860W" $ + readLatLongE "545907S1801860W" `shouldBe` Left "couldn't read geo pos 545907S1801860W" + describe "Showing geographic positions" $ do + it "shows the N/E position formatted in DMS with symbols" $ + show (decimalLatLong 55.60583333 13.00055556) `shouldBe` "55°36'21.000\"N,13°0'2.000\"E" + it "shows the S/E position formatted in DMS with symbols" $ + show (decimalLatLong (-1.28305556) 36.81666) `shouldBe` "1°16'59.000\"S,36°48'59.976\"E" + it "shows the N/W position formatted in DMS with symbols" $ + show (decimalLatLong 47.60611 (-122.33194)) `shouldBe` + "47°36'21.996\"N,122°19'54.984\"W" + it "shows the S/W position formatted in DMS with symbols" $ + show (decimalLatLong (-54.80194) (-68.30305)) `shouldBe` + "54°48'6.984\"S,68°18'10.980\"W"
test/Data/Geo/Jord/LengthSpec.hs view
@@ -1,35 +1,40 @@-module Data.Geo.Jord.LengthSpec- ( spec- ) where--import Data.Geo.Jord-import Test.Hspec--spec :: Spec-spec = do- describe "Reading valid lengths" $ do- it "reads -15.2m" $ readLength "-15.2m" `shouldBe` metres (-15.2)- it "reads 154km" $ readLength "154km" `shouldBe` kilometres 154- it "reads 1000Nm" $ readLength "1000Nm" `shouldBe` nauticalMiles 1000- it "reads 25000ft" $ readLength "25000ft" `shouldBe` feet 25000- describe "Reading invalid lengths" $ do- it "fails to read 5" $ readLengthE "5" `shouldBe` Left "couldn't read length 5"- it "fails to read 5nmi" $ readLengthE "5nmi" `shouldBe` Left "couldn't read length 5nmi"- describe "Showing lengths" $ do- it "shows length in metres when <= 10000 m" $ show (metres 5) `shouldBe` "5.0m"- it "shows length in kilometres when > 10000 m" $- show (kilometres 1000) `shouldBe` "1000.0km"- describe "Converting lengths" $ do- it "converts metres to kilometres" $ toKilometres (metres 1000) `shouldBe` 1.0- it "converts metres to nautical miles" $- toNauticalMiles (metres 1000) `shouldBe` 0.5399568034557235- it "converts kilometres to nautical miles" $- toNauticalMiles (kilometres 1000) `shouldBe` 539.9568034557235- it "converts nautical miles to metres" $ toMetres (nauticalMiles 10.5) `shouldBe` 19446- it "converts nautical miles to kilometres" $- toKilometres (nauticalMiles 10.5) `shouldBe` 19.446- it "converts feet to metres" $ toMetres (feet 25000) `shouldBe` 7620- it "converts metres to feet" $ toFeet (metres 7620) `shouldBe` 25000- describe "Adding/Subtracting lengths" $ do- it "adds lengths" $ add (kilometres 1000) (metres 1000) `shouldBe` metres 1001000- it "subtracts lengths" $ sub (metres 1000) (nauticalMiles 10.5) `shouldBe` metres (-18446)+module Data.Geo.Jord.LengthSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid lengths" $ do + it "reads -15.2m" $ readLength "-15.2m" `shouldBe` metres (-15.2) + it "reads 154km" $ readLength "154km" `shouldBe` kilometres 154 + it "reads 1000Nm" $ readLength "1000Nm" `shouldBe` nauticalMiles 1000 + it "reads 25000ft" $ readLength "25000ft" `shouldBe` feet 25000 + describe "Reading invalid lengths" $ do + it "fails to read 5" $ readLengthE "5" `shouldBe` Left "couldn't read length 5" + it "fails to read 5nmi" $ readLengthE "5nmi" `shouldBe` Left "couldn't read length 5nmi" + describe "Showing lengths" $ do + it "shows length in metres when <= 10000 m" $ show (metres 5) `shouldBe` "5.0m" + it "shows length in kilometres when > 10000 m" $ + show (kilometres 1000) `shouldBe` "1000.0km" + describe "Converting lengths" $ do + it "converts metres to kilometres" $ toKilometres (metres 1000) `shouldBe` 1.0 + it "converts metres to nautical miles" $ + toNauticalMiles (metres 1000) `shouldBe` 0.5399568034557235 + it "converts kilometres to nautical miles" $ + toNauticalMiles (kilometres 1000) `shouldBe` 539.9568034557235 + it "converts nautical miles to metres" $ toMetres (nauticalMiles 10.5) `shouldBe` 19446 + it "converts nautical miles to kilometres" $ + toKilometres (nauticalMiles 10.5) `shouldBe` 19.446 + it "converts feet to metres" $ toMetres (feet 25000) `shouldBe` 7620 + it "converts metres to feet" $ toFeet (metres 7620) `shouldBe` 25000 + describe "Resolution" $ do + it "handles 1 kilometre" $ toKilometres (kilometres 1) `shouldBe` 1 + it "handles 1 metre" $ toMetres (metres 1) `shouldBe` 1 + it "handles 1 nautical mile" $ toNauticalMiles (nauticalMiles 1) `shouldBe` 1 + it "handles 1 foot" $ toFeet (feet 1) `shouldBe` 1 + describe "Adding/Subtracting lengths" $ do + it "adds lengths" $ add (kilometres 1000) (metres 1000) `shouldBe` metres 1001000 + it "subtracts lengths" $ sub (metres 1000) (nauticalMiles 10.5) `shouldBe` metres (-18446)
+ test/Data/Geo/Jord/SpeedSpec.hs view
@@ -0,0 +1,41 @@+module Data.Geo.Jord.SpeedSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid speeds" $ do + it "reads -15.2m/s" $ readSpeed "-15.2m/s" `shouldBe` metresPerSecond (-15.2) + it "reads 154km/h" $ readSpeed "154km/h" `shouldBe` kilometresPerHour 154 + it "reads 200mph" $ readSpeed "200mph" `shouldBe` milesPerHour 200 + it "reads 400kt" $ readSpeed "400kt" `shouldBe` knots 400 + it "reads 1ft/s" $ readSpeed "1ft/s" `shouldBe` feetPerSecond 1 + describe "Reading invalid speeds" $ do + it "fails to read 5" $ readSpeedE "5" `shouldBe` Left "couldn't read speed 5" + it "fails to read 5mps" $ readSpeedE "5mps" `shouldBe` Left "couldn't read speed 5mps" + describe "Showing speeds" $ + it "shows speed in kilometres per hour" $ + show (kilometresPerHour 154) `shouldBe` "154.0km/h" + describe "Converting speeds" $ do + it "converts metres per seconds to kilometres per hour" $ + toKilometresPerHour (metresPerSecond 100) `shouldBe` 360.0 + it "converts metres per seconds to miles per hour" $ + toMilesPerHour (metresPerSecond 100) `shouldBe` 223.69362920544023 + it "converts kilometres per hour to knots" $ + toKnots (kilometresPerHour 1000) `shouldBe` 539.9568034557235 + it "converts feet per second to kilometres per hour" $ + toKilometresPerHour (feetPerSecond 1) `shouldBe` 1.09728 + describe "Resolution" $ do + it "handles 1 km/h" $ toKilometresPerHour (kilometresPerHour 1) `shouldBe` 1 + it "handles 1 1m/s" $ toMetresPerSecond (metresPerSecond 1) `shouldBe` 1 + it "handles 1 1mph" $ toMilesPerHour (milesPerHour 1) `shouldBe` 1 + it "handles 1 knot" $ toKnots (knots 1) `shouldBe` 1 + it "handles 1 fp/s" $ toFeetPerSecond (feetPerSecond 1) `shouldBe` 1 + describe "Adding/Subtracting speeds" $ do + it "adds speeds" $ + add (kilometresPerHour 1000) (metresPerSecond 1000) `shouldBe` kilometresPerHour 4600 + it "subtracts lengths" $ + sub (metresPerSecond 1000) (knots 10.5) `shouldBe` kilometresPerHour 3580.554
− test/Data/Geo/Jord/TransformSpec.hs
@@ -1,69 +0,0 @@-module Data.Geo.Jord.TransformSpec- ( spec- ) where--import Data.Geo.Jord-import Test.Hspec--spec :: Spec-spec = do- describe "Ellipsoidal transformation between coordinates systems" $ do- it "transforms NVector position to ECEF position" $ do- let p = nvector 0.5 0.5 0.7071- toEcef p wgs84 `shouldBe` ecefMetres 3194434.411 3194434.411 4487326.819- it "transforms angular position to ECEF position" $ do- let refAngular =- [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834)- , decimalLatLongHeight 45.0 45.0 zero- , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)- ]- let refEcefs =- [ ecefMetres 5733855.775 (-6370998.38) 7008137.511- , ecefMetres 3194419.145 3194419.145 4487348.409- , ecefMetres 4200996.77 172460.321 4780102.808- ]- mapM_ (\(a, e) -> toEcef a wgs84 `shouldBe` e) (zip refAngular refEcefs)- it "transforms ECEF position to angular position" $ do- let refAngular =- [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834050887)- , decimalLatLongHeight 45.0 45.0 zero- , decimalLatLongHeight 48.8562 2.3508 (metres 67.36990469945641)- ]- let refEcefs =- [ ecefMetres 5733855.774881717 (-6370998.380260889) 7008137.510624695- , ecefMetres 3194419.145121972 3194419.145121971 4487348.408606014- , ecefMetres 4200996.769831858 172460.320727757 4780102.807914356- ]- mapM_ (\(a, e) -> fromEcef e wgs84 `shouldBe` a) (zip refAngular refEcefs)- describe "Spherical transformation between coordinates systems" $ do- it "transforms NVector position to ECEF position" $ do- let p = nvector 0.5 0.5 0.7071- toEcef p s84 `shouldBe` ecefMetres 3185519.660103391 3185519.660103391 4504961.903318216- it "transforms angular position to ECEF position" $ do- let refAngular =- [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834)- , decimalLatLongHeight 45.0 45.0 zero- , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)- , latLongHeight (readLatLong "531914N0014347W") (metres 15000.0)- , decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0)- ]- let refEcefs =- [ ecefMetres 5725717.354041086 (-6361955.622990872) 7025277.913631903- , ecefMetres 3185504.385500001 3185504.3855 4504983.504973072- , ecefMetres 4188328.8912726147 171940.27595767862 4797806.669141033- , ecefMetres 3812864.094233316 (-115142.863124558) 5121515.160893968- , ecefMetres 3826406.4642097903 8900.535428827865 5112694.238306408- ]- mapM_ (\(a, e) -> toEcef a s84 `shouldBe` e) (zip refAngular refEcefs)- it "transforms ECEF position to angular position" $ do- let refAngular =- [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834217537)- , decimalLatLongHeight 45.0 45.0 (metres 1e-3)- , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)- ]- let refEcefs =- [ ecefMetres 5725717.354 (-6361955.623) 7025277.914- , ecefMetres 3185504.386 3185504.386 4504983.505- , ecefMetres 4188328.891 171940.276 4797806.669- ]- mapM_ (\(a, e) -> fromEcef e s84 `shouldBe` a) (zip refAngular refEcefs)
+ test/Data/Geo/Jord/TransformationSpec.hs view
@@ -0,0 +1,69 @@+module Data.Geo.Jord.TransformationSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Ellipsoidal transformation between coordinates systems" $ do + it "transforms NVector position to ECEF position" $ do + let p = nvector 0.5 0.5 0.7071 + toEcef p wgs84 `shouldBe` ecefMetres 3194434.411 3194434.411 4487326.8195 + it "transforms angular position to ECEF position" $ do + let refAngular = + [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834) + , decimalLatLongHeight 45.0 45.0 zero + , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099) + ] + let refEcefs = + [ ecefMetres 5733855.7748 (-6370998.3802) 7008137.5108 + , ecefMetres 3194419.1451 3194419.1451 4487348.4088 + , ecefMetres 4200996.7697 172460.3207 4780102.808 + ] + mapM_ (\(a, e) -> toEcef a wgs84 `shouldBe` e) (zip refAngular refEcefs) + it "transforms ECEF position to angular position" $ do + let refAngular = + [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.8339) + , decimalLatLongHeight 45.0 45.0 (metres (-0.0001)) + , decimalLatLongHeight 48.8562 2.3508 (metres 67.3697) + ] + let refEcefs = + [ ecefMetres 5733855.774881717 (-6370998.380260889) 7008137.510624695 + , ecefMetres 3194419.145121972 3194419.145121971 4487348.408606014 + , ecefMetres 4200996.769831858 172460.320727757 4780102.807914356 + ] + mapM_ (\(a, e) -> fromEcef e wgs84 `shouldBe` a) (zip refAngular refEcefs) + describe "Spherical transformation between coordinates systems" $ do + it "transforms NVector position to ECEF position" $ do + let p = nvector 0.5 0.5 0.7071 + toEcef p s84 `shouldBe` ecefMetres 3185519.6603 3185519.6603 4504961.9036 + it "transforms angular position to ECEF position" $ do + let refAngular = + [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834) + , decimalLatLongHeight 45.0 45.0 zero + , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099) + , latLongHeight (readLatLong "531914N0014347W") (metres 15000.0) + , decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0) + ] + let refEcefs = + [ ecefMetres 5725717.3542 (-6361955.6232) 7025277.9139 + , ecefMetres 3185504.3857 3185504.3857 4504983.5053 + , ecefMetres 4188328.8913 171940.276 4797806.6692 + , ecefMetres 3812864.0945 (-115142.8631) 5121515.1612 + , ecefMetres 3826406.4644 8900.5354 5112694.2386 + ] + mapM_ (\(a, e) -> toEcef a s84 `shouldBe` e) (zip refAngular refEcefs) + it "transforms ECEF position to angular position" $ do + let refAngular = + [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.8338) + , decimalLatLongHeight 45.0 45.0 (metres 1e-4) + , decimalLatLongHeight 48.8562 2.3508 (metres 67.3693) + ] + let refEcefs = + [ ecefMetres 5725717.354 (-6361955.623) 7025277.914 + , ecefMetres 3185504.386 3185504.386 4504983.505 + , ecefMetres 4188328.891 171940.276 4797806.669 + ] + mapM_ (\(a, e) -> fromEcef e s84 `shouldBe` a) (zip refAngular refEcefs)