jord 0.4.2.0 → 0.5.0.0
raw patch · 25 files changed
+1581/−1194 lines, 25 filesdep +criterionnew-component:exe:jord-benchmarks
Dependencies added: criterion
Files
- ChangeLog.md +46/−36
- app/Eval.hs +45/−8
- app/Main.hs +206/−201
- app/Show.hs +96/−90
- app/State.hs +1/−0
- benchmarks/GeodeticsBG.hs +37/−0
- benchmarks/KinematicsBG.hs +45/−0
- benchmarks/Main.hs +9/−0
- benchmarks/TransformationBG.hs +36/−0
- jord.cabal +19/−3
- src/Data/Geo/Jord.hs +1/−1
- src/Data/Geo/Jord/Angle.hs +14/−11
- src/Data/Geo/Jord/Duration.hs +2/−2
- src/Data/Geo/Jord/Geodetics.hs +196/−92
- src/Data/Geo/Jord/Internal.hs +37/−0
- src/Data/Geo/Jord/Kinematics.hs +143/−140
- src/Data/Geo/Jord/LatLong.hs +1/−1
- src/Data/Geo/Jord/Length.hs +1/−1
- src/Data/Geo/Jord/Parse.hs +0/−52
- src/Data/Geo/Jord/Parser.hs +52/−0
- src/Data/Geo/Jord/Speed.hs +1/−1
- src/Data/Geo/Jord/Transformation.hs +141/−141
- test/Data/Geo/Jord/GeodeticsSpec.hs +251/−213
- test/Data/Geo/Jord/KinematicsSpec.hs +199/−199
- test/Data/Geo/Jord/SpeedSpec.hs +2/−2
ChangeLog.md view
@@ -1,36 +1,46 @@-### 0.4.2.0 - -- Fixed intercept -- jord-exe renamed jord-repl - -### 0.4.1.0 - -- Fixed interceptBySpeed -- Nautical miles symbol is "nm" -- REPL: intercept for intercept, interceptBySpeed and interceptByTime -- REPL: show length and speed in user selected unit - -### 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 -- Added Frames (Body, Local, North East Down) -- Added delta and target from position(s), frame and earth model -- Added earth models (WGS84, WGS72, GRS80 and derived spherical models) -- Builds against LTS 12.2 (GHC 8.4.3) and LTS 11.18 (GHC 8.2.2) - -### 0.2.0.0 - -- GeoPos -> LatLong -- Split Position from GreatCircle -- require base >= 4.9 - -### 0.1.0.0 - -- Initial version +### 0.5.0.0++- Added Benchmarks+- Added GreatArc+- Added GreatArc from tuple of positions+- Added GreatArc from GreatCircle+- Added GreatArc from Track and Duration+- Added alongTrackDistance+- Added GreatArcs intersection++### 0.4.2.0++- Fixed intercept+- jord-exe renamed jord-repl++### 0.4.1.0++- Fixed interceptBySpeed+- Nautical miles symbol is "nm"+- REPL: intercept for intercept, interceptBySpeed and interceptByTime+- REPL: show length and speed in user selected unit++### 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+- Added Frames (Body, Local, North East Down)+- Added delta and target from position(s), frame and earth model+- Added earth models (WGS84, WGS72, GRS80 and derived spherical models)+- Builds against LTS 12.2 (GHC 8.4.3) and LTS 11.18 (GHC 8.2.2)++### 0.2.0.0++- GeoPos -> LatLong+- Split Position from GreatCircle+- require base >= 4.9++### 0.1.0.0++- Initial version
app/Eval.hs view
@@ -96,7 +96,8 @@ -- | All supported functions. functions :: [String] functions = - [ "antipode" + [ "alongTrackDistance" + , "antipode" , "crossTrackDistance" , "cpa" , "delta" @@ -109,12 +110,14 @@ , "finalBearing" , "fromEcef" , "geo" + , "greatArc" , "greatCircle" , "initialBearing" , "intercept" , "interpolate" + , "intersection" , "intersections" - , "insideSurface" + , "isInsideSurface" , "mean" , "ned" , "nedBetween" @@ -143,6 +146,10 @@ Just (Gp geo) -> Right (Np (toNVector geo)) Just v -> Right v Nothing -> tryRead p +evalExpr (AlongTrackDistance a b) state = + case [evalExpr a state, evalExpr b state] of + [Right (Np p), Right (Ga ga)] -> Right (Len (alongTrackDistance84 p ga)) + r -> Left ("Call error: alongTrackDistance84 " ++ showErr r state) evalExpr (Antipode a) state = case evalExpr a state of (Right (Np p)) -> Right (Np (antipode p)) @@ -221,10 +228,16 @@ r -> Left ("Call error: geo " ++ showErr r state) where vs = map (`evalExpr` state) as +evalExpr (GreatArcE as) state = + case fmap (`evalExpr` state) as of + [Right (Np p1), Right (Np p2)] -> bimap id Ga (greatArcE (p1, p2)) + [Right (Trk t), Right (Dur d)] -> bimap id Ga (greatArcE (t, d)) + r -> Left ("Call error: greatArc " ++ showErr r state) evalExpr (GreatCircleE as) state = case fmap (`evalExpr` state) 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 (Ga ga)] -> bimap id Gc (greatCircleE ga) [Right (Trk t)] -> bimap id Gc (greatCircleE t) r -> Left ("Call error: greatCircle " ++ showErr r state) evalExpr (InitialBearing a b) state = @@ -248,6 +261,14 @@ case [evalExpr a state, evalExpr b state] of [Right (Np p1), Right (Np p2)] -> Right (Np (interpolate p1 p2 c)) r -> Left ("Call error: interpolate " ++ showErr r state) +evalExpr (Intersection a b) state = + case [evalExpr a state, evalExpr b state] of + [Right (Ga ga1), Right (Ga ga2)] -> + maybe + (Left "no great arcs intersection") + (Right . Np) + (intersection ga1 ga2 :: Maybe (AngularPosition NVector)) + r -> Left ("Call error: intersection " ++ showErr r state) evalExpr (Intersections a b) state = case [evalExpr a state, evalExpr b state] of [Right (Gc gc1), Right (Gc gc2)] -> @@ -256,12 +277,12 @@ (\is -> Right (Vals [Np (fst is), Np (snd is)])) (intersections gc1 gc2 :: Maybe (AngularPosition NVector, AngularPosition NVector)) r -> Left ("Call error: intersections " ++ showErr r state) -evalExpr (InsideSurface as) state = +evalExpr (IsInsideSurface as) state = let m = map (`evalExpr` state) as 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 state) + then Right (Bool (isInsideSurface (head ps) (tail ps))) + else Left ("Call error: isInsideSurface " ++ showErr m state) evalExpr (Mean as) state = let m = map (`evalExpr` state) as ps = [p | Right (Np p) <- m] @@ -440,6 +461,8 @@ ------------------------------------- data Expr = Param String + | AlongTrackDistance Expr + Expr | Antipode Expr | ClosestPointOfApproach Expr Expr @@ -468,6 +491,7 @@ | FromEcef Expr String | Geo [Expr] + | GreatArcE [Expr] | GreatCircleE [Expr] | InitialBearing Expr Expr @@ -475,9 +499,11 @@ | Interpolate Expr Expr Double + | Intersection Expr + Expr | Intersections Expr Expr - | InsideSurface [Expr] + | IsInsideSurface [Expr] | Mean [Expr] | NedBetween Expr Expr @@ -505,6 +531,10 @@ deriving (Show) transform :: (MonadFail m) => Ast -> m Expr +transform (Call "alongTrackDistance" [e1, e2]) = do + p <- transform e1 + ga <- transform e2 + return (AlongTrackDistance p ga) transform (Call "antipode" [e]) = fmap Antipode (transform e) transform (Call "cpa" [e1, e2]) = do t1 <- transform e1 @@ -559,6 +589,9 @@ transform (Call "geo" e) = do ps <- mapM transform e return (Geo ps) +transform (Call "greatArc" e) = do + ps <- mapM transform e + return (GreatArcE ps) transform (Call "greatCircle" e) = do ps <- mapM transform e return (GreatCircleE ps) @@ -576,13 +609,17 @@ if d >= 0.0 && d <= 1.0 then return (Interpolate p1 p2 d) else fail "Semantic error: interpolate expects [0..1] as last argument" +transform (Call "intersection" [e1, e2]) = do + ga1 <- transform e1 + ga2 <- transform e2 + return (Intersection ga1 ga2) transform (Call "intersections" [e1, e2]) = do gc1 <- transform e1 gc2 <- transform e2 return (Intersections gc1 gc2) -transform (Call "insideSurface" e) = do +transform (Call "isInsideSurface" e) = do ps <- mapM transform e - return (InsideSurface ps) + return (IsInsideSurface ps) transform (Call "mean" e) = do ps <- mapM transform e return (Mean ps)
app/Main.hs view
@@ -8,206 +8,211 @@ -- -- REPL around "Jord". -- -module Main where - -import Data.Geo.Jord -import Data.List ((\\), dropWhileEnd, isPrefixOf) -import Eval -import Prelude hiding (lookup) -import Show -import State -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 Show+import State+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 emptyState - 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 -> State -> (Either String String, State) -evalS s state - | null s = (Right "", state) - | head s == ':' = evalC w state - | (v:"=":e) <- w = - if v `elem` functions - then (Left (v ++ " is a reserved keyword"), state) - else let r = eval (unwords e) state - state' = save r v state - in (showR r state', state') - | otherwise = (showR (eval s state) state, state) - where - w = words s - -evalC :: [String] -> State -> (Either String String, State) -evalC [":help"] state = (Right (help state), state) -evalC [":?"] state = (Right (help state), state) -evalC [":show", v] state = (evalShow v state, state) -evalC [":delete", v] state = (Right ("deleted var: " ++ v), delete v state) -evalC [":units", u1, u2] state = evalUnits [u1, u2] state -evalC [":units", u] state = evalUnits [u] state -evalC [":units"] state = showUnits state -evalC [":reset"] _ = (Right "REPL reset ", emptyState) -evalC c state = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), state) - -evalShow :: String -> State -> Either String String -evalShow n state = - maybe (Left ("Unbound variable: " ++ n)) (\v -> Right (showVar n v state)) (lookup n state) - -evalUnits :: [String] -> State -> (Either String String, State) -evalUnits us s = showUnits (setUnits us s) - -showUnits :: State -> (Either String String, State) -showUnits s = (Right ("Units:\n length = " ++ lengthUnit s ++ "\n speed = " ++ speedUnit s), s) - -help :: State -> String -help s = - "\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} show {var}\n" ++ - " :delete {var} delete {var}\n" ++ - " :units length speed set length and speed units used for display\n" ++ - " see supported length and speed format\n" ++ - " currently: length = " ++ - lengthUnit s ++ - "; speed = " ++ - speedUnit s ++ - "\n" ++ - " :units show length and speed units used for display\n" ++ - " :reset reset REPL to default state (including deleting all variables)\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: " ++ - showLength r84 s ++ - "\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: " ++ - showLength r84 s ++ - "\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" ++ - " intercept track pos spd time needed by interceptor at pos and travelling at spd to intercept target\n" ++ - " intercept 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 -> State -> State -save (Right v) k state = insert k v state -save _ _ state = state - -showR :: Result -> State -> Either String String -showR (Left err) _ = Left err -showR (Right v) s = Right (showV v s) - -showVar :: String -> Value -> State -> String -showVar n v s = n ++ "=" ++ showV v s + 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 emptyState+ 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 -> State -> (Either String String, State)+evalS s state+ | null s = (Right "", state)+ | head s == ':' = evalC w state+ | (v:"=":e) <- w =+ if v `elem` functions+ then (Left (v ++ " is a reserved keyword"), state)+ else let r = eval (unwords e) state+ state' = save r v state+ in (showR r state', state')+ | otherwise = (showR (eval s state) state, state)+ where+ w = words s++evalC :: [String] -> State -> (Either String String, State)+evalC [":help"] state = (Right (help state), state)+evalC [":?"] state = (Right (help state), state)+evalC [":show", v] state = (evalShow v state, state)+evalC [":delete", v] state = (Right ("deleted var: " ++ v), delete v state)+evalC [":units", u1, u2] state = evalUnits [u1, u2] state+evalC [":units", u] state = evalUnits [u] state+evalC [":units"] state = showUnits state+evalC [":reset"] _ = (Right "REPL reset ", emptyState)+evalC c state = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), state)++evalShow :: String -> State -> Either String String+evalShow n state =+ maybe (Left ("Unbound variable: " ++ n)) (\v -> Right (showVar n v state)) (lookup n state)++evalUnits :: [String] -> State -> (Either String String, State)+evalUnits us s = showUnits (setUnits us s)++showUnits :: State -> (Either String String, State)+showUnits s = (Right ("Units:\n length = " ++ lengthUnit s ++ "\n speed = " ++ speedUnit s), s)++help :: State -> String+help s =+ "\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} show {var}\n" +++ " :delete {var} delete {var}\n" +++ " :units length speed set length and speed units used for display\n" +++ " see supported length and speed format\n" +++ " currently: length = " +++ lengthUnit s +++ "; speed = " +++ speedUnit s +++ "\n" +++ " :units show length and speed units used for display\n" +++ " :reset reset REPL to default state (including deleting all variables)\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: " +++ showLength r84 s +++ "\n" +++ "\n alongTrackDistance pos ga signed distance of pos to perpendicular of great arc ga\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: " +++ showLength r84 s +++ "\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" +++ " intercept track pos spd time needed by interceptor at pos and travelling at spd to intercept target\n" +++ " intercept 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" +++ " greatArc pos1 pos2 great arc passing by pos1 and pos2\n" +++ " greatArc track dur great arc from track and duration\n" +++ " greatCircle pos1 pos2 great circle passing by pos1 and pos2\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 ga great circle from great arc\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 -> State -> State+save (Right v) k state = insert k v state+save _ _ state = state++showR :: Result -> State -> Either String String+showR (Left err) _ = Left err+showR (Right v) s = Right (showV v s)++showVar :: String -> Value -> State -> String+showVar n v s = n ++ "=" ++ showV v s
app/Show.hs view
@@ -8,94 +8,100 @@ -- -- Show 'Value's. -- -module Show - ( showV - ) where - -import Data.Geo.Jord -import Data.List (intercalate) -import State - +module Show+ ( showV+ ) where++import Data.Geo.Jord+import Data.List (intercalate)+import State+ -- | show value. -showV :: Value -> State -> String -showV (Ang a) _ = "angle: " ++ showAng a -showV (Bool b) _ = show b -showV (Cpa c) s = - "closest point of approach:" ++ - "\n time : " ++ - show (cpaTime c) ++ - "\n distance: " ++ - showLength (cpaDistance c) s ++ - "\n pos1 : " ++ - showLl (fromNVector . cpaPosition1 $ c :: LatLong) ++ - "\n pos2 : " ++ showLl (fromNVector . cpaPosition2 $ c :: LatLong) -showV (Dlt d) s = - "Delta:" ++ - "\n x: " ++ - showLength (dx d) s ++ - "\n y: " ++ showLength (dy d) s ++ "\n z: " ++ showLength (dz d) s -showV (Dur d) _ = "duration: " ++ show d -showV (Double d) _ = show d -showV (Em m) _ = "Earth model: " ++ show m -showV (Ep p) s = - "ECEF:" ++ - "\n x: " ++ - showLength (ex p) s ++ - "\n y: " ++ showLength (ey p) s ++ "\n z: " ++ showLength (ez p) s -showV (FrmB y p r) _ = - "Body (vehicle) frame:" ++ - "\n yaw : " ++ - showAng y ++ "\n pitch: " ++ showAng p ++ "\n roll : " ++ showAng r -showV (FrmL w) _ = "Local frame:" ++ "\n wander azimuth: " ++ showAng w -showV FrmN _ = "North, East, Down frame" -showV (Gc gc) _ = "great circle: " ++ show gc -showV (Gp g) s = "latlong: " ++ showLl ll ++ "\n height : " ++ showLength h s - where - ll = pos g - h = height g -showV (Intp i) s = - "intercept:" ++ - "\n time : " ++ - show (interceptTime i) ++ - "\n distance : " ++ - showLength (interceptDistance i) s ++ - "\n pos : " ++ - showLl (fromNVector . interceptPosition $ i :: LatLong) ++ - "\n interceptor speed : " ++ - showSpeed (interceptorSpeed i) s ++ - "\n interceptor bearing: " ++ showAng (interceptorBearing i) -showV (Len l) s = "length: " ++ showLength l s -showV (Ned d) s = - "NED:" ++ - "\n north: " ++ - showLength (north d) s ++ - "\n east : " ++ showLength (east d) s ++ "\n down : " ++ showLength (down d) s -showV (Np nv) s = - "n-vector: " ++ - show x ++ ", " ++ show y ++ ", " ++ show z ++ "\n height : " ++ showLength h s - where - v = vec (pos nv) - x = vx v - y = vy v - z = vz v - h = height nv -showV (Trk t) s = - "track:" ++ - "\n position: " ++ - showLl (fromNVector . trackPos $ t :: LatLong) ++ - "\n height : " ++ - showLength (height . trackPos $ t) s ++ - "\n bearing : " ++ - showAng (trackBearing t) ++ "\n speed : " ++ showSpeed (trackSpeed t) s -showV (Spd spd) s = "speed: " ++ showSpeed spd s -showV (Vals []) _ = "empty" -showV (Vals vs) s = "\n " ++ intercalate "\n\n " (map (`showV` s) 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)) ++ ")" +showV :: Value -> State -> String+showV (Ang a) _ = "angle: " ++ showAng a+showV (Bool b) _ = show b+showV (Cpa c) s =+ "closest point of approach:" +++ "\n time : " +++ show (cpaTime c) +++ "\n distance: " +++ showLength (cpaDistance c) s +++ "\n pos1 : " +++ showLl (fromNVector . cpaPosition1 $ c :: LatLong) +++ "\n pos2 : " ++ showLl (fromNVector . cpaPosition2 $ c :: LatLong)+showV (Dlt d) s =+ "Delta:" +++ "\n x: " +++ showLength (dx d) s +++ "\n y: " ++ showLength (dy d) s ++ "\n z: " ++ showLength (dz d) s+showV (Dur d) _ = "duration: " ++ show d+showV (Double d) _ = show d+showV (Em m) _ = "Earth model: " ++ show m+showV (Ep p) s =+ "ECEF:" +++ "\n x: " +++ showLength (ex p) s +++ "\n y: " ++ showLength (ey p) s ++ "\n z: " ++ showLength (ez p) s+showV (FrmB y p r) _ =+ "Body (vehicle) frame:" +++ "\n yaw : " +++ showAng y ++ "\n pitch: " ++ showAng p ++ "\n roll : " ++ showAng r+showV (FrmL w) _ = "Local frame:" ++ "\n wander azimuth: " ++ showAng w+showV FrmN _ = "North, East, Down frame"+showV (Ga ga) _ =+ "great arc: passing by " +++ (showLl . nvectorToLatLong . gaStart $ ga) +++ " and " ++ (showLl . nvectorToLatLong . gaEnd $ ga)+showV (Gc gc) _ =+ "great circle: passing by " +++ (showLl . nvectorToLatLong . gcPos $ gc) ++ " heading on " ++ (showAng . gcBearing $ gc)+showV (Gp g) s = "latlong: " ++ showLl ll ++ "\n height : " ++ showLength h s+ where+ ll = pos g+ h = height g+showV (Intp i) s =+ "intercept:" +++ "\n time : " +++ show (interceptTime i) +++ "\n distance : " +++ showLength (interceptDistance i) s +++ "\n pos : " +++ showLl (fromNVector . interceptPosition $ i :: LatLong) +++ "\n interceptor speed : " +++ showSpeed (interceptorSpeed i) s +++ "\n interceptor bearing: " ++ showAng (interceptorBearing i)+showV (Len l) s = "length: " ++ showLength l s+showV (Ned d) s =+ "NED:" +++ "\n north: " +++ showLength (north d) s +++ "\n east : " ++ showLength (east d) s ++ "\n down : " ++ showLength (down d) s+showV (Np nv) s =+ "n-vector: " +++ show x ++ ", " ++ show y ++ ", " ++ show z ++ "\n height : " ++ showLength h s+ where+ v = vec (pos nv)+ x = vx v+ y = vy v+ z = vz v+ h = height nv+showV (Trk t) s =+ "track:" +++ "\n position: " +++ showLl (fromNVector . trackPos $ t :: LatLong) +++ "\n height : " +++ showLength (height . trackPos $ t) s +++ "\n bearing : " +++ showAng (trackBearing t) ++ "\n speed : " ++ showSpeed (trackSpeed t) s+showV (Spd spd) s = "speed: " ++ showSpeed spd s+showV (Vals []) _ = "empty"+showV (Vals vs) s = "\n " ++ intercalate "\n\n " (map (`showV` s) 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)) ++ ")"
app/State.hs view
@@ -49,6 +49,7 @@ Angle -- ^ yaw, pitch and roll of Body frame | FrmL Angle -- ^ wander azimuth of Local frame | FrmN -- ^ North, east, down frame + | Ga GreatArc -- ^ great arc | Gc GreatCircle -- ^ great circle | Gp (AngularPosition LatLong) -- ^ latitude, longitude and height | Intp (Intercept (AngularPosition NVector)) -- ^ Intercept
+ benchmarks/GeodeticsBG.hs view
@@ -0,0 +1,37 @@+module GeodeticsBG + ( bggeodetics + ) where + +import Criterion.Types +import Data.Geo.Jord + +bggeodetics :: Benchmark +bggeodetics = + bgroup + "Geodetics" + [ bench "angularDistance" $ whnf (angularDistance nv1 nv2) (Just nv3) + , bench "crossTrackDistance" $ whnf (crossTrackDistance84 nv3) gc + , bench "destination" $ whnf (destination84 nv1 a) l + , bench "finalBearing" $ whnf (finalBearing nv1) nv2 + , bench "initialBearing" $ whnf (initialBearing nv1) nv2 + , bench "interpolate" $ whnf (interpolate nv1 nv2) 0.5 + , bench "surfaceDistance" $ whnf (surfaceDistance84 nv1) nv2 + ] + +nv1 :: NVector +nv1 = nvector 0.5504083453140064 0.12711022980808237 0.8251627978083076 + +nv2 :: NVector +nv2 = nvector 0.484947835927087 0.1582112780286092 0.860113241343365 + +nv3 :: NVector +nv3 = nvector 0.5225962210695282 0.11083913756305296 0.8453448262739457 + +gc :: GreatCircle +gc = greatCircle (nv1, nv2) + +a :: Angle +a = decimalDegrees 45.0 + +l :: Length +l = kilometres 5000
+ benchmarks/KinematicsBG.hs view
@@ -0,0 +1,45 @@+module KinematicsBG + ( bgkinematics + ) where + +import Criterion.Types +import Data.Geo.Jord + +bgkinematics :: Benchmark +bgkinematics = + bgroup + "Kinematics" + [ bgroup + "CPA" + [ bench "in the past" $ whnf (cpa84 t1) t2 + , bench "in the future" $ whnf (cpa84 t3) t4 + , bench "same positions" $ whnf (cpa84 t1') t1 + ] + , bgroup + "intercept" + [ bench "min speed" $ whnf (intercept84 t5) ip1 + , bench "by speed" $ whnf (interceptBySpeed84 t5 ip1) (knots 700) + , bench "by time" $ whnf (interceptByTime84 t5 ip1) (seconds 2700) + ] + ] + +t1 :: Track NVector +t1 = Track (latLongToNVector (decimalLatLong 20 (-60))) (decimalDegrees 10) (knots 15) + +t1' :: Track NVector +t1' = Track (latLongToNVector (decimalLatLong 20 (-60))) (decimalDegrees 10) (knots 15) + +t2 :: Track NVector +t2 = Track (latLongToNVector (decimalLatLong 34 (-50))) (decimalDegrees 220) (knots 300) + +t3 :: Track NVector +t3 = Track (latLongToNVector (decimalLatLong 30 30)) (decimalDegrees 45) (knots 400) + +t4 :: Track NVector +t4 = Track (latLongToNVector (decimalLatLong 30.01 30)) (decimalDegrees 315) (knots 400) + +t5 :: Track NVector +t5 = Track (latLongToNVector (decimalLatLong 34 (-50))) (decimalDegrees 220) (knots 600) + +ip1 :: NVector +ip1 = latLongToNVector (decimalLatLong 20 (-60))
+ benchmarks/Main.hs view
@@ -0,0 +1,9 @@+module Main where + +import Criterion.Main +import GeodeticsBG +import KinematicsBG +import TransformationBG + +main :: IO () +main = defaultMain [bggeodetics, bgkinematics, bgtransformation]
+ benchmarks/TransformationBG.hs view
@@ -0,0 +1,36 @@+module TransformationBG + ( bgtransformation + ) where + +import Criterion.Types +import Data.Geo.Jord + +bgtransformation :: Benchmark +bgtransformation = + bgroup + "Transformation" + [ bench "latLongToNVector" $ whnf latLongToNVector ll + , bench "nvectorToLatLong" $ whnf nvectorToLatLong nv + , bgroup + "Ellipsoidal" + [ bench "ecefToNVector" $ whnf (`ecefToNVector` wgs84) ep + , bench "nvectorToEcef" $ whnf (`nvectorToEcef` wgs84) ap + ] + , bgroup + "Spherical" + [ bench "ecefToNVector" $ whnf (`ecefToNVector` s84) ep + , bench "nvectorToEcef" $ whnf (`nvectorToEcef` s84) ap + ] + ] + +ll :: LatLong +ll = decimalLatLong 55.6050 13.0038 + +ap :: AngularPosition NVector +ap = AngularPosition nv (metres 15000.0) + +nv :: NVector +nv = nvector 0.5 0.5 0.7071 + +ep :: EcefPosition +ep = ecefMetres 5733855.7748 (-6370998.3802) 7008137.5108
jord.cabal view
@@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 8bc0043cf9b7f8b4b203c71df609cfa6cc47735768a22d31029d559d8f49ad70+-- hash: 2f1ced38bca95f83cdbb2509dc23645a34c628652e960e7a84e9e1b7b6d4ca8c name: jord-version: 0.4.2.0+version: 0.5.0.0 synopsis: Geographical Position Calculations description: Please see the README on GitHub at <https://github.com/ofmooseandmen/jord#readme> category: Geography@@ -47,12 +47,28 @@ Data.Geo.Jord.Transformation Data.Geo.Jord.Vector3d other-modules:- Data.Geo.Jord.Parse+ Data.Geo.Jord.Internal+ Data.Geo.Jord.Parser hs-source-dirs: src ghc-options: -Wall build-depends: base >=4.9 && <5+ default-language: Haskell2010++executable jord-benchmarks+ main-is: Main.hs+ other-modules:+ GeodeticsBG+ KinematicsBG+ TransformationBG+ hs-source-dirs:+ benchmarks+ ghc-options: -Wall+ build-depends:+ base >=4.9 && <5+ , criterion+ , jord default-language: Haskell2010 executable jord-repl
src/Data/Geo/Jord.hs view
@@ -55,4 +55,4 @@ -- | version. jordVersion :: String -jordVersion = "0.4.2.0" +jordVersion = "0.5.0.0"
src/Data/Geo/Jord/Angle.hs view
@@ -17,6 +17,7 @@ , dms , dmsE , dmsF + , radians -- * Calculations , arcLength , central @@ -34,7 +35,9 @@ , getMinutes , getSeconds , getMilliseconds + -- * Conversions , toDecimalDegrees + , toRadians -- * Read , angle , readAngle @@ -46,7 +49,7 @@ import Control.Monad.Fail import Data.Fixed import Data.Geo.Jord.Length -import Data.Geo.Jord.Parse +import Data.Geo.Jord.Parser import Data.Geo.Jord.Quantity import Data.Maybe import Prelude hiding (fail, length) @@ -91,7 +94,7 @@ decimalDegrees :: Double -> Angle decimalDegrees dec = Angle (round (dec * 3600000.0)) --- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. +-- | 'Angle' from the 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. @@ -103,7 +106,7 @@ 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. +-- | 'Angle' from the 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. @@ -121,7 +124,7 @@ (fromIntegral millis / 3600000.0 :: Double)) (signum degs) --- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. +-- | 'Angle' from the 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. @@ -133,13 +136,17 @@ where e = dmsE degs mins secs millis +-- | 'Angle' from the given radians. +radians :: Double -> Angle +radians r = decimalDegrees (r / pi * 180.0) + -- | @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) -- | @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 s r = radians (toMetres s / toMetres r) -- | Returns the given 'Angle' negated. negate' :: Angle -> Angle @@ -161,11 +168,11 @@ -- | @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' y x = radians (atan2 y x) -- | @asin' a@ computes arc sine of @a@. asin' :: Double -> Angle -asin' a = fromRadians (asin a) +asin' a = radians (asin a) -- | @cos' a@ returns the cosinus of @a@. cos' :: Angle -> Double @@ -174,10 +181,6 @@ -- | @sin' a@ returns the sinus of @a@. sin' :: Angle -> Double sin' a = sin (toRadians a) - --- | radians to degrees. -fromRadians :: Double -> Angle -fromRadians r = decimalDegrees (r / pi * 180.0) -- | degrees to radians. toRadians :: Angle -> Double
src/Data/Geo/Jord/Duration.hs view
@@ -19,7 +19,7 @@ , minutes , seconds , hms - -- * Accessors + -- * Conversions , toHours , toMinutes , toSeconds @@ -30,7 +30,7 @@ ) where import Control.Monad.Fail -import Data.Geo.Jord.Parse +import Data.Geo.Jord.Parser import Data.Geo.Jord.Quantity import Prelude hiding (fail) import Text.ParserCombinators.ReadP
src/Data/Geo/Jord/Geodetics.hs view
@@ -18,7 +18,16 @@ -- * The 'GreatCircle' type GreatCircle , IsGreatCircle(..)+ , gcPos+ , gcBearing+ -- * The 'GreatArc' type+ , GreatArc+ , IsGreatArc(..)+ , gaStart+ , gaEnd -- * Calculations + , alongTrackDistance+ , alongTrackDistance84 , angularDistance , antipode , crossTrackDistance@@ -28,8 +37,10 @@ , finalBearing , initialBearing , interpolate+ , intersection , intersections- , insideSurface+ , isBetween+ , isInsideSurface , mean , surfaceDistance , surfaceDistance84@@ -40,14 +51,14 @@ 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.Internal (nvec, sad) 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 Data.Maybe (fromJust, fromMaybe, isNothing) import Prelude hiding (fail) -- | A circle on the __surface__ of the Earth which lies in a plane passing through @@ -57,22 +68,17 @@ -- It is internally represented as its normal vector - i.e. the normal vector -- to the plane containing the great circle. -- --- See 'greatCircle', 'greatCircleE', 'greatCircleF' or 'greatCircleBearing' constructors. --- +-- see 'IsGreatCircle'. data GreatCircle = GreatCircle- { normal :: Vector3d- , dscr :: String- } deriving (Eq)--instance Show GreatCircle where- show = dscr+ { gcNormal :: Vector3d -- ^ normal vector to the plane containing the great circle+ , gcPos :: NVector -- ^ position (/n/-vector) on the great circle+ , gcBearing :: Angle -- ^ bearing from 'gcPos'.+ } deriving (Eq, Show) -- | Class for data from which a 'GreatCircle' can be computed. -class (Show a) =>- IsGreatCircle a- where+class 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)+ greatCircle a = fromMaybe (error "Could not make 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 =@@ -82,6 +88,29 @@ where e = greatCircleE a +-- | A closed segment of 'GreatCircle'. It represent the shortest path on the __surface__ of the Earth+-- between the two positions.+--+-- see 'IsGreatArc'.+data GreatArc = GreatArc+ { gaNormal :: Vector3d -- ^ normal vector to the plane containing the great circle+ , gaStart :: NVector -- ^ start position (/n/-vector) of the great arc+ , gaEnd :: NVector -- ^ end position (/n/-vector) of the great arc+ } deriving (Eq, Show)++-- | Class for data from which a 'GreatArc' can be computed.+class IsGreatArc a where+ greatArc :: a -> GreatArc -- ^ 'GreatCircle' from @a@, if 'greatArcE' returns a 'Left', this function 'error's.+ greatArc a = fromMaybe (error "Could not make a Great Arc") (greatArcF a)+ greatArcE :: a -> Either String GreatArc -- ^ 'GreatArc' from @a@, A 'Left' indicates an error.+ greatArcF :: (MonadFail m) => a -> m GreatArc -- ^ 'GreatArc' from @a@, if 'greatArcE' returns a 'Left', this function 'fail's.+ greatArcF a =+ case e of+ Left err -> fail err+ Right ga -> return ga+ where+ e = greatArcE a+ -- | 'GreatCircle' passing by both given positions'. A 'Left' indicates that given positions are -- equal or antipodal. -- @@ -90,51 +119,93 @@ -- 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+instance NTransform 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)))+ | isNothing b = Left "Invalid Great Circle: positions are equal"+ | otherwise = Right (GreatCircle (vcross v1 v2) nv1 (fromJust b)) where- v1 = vector3d p1- v2 = vector3d p2+ nv1 = pos . toNVector $ p1+ v1 = vec nv1+ nv2 = pos . toNVector $ p2+ v2 = vec nv2+ b = initialBearing nv1 nv2 -- | 'GreatCircle' passing by the given position and heading on given bearing. -- -- @ -- 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))+instance NTransform a => IsGreatCircle (a, Angle) where+ greatCircleE (p, b) = Right (GreatCircle (vsub n' e') nv b) where- v = vector3d p+ nv = pos . toNVector $ p+ v = nvec nv 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' from given 'GreatArc'.+instance IsGreatCircle GreatArc where+ greatCircleE (GreatArc n s e) =+ case initialBearing s e of+ Nothing -> Left "Could not computed initial bearing"+ (Just b) -> Right (GreatCircle n s b)++-- | 'GreatArc' passing by both given positions'. A 'Left' indicates that given positions are+-- equal or antipodal.+--+-- @+-- let p1 = decimalLatLongHeight 45.0 (-143.5) (metres 1500)+-- let p2 = decimalLatLongHeight 46.0 14.5 (metres 3000)+-- greatArc (p1, p2) -- heights are ignored, great arc are always at earth surface.+-- @+instance NTransform a => IsGreatArc (a, a) where+ greatArcE ps@(p1, p2) =+ case greatCircleE ps of+ Left e -> Left e+ Right gcv ->+ Right (GreatArc (gcNormal gcv) (pos . toNVector $ p1) (pos . toNVector $ p2))++-- | @alongTrackDistance p ga r@ how far position @p@ is along a path described+-- by great arc @ga@: if a perpendicular is drawn from @p@ to the great arc, the+-- along-track distance is the signed distance from the start point to where the+-- perpendicular crosses the path.+--+-- @+-- let p = decimalLatLong 53.2611 (-0.7972)+-- let ga = greatArc (decimalLatLong 53.3206 (-1.7297)) (decimalLatLong 53.1887 0.1334)+-- alongTrackDistance p ga r84 -- 62.3315757 kilometres+-- @+alongTrackDistance :: (NTransform a) => a -> GreatArc -> Length -> Length+alongTrackDistance p (GreatArc n s _) =+ arcLength (sad' (nvec s) (vcross (vcross n (nvec p)) n) (Just n))++-- | 'alongTrackDistance' using the mean radius of the WGS84 reference ellipsoid.+alongTrackDistance84 :: (NTransform a) => a -> GreatArc -> Length+alongTrackDistance84 p ga = alongTrackDistance p ga r84+ -- | @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+angularDistance p1 p2 n = sad' v1 v2 vn where- v1 = vector3d p1- v2 = vector3d p2- vn = fmap vector3d n+ v1 = nvec p1+ v2 = nvec p2+ vn = fmap nvec 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)+antipode p = fromNVector (angular (vscale (nvec 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@. +-- | @crossTrackDistance p gc r@ 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: @@ -142,15 +213,15 @@ -- @ -- 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 r84 = (- crossTrackDistance p gc2 r84) -- -- 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 p (GreatCircle n _ _) =+ arcLength (sub (sad' n (nvec p) Nothing) (decimalDegrees 90)) -- | 'crossTrackDistance' using the mean radius of the WGS84 reference ellipsoid. crossTrackDistance84 :: (NTransform a) => a -> GreatCircle -> Length@@ -201,10 +272,10 @@ 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))+ | otherwise = Just (normalise (sad' gc1 gc2 (Just v1)) (decimalDegrees 360)) where- v1 = vector3d p1- v2 = vector3d p2+ v1 = nvec p1+ v2 = nvec p2 gc1 = vcross v1 v2 -- great circle through p1 & p2 gc2 = vcross v1 (vec northPole) -- great circle through p1 & north pole @@ -238,41 +309,27 @@ 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@). --- --- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles --- will sum to less than 360°, and exterior point angles will be small but non-zero. --- --- Always returns 'False' if @ps@ does not at least defines a triangle. --- --- @ --- 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 intersection between the two given 'GreatArc's.+--+-- see also 'intersections'+--+-- @+-- let spd = kilometresPerHour 1000+-- let t1 = Track (decimalLatLong 51.885 0.235) (decimalDegrees 108.63) spd+-- let t2 = Track (decimalLatLong 49.008 2.549) (decimalDegrees 32.72) spd+-- let oneHour = hours 1+-- let ga1 = greatArc (t1, oneHour)+-- let ga2 = greatArc (t2, oneHour)+-- intersection ga1 ga2 = Just (decimalLatLong 50.9017225 4.494278333333333)+-- @+intersection :: (NTransform a) => GreatArc -> GreatArc -> Maybe a+intersection ga1@(GreatArc n1 _ _) ga2@(GreatArc n2 _ _) =+ case intersections' n1 n2 of+ Nothing -> Nothing+ (Just (i1, i2))+ | isBetween i1 ga1 && isBetween i1 ga2 -> Just i1+ | isBetween i2 ga1 && isBetween i2 ga2 -> Just i2+ | otherwise -> Nothing -- | Computes the intersections between the two given 'GreatCircle's. -- Two 'GreatCircle's intersect exactly twice unless there are equal (regardless of orientation), @@ -286,13 +343,64 @@ -- 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)+intersections (GreatCircle n1 _ _) (GreatCircle n2 _ _) = intersections' n1 n2++-- | @isBetween p ga@ determines whether position @p@ is between start and end points+-- of great arc @ga@.+-- If @p@ is not on the great arc, returns whether @p@ is within the area bound+-- by perpendiculars to the great arc at each point (in the same hemisphere).+--+isBetween :: (NTransform a) => a -> GreatArc -> Bool+isBetween p (GreatArc _ s e) = between && hemisphere where- i = vcross (normal gc1) (normal gc2)+ v0 = nvec p+ v1 = nvec s+ v2 = nvec e+ v10 = vsub v0 v1+ v12 = vsub v2 v1+ v20 = vsub v0 v2+ v21 = vsub v1 v2+ e1 = vdot v10 v12 -- p is on e side of s+ e2 = vdot v20 v21 -- p is on s side of e+ between = e1 >= 0 && e2 >= 0+ hemisphere = vdot v0 v1 >= 0 && vdot v0 v2 >= 0 +-- | @isInsideSurface 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@).+--+-- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles+-- will sum to less than 360°, and exterior point angles will be small but non-zero.+--+-- Always returns 'False' if @ps@ does not at least defines a triangle.+--+-- @+-- 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+-- isInsideSurface hoor polygon = True+-- isInsideSurface hassleholm polygon = False+-- @+isInsideSurface :: (Eq a, NTransform a) => a -> [a] -> Bool+isInsideSurface p ps+ | null ps = False+ | head ps == last ps = isInsideSurface p (init ps)+ | length ps < 3 = False+ | otherwise =+ let aSum =+ foldl+ (\a v' -> add a (uncurry sad' v' (Just v)))+ (decimalDegrees 0)+ (egdes (map (vsub v) vs))+ in abs (toDecimalDegrees aSum) > 180.0+ where+ v = nvec p+ vs = fmap nvec ps+ -- | @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 @@ -314,7 +422,7 @@ then Just (fromNVector (angular (vunit (foldl vadd vzero vs)) zero)) else Nothing where- vs = fmap vector3d ps+ vs = fmap nvec ps ts = filter (\l -> length l == 2) (subsequences vs) antipodals = filter (\t -> (realToFrac (vnorm (vadd (head t) (last t)) :: Double) :: Nano) == 0) ts@@ -327,15 +435,9 @@ 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+-- | Signed angular distance - see 'sad'.+sad' :: Vector3d -> Vector3d -> Maybe Vector3d -> Angle+sad' v1 v2 n = radians (sad v1 v2 n) -- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] egdes :: [Vector3d] -> [(Vector3d, Vector3d)]@@ -348,11 +450,13 @@ 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+intersections' :: (NTransform a) => Vector3d -> Vector3d -> Maybe (a, a)+intersections' n1 n2+ | (vnorm i :: Double) == 0.0 = Nothing+ | otherwise+ , let ni = fromNVector (angular (vunit i) zero) = Just (ni, antipode ni)+ where+ i = vcross n1 n2
+ src/Data/Geo/Jord/Internal.hs view
@@ -0,0 +1,37 @@+-- |+-- Module: Data.Geo.Jord.Internal+-- Copyright: (c) 2018 Cedric Liegeois+-- License: BSD3+-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr>+-- Stability: experimental+-- Portability: portable+--+-- internal functions.+--+module Data.Geo.Jord.Internal+ ( ad+ , nvec+ , sad+ ) where++import Data.Geo.Jord.AngularPosition (pos)+import Data.Geo.Jord.Transformation (NTransform(..))+import Data.Geo.Jord.Vector3d++-- | angle in __radians__ between 2 /n/-vectors (as 'Vector3d').+ad :: Vector3d -> Vector3d -> Double+ad v1 v2 = sad v1 v2 Nothing++-- | /n/-vector (as a 'Vector3d') from given position.+nvec :: (NTransform a) => a -> Vector3d+nvec = vec . pos . toNVector++-- | Signed angle in __radians__ between 2 /n/-vectors (as 'Vector3d').+-- If @n@ is 'Nothing', the angle is always in [0..pi], otherwise it is in [-pi, +pi],+-- signed + if @v1@ is clockwise looking along @n@, - in opposite direction.+sad :: Vector3d -> Vector3d -> Maybe Vector3d -> Double+sad 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
src/Data/Geo/Jord/Kinematics.hs view
@@ -1,37 +1,39 @@--- | --- 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> --- +{-# LANGUAGE FlexibleInstances #-}++-- |+-- 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. + -- * The 'Course' type. , Course- -- * The 'Cpa' type. + -- * The 'Cpa' type. , Cpa , cpaTime , cpaDistance , cpaPosition1 , cpaPosition2- -- * The 'Intercept' type. + -- * The 'Intercept' type. , Intercept , interceptTime , interceptDistance , interceptPosition , interceptorBearing , interceptorSpeed- -- * Calculations + -- * Calculations , course , position , position84@@ -51,6 +53,7 @@ import Data.Geo.Jord.Duration import Data.Geo.Jord.Earth import Data.Geo.Jord.Geodetics+import Data.Geo.Jord.Internal(ad, nvec) import Data.Geo.Jord.LatLong import Data.Geo.Jord.Length import Data.Geo.Jord.NVector@@ -60,18 +63,26 @@ import Data.Geo.Jord.Vector3d import Data.Maybe (isNothing) --- | 'Track' represents the state of a vehicle by its current position, bearing and speed. +-- | '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. + { 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+-- | 'GreatCircle' from track.+instance NTransform 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. +-- | 'GreatArc' from track and duration using the mean radius of the WGS84 reference ellipsoid.+instance NTransform a => IsGreatArc (Track a, Duration) where+ greatArcE (t, d) = greatArcE (t, d, r84)++-- | 'GreatArc' from track, duration and earth mean radius.+instance NTransform a => IsGreatArc (Track a, Duration, Length) where+ greatArcE (t, d, r) = greatArcE (trackPos t, position t d r)++-- | 'Course' represents the cardinal direction in which the vehicle is to be steered. newtype Course = Course Vector3d deriving (Eq, Show)@@ -79,24 +90,24 @@ 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. +-- | 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. + { 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. +-- | 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. + { 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 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@@ -105,37 +116,37 @@ 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 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. +-- | '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 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)@@ -149,42 +160,42 @@ cp2 = position' p2 s2 c2 t r d = surfaceDistance cp1 cp2 r --- | 'cpa' using the mean radius of the WGS84 reference ellipsoid. +-- | '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 "behind" the target --- --- @ --- 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.633367756059) --- fmap interceptTime i = Just (seconds 5993.831) --- @ +-- | @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 "behind" the target+--+-- @+-- 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.633367756059)+-- fmap interceptTime i = Just (seconds 5993.831)+-- @ intercept :: (Eq a, NTransform a) => Track a -> a -> Length -> Maybe (Intercept a) intercept t p r = interceptByTime t p (seconds (timeToIntercept t p r)) r --- | 'intercept' using the mean radius of the WGS84 reference ellipsoid. +-- | '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 returned by 'intercept' +-- | @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 returned by 'intercept' interceptBySpeed :: (Eq a, NTransform a) => Track a -> a -> Speed -> Length -> Maybe (Intercept a) interceptBySpeed t p s r | isNothing minInt = Nothing@@ -193,29 +204,29 @@ where minInt = intercept t p r --- | 'interceptBySpeed' using the mean radius of the WGS84 reference ellipsoid. +-- | '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) --- @ --- --- Note: contrary to 'intercept' and 'interceptBySpeed' this function handles --- cases where the interceptor has to catch up the target. +-- | @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)+-- @+--+-- Note: contrary to 'intercept' and 'interceptBySpeed' this function handles+-- cases where the interceptor has to catch up the target. interceptByTime :: (Eq a, NTransform a) => Track a -> a -> Duration -> Length -> Maybe (Intercept a) interceptByTime t p d r | toMilliseconds d <= 0 = Nothing@@ -227,11 +238,11 @@ ib = initialBearing p ipos <|> initialBearing p (trackPos t) is = metresPerSecond (toMetres idist / toSeconds d) --- | 'interceptByTime' using the mean radius of the WGS84 reference ellipsoid. +-- | '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 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@@ -240,47 +251,47 @@ h0 = height nv0 v1 = position'' v0 s (vec c) sec r --- | position from speed course and seconds. +-- | position from speed course and seconds. position'' :: Vector3d -> Speed -> Vector3d -> Double -> Length -> Vector3d position'' v0 s c sec r = v1 where a = toMetresPerSecond s / toMetres r * sec v1 = vadd (vscale v0 (cos a)) (vscale c (sin a)) --- | time to CPA. +-- | 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 = vec3d p1+ v10 = nvec p1 c10 = vec c1 rm = toMetres r w1 = toMetresPerSecond s1 / rm- v20 = vec . pos . toNVector $ p2+ v20 = nvec p2 c20 = vec c2 w2 = toMetresPerSecond s2 / rm --- | time to intercept with minimum speed +-- | time to intercept with minimum speed timeToIntercept :: (NTransform a) => Track a -> a -> Length -> Double timeToIntercept (Track p2 b2 s2) p1 r = intMinNrRec v10v20 v10c2 w2 (sep v10 v20 c2 s2 r) t0 0 where- v10 = vec3d p1- v20 = vec3d p2+ v10 = nvec p1+ v20 = nvec p2 c2 = vec (course p2 b2) v10v20 = vdot v10 v20 v10c2 = vdot v10 c2 s2mps = toMetresPerSecond s2 rm = toMetres r w2 = s2mps / rm- s0 = ad v10 v20 -- initial angular distance between target and interceptor - t0 = rm * s0 / s2mps -- assume target is travelling towards interceptor + s0 = ad v10 v20 -- initial angular distance between target and interceptor+ t0 = rm * s0 / s2mps -- assume target is travelling towards interceptor --- | time to intercept with speed. +-- | time to intercept with speed. timeToInterceptSpeed :: (NTransform a) => Track a -> a -> Speed -> Length -> Double timeToInterceptSpeed (Track p2 b2 s2) p1 s1 r = intSpdNrRec v10v20 v10c2 w1 w2 (sep v10 v20 c2 s2 r) t0 0 where- v10 = vec3d p1- v20 = vec3d p2+ v10 = nvec p1+ v20 = nvec p2 c2 = vec (course p2 b2) v10v20 = vdot v10 v20 v10c2 = vdot v10 c2@@ -353,21 +364,21 @@ c = cpaC v10 c10 w1 v20 c20 w2 d = cpaD v10 c10 w1 v20 c20 w2 --- | Newton-Raphson for CPA time. +-- | Newton-Raphson for CPA time. 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 + | i == 50 = -1.0 -- no convergence | abs fi < 1e-11 = 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. +-- | Newton-Raphson for min speed intercept. intMinNrRec :: Double -> Double -> Double -> (Double -> Double) -> Double -> Int -> Double intMinNrRec v10v20 v10c2 w2 st ti i- | i == 50 = -1.0 -- no convergence + | i == 50 = -1.0 -- no convergence | abs fi < 1e-11 = ti1 | otherwise = intMinNrRec v10v20 v10c2 w2 st ti1 (i + 1) where@@ -385,10 +396,10 @@ fi = f / df ti1 = ti - fi --- | Newton-Raphson for speed intercept. +-- | Newton-Raphson for speed intercept. intSpdNrRec :: Double -> Double -> Double -> Double -> (Double -> Double) -> Double -> Int -> Double intSpdNrRec v10v20 v10c2 w1 w2 st ti i- | i == 50 = -1.0 -- no convergence + | i == 50 = -1.0 -- no convergence | abs fi < 1e-11 = ti1 | otherwise = intSpdNrRec v10v20 v10c2 w1 w2 st ti1 (i + 1) where@@ -400,16 +411,8 @@ df = (dsdt - (si / ti)) / ti fi = f / df ti1 = ti - fi- --- | angular separation in radians at ti between v10 and track with initial position v20, --- course c2 and speed s2. ++-- | angular separation in radians at ti between v10 and track with initial position v20,+-- course c2 and speed s2. sep :: Vector3d -> Vector3d -> Vector3d -> Speed -> Length -> Double -> Double sep v10 v20 c2 s2 r ti = ad v10 (position'' v20 s2 c2 ti r)---- | 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)--vec3d :: (NTransform a) => a -> Vector3d-vec3d = vec . pos . toNVector
src/Data/Geo/Jord/LatLong.hs view
@@ -33,7 +33,7 @@ import Control.Monad.Fail import Data.Char import Data.Geo.Jord.Angle -import Data.Geo.Jord.Parse +import Data.Geo.Jord.Parser import Data.Maybe import Prelude hiding (fail) import Text.ParserCombinators.ReadP
src/Data/Geo/Jord/Length.hs view
@@ -30,7 +30,7 @@ import Control.Applicative import Control.Monad.Fail -import Data.Geo.Jord.Parse +import Data.Geo.Jord.Parser import Data.Geo.Jord.Quantity import Prelude hiding (fail, length) import Text.ParserCombinators.ReadP
− src/Data/Geo/Jord/Parse.hs
@@ -1,52 +0,0 @@--- | --- 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/Parser.hs view
@@ -0,0 +1,52 @@+-- | +-- Module: Data.Geo.Jord.Parser +-- 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.Parser + ( 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/Speed.hs view
@@ -32,7 +32,7 @@ import Control.Applicative import Control.Monad.Fail -import Data.Geo.Jord.Parse +import Data.Geo.Jord.Parser import Data.Geo.Jord.Quantity import Prelude hiding (fail) import Text.ParserCombinators.ReadP
src/Data/Geo/Jord/Transformation.hs view
@@ -1,5 +1,5 @@-{-# LANGUAGE FlexibleInstances #-} - +{-# LANGUAGE FlexibleInstances #-}+ -- | -- Module: Data.Geo.Jord.Transformation -- Copyright: (c) 2018 Cedric Liegeois @@ -15,178 +15,178 @@ -- -- 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 - +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 position and /n/-vector and height. -class NTransform a where - toNVector :: a -> AngularPosition NVector -- ^ position to 'AngularPosition' of 'NVector'. +class NTransform a where+ toNVector :: a -> AngularPosition NVector -- ^ data to 'AngularPosition' of 'NVector'. fromNVector :: AngularPosition NVector -> a -- ^ 'AngularPosition' of 'NVector' and height to position. - + -- | 'NVector' to, from 'AngularPosition' of 'NVector'. -instance NTransform NVector where - toNVector nv = AngularPosition nv zero - fromNVector = pos - +instance NTransform NVector where+ toNVector nv = AngularPosition nv zero+ fromNVector = pos+ -- | 'LatLong' to, from 'AngularPosition' of 'NVector'. -instance NTransform LatLong where - toNVector ll = AngularPosition (latLongToNVector ll) zero - fromNVector = nvectorToLatLong . pos - +instance NTransform LatLong where+ toNVector ll = AngularPosition (latLongToNVector ll) zero+ fromNVector = nvectorToLatLong . pos+ -- | 'NTransform' identity. -instance NTransform (AngularPosition NVector) where - toNVector = id - fromNVector = id +instance NTransform (AngularPosition NVector) where+ toNVector = id+ fromNVector = id -- | 'AngularPosition' of 'LatLong' to, from 'AngularPosition' of 'NVector'. -instance NTransform (AngularPosition LatLong) where - toNVector (AngularPosition ll h) = AngularPosition (latLongToNVector ll) h - fromNVector (AngularPosition nv h) = AngularPosition (nvectorToLatLong nv) h - +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 +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' to, from 'EcefPosition'. -instance ETransform NVector where - fromEcef p e = pos (ecefToNVector p e) - toEcef v = nvectorToEcef (nvectorHeight v zero) - +instance ETransform NVector where+ fromEcef p e = pos (ecefToNVector p e)+ toEcef v = nvectorToEcef (nvectorHeight v zero)+ -- | 'LatLong' to, from 'EcefPosition'. -instance ETransform LatLong where - fromEcef p e = fromNVector (nvectorHeight (fromEcef p e :: NVector) zero) - toEcef = toEcef . toNVector - +instance ETransform LatLong where+ fromEcef p e = fromNVector (nvectorHeight (fromEcef p e :: NVector) zero)+ toEcef = toEcef . toNVector+ -- | 'AngularPosition' of 'NVector' to, from 'EcefPosition'. -instance ETransform (AngularPosition NVector) where - fromEcef = ecefToNVector - toEcef = nvectorToEcef - +instance ETransform (AngularPosition NVector) where+ fromEcef = ecefToNVector+ toEcef = nvectorToEcef+ -- | 'AngularPosition' of 'LatLong' to, from 'EcefPosition'. -instance ETransform (AngularPosition LatLong) where - fromEcef p e = fromNVector (ecefToNVector p e) - toEcef = nvectorToEcef . toNVector - +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 - +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) - +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 - +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 +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) +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 - +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 +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') +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) - +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) +geodeticHeight :: EcefPosition -> Earth -> Length+geodeticHeight p e = height (ecefToNVector p e)
test/Data/Geo/Jord/GeodeticsSpec.hs view
@@ -1,213 +1,251 @@-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 +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 "alongTrackDistance" $ do+ it "returns a positive length when position is ahead start of great arc" $ do+ let p = decimalLatLong 53.2611 (-0.7972)+ let ga = greatArc (decimalLatLong 53.3206 (-1.7297), decimalLatLong 53.1887 0.1334)+ alongTrackDistance p ga r84 `shouldBe` kilometres 62.3315757+ it "returns a negative length when position is ahead start of great arc" $ do+ let p = decimalLatLong 53.3206 (-1.7297)+ let ga = greatArc (decimalLatLong 53.2611 (-0.7972), decimalLatLong 53.1887 0.1334)+ alongTrackDistance p ga r84 `shouldBe` kilometres (-62.3293209)+ it "returns a 0 when position is start of great arc" $ do+ let p = decimalLatLong 53.2611 (-0.7972)+ let ga = greatArc (p, decimalLatLong 53.1887 0.1334)+ alongTrackDistance p ga r84 `shouldBe` zero+ 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 "intersection" $ do+ it "returns nothing if both great arc are equals" $ do+ let ga = greatArc (decimalLatLong 51.885 0.235, decimalLatLong 52.885 1.235)+ (intersection ga ga :: Maybe LatLong) `shouldBe` Nothing+ it "returns nothing if both great arc are equals (opposite orientation)" $ do+ let ga1 = greatArc (decimalLatLong 51.885 0.235, decimalLatLong 52.885 1.235)+ let ga2 = greatArc (decimalLatLong 51.885 0.235, decimalLatLong 52.885 1.235)+ (intersection ga1 ga2 :: Maybe LatLong) `shouldBe` Nothing+ it "returns nothing if great circle intersection is outside either great arc" $ do+ let ga1 = greatArc (decimalLatLong 0 0, decimalLatLong 0 10)+ let ga2 = greatArc (decimalLatLong (-5) 5, decimalLatLong (-1) 5)+ (intersection ga1 ga2 :: Maybe LatLong) `shouldBe` Nothing+ it "returns nothing if great circle intersection is outside both great arcs" $ do+ let ga1 = greatArc (decimalLatLong 0 (-10), decimalLatLong 0 (-1))+ let ga2 = greatArc (decimalLatLong (-5) 5, decimalLatLong (-1) 5)+ (intersection ga1 ga2 :: Maybe LatLong) `shouldBe` Nothing+ it "returns the point where the two great arcs intersect" $ do+ let spd = kilometresPerHour 1000+ let t1 = Track (decimalLatLong 51.885 0.235) (decimalDegrees 108.63) spd+ let t2 = Track (decimalLatLong 49.008 2.549) (decimalDegrees 32.72) spd+ let oneHour = hours 1+ let ga1 = greatArc (t1, oneHour)+ let ga2 = greatArc (t2, oneHour)+ (intersection ga1 ga2 :: Maybe LatLong) `shouldBe`+ Just (decimalLatLong 50.9017225 4.494278333333333)+ 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 intersect" $ 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 "isInsideSurface" $ 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" $ isInsideSurface p1 [] `shouldBe` False+ it "return False if polygon does not define at least a triangle" $+ isInsideSurface p1 [p1, p2] `shouldBe` False+ it "returns True if point is inside polygon" $ do+ let polygon = [p1, p2, p4, p3]+ isInsideSurface p5 polygon `shouldBe` True+ it "returns False if point is inside polygon" $ do+ let polygon = [p1, p2, p4, p3]+ let p = antipode p5+ isInsideSurface p polygon `shouldBe` False+ it "returns False if point is a vertex of the polygon" $ do+ let polygon = [p1, p2, p4, p3]+ isInsideSurface p1 polygon `shouldBe` False+ it "handles closed polygons" $ do+ let polygon = [p1, p2, p4, p3, p1]+ isInsideSurface 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+ isInsideSurface hoor polygon `shouldBe` True+ isInsideSurface hassleholm polygon `shouldBe` False+ 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
@@ -1,199 +1,199 @@-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 behing target" $ do - let t = Track (decimalLatLong 45 67) (decimalDegrees 54) (knots 400) - let ip = decimalLatLong 44 66 - 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 Nothing if interceptor is behing target" $ do - let t = Track (decimalLatLong 45 67) (decimalDegrees 181) (knots 400) - let ip = decimalLatLong 44 66 - let i = intercept84 t ip - fmap interceptorSpeed i `shouldBe` Just (knots 228.5538171521) - fmap interceptTime i `shouldBe` Just (seconds 808.770) - let interceptor = - Track - ip - (fromJust (fmap interceptorBearing i)) - (fromJust (fmap interceptorSpeed i)) - fmap interceptPosition i `shouldBe` - Just (position84 interceptor (fromJust (fmap interceptTime i))) - 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.633367756059) - fmap interceptTime i `shouldBe` Just (seconds 5993.831) - 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 "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 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) - it "returns the same as intercept when called with minimum speed" $ do - let t = Track (decimalLatLong 45 50) (decimalDegrees 54) (knots 500) - let ip = decimalLatLong 70 30 - let mi = intercept84 t ip - let i = interceptBySpeed84 t ip (fromJust (fmap interceptorSpeed mi)) - fmap interceptTime i `shouldBe` fmap interceptTime mi - 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" $ - -- 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. - do - 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 +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 behing target" $ do+ let t = Track (decimalLatLong 45 67) (decimalDegrees 54) (knots 400)+ let ip = decimalLatLong 44 66+ 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 Nothing if interceptor is behing target" $ do+ let t = Track (decimalLatLong 45 67) (decimalDegrees 181) (knots 400)+ let ip = decimalLatLong 44 66+ let i = intercept84 t ip+ fmap interceptorSpeed i `shouldBe` Just (knots 228.5538171521)+ fmap interceptTime i `shouldBe` Just (seconds 808.770)+ let interceptor =+ Track+ ip+ (fromJust (fmap interceptorBearing i))+ (fromJust (fmap interceptorSpeed i))+ fmap interceptPosition i `shouldBe`+ Just (position84 interceptor (fromJust (fmap interceptTime i)))+ 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.633367756059)+ fmap interceptTime i `shouldBe` Just (seconds 5993.831)+ 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 "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 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)+ it "returns the same as intercept when called with minimum speed" $ do+ let t = Track (decimalLatLong 45 50) (decimalDegrees 54) (knots 500)+ let ip = decimalLatLong 70 30+ let mi = intercept84 t ip+ let i = interceptBySpeed84 t ip (fromJust (fmap interceptorSpeed mi))+ fmap interceptTime i `shouldBe` fmap interceptTime mi+ 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" $+ -- 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.+ do+ 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/SpeedSpec.hs view
@@ -30,8 +30,8 @@ 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 m/s" $ toMetresPerSecond (metresPerSecond 1) `shouldBe` 1 + it "handles 1 mph" $ 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