diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/app/Eval.hs b/app/Eval.hs --- a/app/Eval.hs +++ b/app/Eval.hs @@ -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) diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -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 diff --git a/app/Show.hs b/app/Show.hs --- a/app/Show.hs +++ b/app/Show.hs @@ -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)) ++ ")" diff --git a/app/State.hs b/app/State.hs --- a/app/State.hs +++ b/app/State.hs @@ -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 diff --git a/benchmarks/GeodeticsBG.hs b/benchmarks/GeodeticsBG.hs new file mode 100644 --- /dev/null +++ b/benchmarks/GeodeticsBG.hs @@ -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 diff --git a/benchmarks/KinematicsBG.hs b/benchmarks/KinematicsBG.hs new file mode 100644 --- /dev/null +++ b/benchmarks/KinematicsBG.hs @@ -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)) diff --git a/benchmarks/Main.hs b/benchmarks/Main.hs new file mode 100644 --- /dev/null +++ b/benchmarks/Main.hs @@ -0,0 +1,9 @@ +module Main where + +import Criterion.Main +import GeodeticsBG +import KinematicsBG +import TransformationBG + +main :: IO () +main = defaultMain [bggeodetics, bgkinematics, bgtransformation] diff --git a/benchmarks/TransformationBG.hs b/benchmarks/TransformationBG.hs new file mode 100644 --- /dev/null +++ b/benchmarks/TransformationBG.hs @@ -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 diff --git a/jord.cabal b/jord.cabal --- a/jord.cabal +++ b/jord.cabal @@ -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 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 diff --git a/src/Data/Geo/Jord.hs b/src/Data/Geo/Jord.hs --- a/src/Data/Geo/Jord.hs +++ b/src/Data/Geo/Jord.hs @@ -55,4 +55,4 @@ -- | version. jordVersion :: String -jordVersion = "0.4.2.0" +jordVersion = "0.5.0.0" diff --git a/src/Data/Geo/Jord/Angle.hs b/src/Data/Geo/Jord/Angle.hs --- a/src/Data/Geo/Jord/Angle.hs +++ b/src/Data/Geo/Jord/Angle.hs @@ -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 diff --git a/src/Data/Geo/Jord/Duration.hs b/src/Data/Geo/Jord/Duration.hs --- a/src/Data/Geo/Jord/Duration.hs +++ b/src/Data/Geo/Jord/Duration.hs @@ -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 diff --git a/src/Data/Geo/Jord/Geodetics.hs b/src/Data/Geo/Jord/Geodetics.hs --- a/src/Data/Geo/Jord/Geodetics.hs +++ b/src/Data/Geo/Jord/Geodetics.hs @@ -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 diff --git a/src/Data/Geo/Jord/Internal.hs b/src/Data/Geo/Jord/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geo/Jord/Internal.hs @@ -0,0 +1,37 @@ +-- | +-- Module: Data.Geo.Jord.Internal +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois +-- 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 diff --git a/src/Data/Geo/Jord/Kinematics.hs b/src/Data/Geo/Jord/Kinematics.hs --- a/src/Data/Geo/Jord/Kinematics.hs +++ b/src/Data/Geo/Jord/Kinematics.hs @@ -1,37 +1,39 @@ --- | --- Module: Data.Geo.Jord.Kinematics --- Copyright: (c) 2018 Cedric Liegeois --- License: BSD3 --- Maintainer: Cedric Liegeois --- 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 --- --- and in --- +{-# LANGUAGE FlexibleInstances #-} + +-- | +-- Module: Data.Geo.Jord.Kinematics +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois +-- 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 +-- +-- and in +-- 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 diff --git a/src/Data/Geo/Jord/LatLong.hs b/src/Data/Geo/Jord/LatLong.hs --- a/src/Data/Geo/Jord/LatLong.hs +++ b/src/Data/Geo/Jord/LatLong.hs @@ -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 diff --git a/src/Data/Geo/Jord/Length.hs b/src/Data/Geo/Jord/Length.hs --- a/src/Data/Geo/Jord/Length.hs +++ b/src/Data/Geo/Jord/Length.hs @@ -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 diff --git a/src/Data/Geo/Jord/Parse.hs b/src/Data/Geo/Jord/Parse.hs deleted file mode 100644 --- a/src/Data/Geo/Jord/Parse.hs +++ /dev/null @@ -1,52 +0,0 @@ --- | --- Module: Data.Geo.Jord.Parse --- Copyright: (c) 2018 Cedric Liegeois --- License: BSD3 --- Maintainer: Cedric Liegeois --- 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 diff --git a/src/Data/Geo/Jord/Parser.hs b/src/Data/Geo/Jord/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geo/Jord/Parser.hs @@ -0,0 +1,52 @@ +-- | +-- Module: Data.Geo.Jord.Parser +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois +-- 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 diff --git a/src/Data/Geo/Jord/Speed.hs b/src/Data/Geo/Jord/Speed.hs --- a/src/Data/Geo/Jord/Speed.hs +++ b/src/Data/Geo/Jord/Speed.hs @@ -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 diff --git a/src/Data/Geo/Jord/Transformation.hs b/src/Data/Geo/Jord/Transformation.hs --- a/src/Data/Geo/Jord/Transformation.hs +++ b/src/Data/Geo/Jord/Transformation.hs @@ -1,5 +1,5 @@ -{-# LANGUAGE FlexibleInstances #-} - +{-# LANGUAGE FlexibleInstances #-} + -- | -- Module: Data.Geo.Jord.Transformation -- Copyright: (c) 2018 Cedric Liegeois @@ -15,178 +15,178 @@ -- -- See -- -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) diff --git a/test/Data/Geo/Jord/GeodeticsSpec.hs b/test/Data/Geo/Jord/GeodeticsSpec.hs --- a/test/Data/Geo/Jord/GeodeticsSpec.hs +++ b/test/Data/Geo/Jord/GeodeticsSpec.hs @@ -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 diff --git a/test/Data/Geo/Jord/KinematicsSpec.hs b/test/Data/Geo/Jord/KinematicsSpec.hs --- a/test/Data/Geo/Jord/KinematicsSpec.hs +++ b/test/Data/Geo/Jord/KinematicsSpec.hs @@ -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 diff --git a/test/Data/Geo/Jord/SpeedSpec.hs b/test/Data/Geo/Jord/SpeedSpec.hs --- a/test/Data/Geo/Jord/SpeedSpec.hs +++ b/test/Data/Geo/Jord/SpeedSpec.hs @@ -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