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