geodetics 0.1.0 → 0.1.2
raw patch · 5 files changed
+129/−105 lines, 5 filesPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
- Geodetics.Geodetic: data (Ellipsoid e) => Geodetic e
+ Geodetics.Geodetic: data Geodetic e
Files
- changelog.md +5/−0
- geodetics.cabal +1/−1
- src/Geodetics/Geodetic.hs +35/−36
- src/Geodetics/LatLongParser.hs +42/−31
- test/Main.hs +46/−37
changelog.md view
@@ -15,3 +15,8 @@ and fix the build on 7.10 with a conditional semigroups dependency Version 0.1.0: Updated for Dimensional 1.3 and GHC 8.6.++Version 0.1.1: Fixed bug #15: for a point p, groundDistance p p returned NaN++Version 0.1.2: Fixed bugs #16 and #17: Unicode PRIME and DOUBLE PRIME now allowed in+ position strings, and the degree symbol is allowed for decimal degrees.
geodetics.cabal view
@@ -1,5 +1,5 @@ name: geodetics-version: 0.1.0+version: 0.1.2 cabal-version: >= 1.10 build-type: Simple author: Paul Johnson <paul@cogito.org.uk>
src/Geodetics/Geodetic.hs view
@@ -33,24 +33,24 @@ -- | Defines a three-D position on or around the Earth using latitude, -- longitude and altitude with respect to a specified ellipsoid, with -- positive directions being North and East. The default "show"--- instance gives position in degrees, minutes and seconds to 5 decimal +-- instance gives position in degrees, minutes and seconds to 5 decimal -- places, which is a -- resolution of about 1m on the Earth's surface. Internally latitude -- and longitude are stored as double precision radians. Convert to -- degrees using e.g. @latitude g /~ degree@.--- +-- -- The functions here deal with altitude by assuming that the local -- height datum is always co-incident with the ellipsoid in use, -- even though the \"mean sea level\" (the usual height datum) can be tens -- of meters above or below the ellipsoid, and two ellipsoids can -- differ by similar amounts. This is because the altitude is -- usually known with reference to a local datum regardless of the--- ellipsoid in use, so it is simpler to preserve the altitude across +-- ellipsoid in use, so it is simpler to preserve the altitude across -- all operations. However if -- you are working with ECEF coordinates from some other source then -- this may give you the wrong results, depending on the altitude -- correction your source has used.--- +-- -- There is no "Eq" instance because comparing two arbitrary -- co-ordinates on the Earth is a non-trivial exercise. Clearly if all -- the parameters are equal on the same ellipsoid then they are indeed@@ -59,7 +59,7 @@ -- physical location. If you want to find out if two co-ordinates are -- the same to within a given tolerance then use "geometricDistance" -- (or its squared variant to avoid an extra @sqrt@ operation).-data (Ellipsoid e) => Geodetic e = Geodetic {+data Geodetic e = Geodetic { latitude, longitude :: Angle Double, geoAlt :: Length Double, ellipsoid :: e@@ -68,43 +68,43 @@ instance (Ellipsoid e) => Show (Geodetic e) where show g = concat [ showAngle (abs $ latitude g), " ", letter "SN" (latitude g), ", ",- showAngle (abs $ longitude g), " ", letter "WE" (longitude g), ", ", + showAngle (abs $ longitude g), " ", letter "WE" (longitude g), ", ", show (altitude g), " ", show (ellipsoid g)] where letter s n = [s !! (if n < _0 then 0 else 1)] --- | Read the latitude and longitude of a ground position and +-- | Read the latitude and longitude of a ground position and -- return a Geodetic position on the specified ellipsoid.--- +-- -- The latitude and longitude may be in any of the following formats.--- The comma between latitude and longitude is optional in all cases. +-- The comma between latitude and longitude is optional in all cases. -- Latitude must always be first.--- +-- -- * Signed decimal degrees: 34.52327, -46.23234--- +-- -- * Decimal degrees NSEW: 34.52327N, 46.23234W ----- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92' --- --- * Degrees, minutes and seconds (units optional): 34° 31' 23.52\" N, 46° 13' 56.43\" W --- +-- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92'+--+-- * Degrees, minutes and seconds (units optional): 34° 31' 23.52\" N, 46° 13' 56.43\" W+-- -- * DDDMMSS format with optional leading zeros: 343123.52N, 0461356.43W readGroundPosition :: (Ellipsoid e) => e -> String -> Maybe (Geodetic e)-readGroundPosition e str = +readGroundPosition e str = case map fst $ filter (null . snd) $ readP_to_S latLong str of [] -> Nothing (lat,long) : _ -> Just $ groundPosition $ Geodetic (lat *~ degree) (long *~ degree) undefined e- - ++ -- | Show an angle as degrees, minutes and seconds to two decimal places. showAngle :: Angle Double -> String showAngle a | isNaN a1 = "NaN" -- Not a Nangle | isInfinite a1 = sgn ++ "Infinity"- | otherwise = concat [sgn, show d, [chr 0xB0, ' '], - show m, "' ", - show s, ".", dstr, "\"" ]+ | otherwise = concat [sgn, show d, [chr 0xB0, ' '],+ show m, "\8242 ",+ show s, ".", dstr, "\8243" ] where a1 = a /~ one sgn = if a < _0 then "-" else ""@@ -114,14 +114,14 @@ (m, s1) = m1 `P.divMod` 6000 -- hundredths of arcsec per arcmin (s, ds) = s1 `P.divMod` 100 dstr = reverse $ take 2 $ reverse (show ds) ++ "00" -- Decimal fraction with zero padding.- + instance (Ellipsoid e) => HasAltitude (Geodetic e) where altitude = geoAlt setAltitude h g = g{geoAlt = h} - - ++ -- | The point on the Earth diametrically opposite the argument, with -- the same altitude. antipode :: (Ellipsoid e) => Geodetic e -> Geodetic e@@ -130,10 +130,10 @@ lat = negate $ latitude g long' = longitude g - 180 *~ degree long | long' < _0 = long' + 360 *~ degree- | otherwise = long' + | otherwise = long' - - ++ -- | Convert a geodetic coordinate into earth centered, relative to the -- ellipsoid in use. geoToEarth :: (Ellipsoid e) => Geodetic e -> ECEF@@ -141,7 +141,7 @@ (n + h) * coslat * coslong, (n + h) * coslat * sinlong, (n * (_1 - eccentricity2 e) + h) * sinlat)- where + where n = normal e $ latitude geo e = ellipsoid geo coslat = cos $ latitude geo@@ -151,7 +151,7 @@ h = altitude geo --- | Convert an earth centred coordinate into a geodetic coordinate on +-- | Convert an earth centred coordinate into a geodetic coordinate on -- the specified geoid. -- -- Uses the closed form solution of H. Vermeille: Direct@@ -199,9 +199,9 @@ h = helmert (ellipsoid g) --- | The absolute distance in a straight line between two geodetic +-- | The absolute distance in a straight line between two geodetic -- points. They must be on the same ellipsoid.--- Note that this is not the geodetic distance taken by following +-- Note that this is not the geodetic distance taken by following -- the curvature of the earth. geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Length Double geometricalDistance g1 g2 = sqrt $ geometricalDistanceSq g1 g2@@ -261,7 +261,7 @@ cosU1 = cos u1 sinU2 = sin u2 cosU2 = cos u2- + nextLambda lambda = (lambda1, (cos2Alpha, delta, sinDelta, cosDelta, cos2DeltaM)) where sinLambda = sin lambda@@ -270,7 +270,7 @@ (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ pos2) cosDelta = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda delta = atan2 sinDelta cosDelta- sinAlpha = cosU1 * cosU2 * sinLambda / sinDelta+ sinAlpha = if sinDelta == _0 then _0 else cosU1 * cosU2 * sinLambda / sinDelta cos2Alpha = _1 - sinAlpha ^ pos2 cos2DeltaM = if cos2Alpha == _0 then _0@@ -285,14 +285,13 @@ -- | Add or subtract multiples of 2*pi so that for all @t@, @-pi < properAngle t < pi@. properAngle :: Angle Double -> Angle Double-properAngle t +properAngle t | r1 <= negate pi = r1 + pi2 | r1 > pi = r1 - pi2- | otherwise = r1 + | otherwise = r1 where pf :: Double -> (Int, Double) pf = properFraction -- Shut up GHC warning about defaulting to Integer. (_,r) = pf (t/pi2 /~ one) r1 = (r *~ one) * pi2 pi2 = pi * _2-
src/Geodetics/LatLongParser.hs view
@@ -21,11 +21,23 @@ import Text.ParserCombinators.ReadP as P - -- | Parse an unsigned Integer value. natural :: ReadP Integer -- Beware arithmetic overflow of Int natural = read <$> munch1 isDigit ++-- | Parse a tick sign for minutes. This accepts either the keyboard \"'\" or the unicode \"Prime\"+-- character U+2032+minuteTick :: ReadP ()+minuteTick = void $ choice [char '\'', char '\8242']+++-- | Parse a double-tick sign for seconds. This accepts either the keyboard \" or the unicode+-- \"Double Prime\" character U+2033.+secondTick :: ReadP ()+secondTick = void $ choice [char '"', char '\8243']++ -- | Parse an unsigned decimal value with optional decimal places but no exponent. decimal :: ReadP Double decimal = do@@ -33,8 +45,8 @@ option (read str1) $ do str2 <- char '.' *> munch1 isDigit return $ read $ str1 ++ '.' : str2- - ++ -- | Read a character indicating the sign of a value. Returns either +1 or -1. signChar :: (Num a) => Char -- ^ Positive sign@@ -43,12 +55,12 @@ signChar pos neg = do c <- char pos +++ char neg return $ if c == pos then 1 else (-1)- - ++ -- | Parse a signed decimal value. signedDecimal :: ReadP Double-signedDecimal = (*) <$> option 1 (signChar '+' '-') <*> decimal - +signedDecimal = (*) <$> option 1 (signChar '+' '-') <*> decimal+ -- | Parse an unsigned angle written using degrees, minutes and seconds separated by spaces. -- All except the last must be integers. degreesMinutesSeconds :: ReadP Double@@ -73,41 +85,41 @@ d <- fromIntegral <$> option 0 (natural <* char '°') guard $ d <= 360 skipSpaces- m <- fromIntegral <$> option 0 (natural <* char '\'')+ m <- fromIntegral <$> option 0 (natural <* minuteTick) guard $ m < 60 skipSpaces- s <- option 0 (decimal <* char '"')+ s <- option 0 (decimal <* secondTick) guard $ s < 60 return $ d + m / 60 + s / 3600 guard $ not $ null s -- Must specify at least one component. return a- --- | Parse an unsigned angle written using degrees and decimal minutes. ++-- | Parse an unsigned angle written using degrees and decimal minutes. degreesDecimalMinutes :: ReadP Double degreesDecimalMinutes = do d <- fromIntegral <$> natural skipSpaces guard $ d <= 360 -- Difference from degreesMinutesSeconds just to shut style checker up. m <- option 0 decimal- guard $ m < 60 + guard $ m < 60 return $ d + m/60- - ++ -- | Parse an unsigned angle written using degrees and decimal minutes with units (° ') degreesDecimalMinutesUnits :: ReadP Double degreesDecimalMinutesUnits = do (s, a) <- gather $ do- d <- fromIntegral <$> option 0 (natural <* char '°')+ d <- fromIntegral <$> option 0 (natural <* char '°') guard $ d <= 360- m <- option 0 (decimal <* char '\'')+ m <- option 0 (decimal <* minuteTick) guard $ m < 60 return $ d + m / 60 guard $ not $ null s -- Must specify at least one component. return a --- | Parse an unsigned angle written in DDDMMSS.ss format. +-- | Parse an unsigned angle written in DDDMMSS.ss format. -- Leading zeros on the degrees and decimal places on the seconds are optional dms7 :: ReadP Double dms7 = do@@ -126,18 +138,18 @@ -- | Parse an unsigned angle, either in decimal degrees or in degrees, minutes and seconds.--- In the latter case the unit indicators are optional. +-- In the latter case the unit indicators are optional. angle :: ReadP Double angle = choice [- decimal, + decimal <* optional (char '°'), degreesMinutesSeconds, degreesMinutesSecondsUnits, degreesDecimalMinutes, degreesDecimalMinutesUnits, dms7 ]- + -- | Parse latitude as an unsigned angle followed by 'N' or 'S' latitudeNS :: ReadP Double latitudeNS = do@@ -147,33 +159,33 @@ sgn <- signChar 'N' 'S' return $ sgn * ul - + -- | Parse longitude as an unsigned angle followed by 'E' or 'W'. longitudeEW :: ReadP Double longitudeEW = do- ul <- angle + ul <- angle guard $ ul <= 180 skipSpaces sgn <- signChar 'E' 'W' return $ sgn * ul- - --- | Parse latitude and longitude as two signed decimal numbers in that order, optionally separated by a comma. +++-- | Parse latitude and longitude as two signed decimal numbers in that order, optionally separated by a comma. -- Longitudes in the western hemisphere may be represented either by negative angles down to -180 -- or by positive angles less than 360. signedLatLong :: ReadP (Double, Double) signedLatLong = do- lat <- signedDecimal+ lat <- signedDecimal <* optional (char '°') guard $ lat >= (-90) guard $ lat <= 90 skipSpaces P.optional $ char ',' >> skipSpaces- long <- signedDecimal+ long <- signedDecimal <* optional (char '°') guard $ long >= (-180) guard $ long < 360- return (lat, if long > 180 then long-180 else long)- - + return (lat, if long > 180 then long-360 else long)++ -- | Parse latitude and longitude in any format. latLong :: ReadP (Double, Double) latLong = latLong1 +++ longLat +++ signedLatLong@@ -190,4 +202,3 @@ P.optional $ char ',' >> skipSpaces lat <- latitudeNS return (lat, long)-
test/Main.hs view
@@ -59,6 +59,7 @@ tests = [ testGroup "Geodetic" [ testProperty "WGS84 and back" prop_WGS84_and_back,+ testProperty "Zero ground distance" prop_zero_ground, testGroup "UK Points" $ map pointTest ukPoints], testGroup "World lines" $ map worldLineTests worldLines, testGroup "Grid" [@@ -123,7 +124,7 @@ closeGrid p1 p2 = check eastings && check northings && check altitude where check f = f p1 - f p2 < 1 *~ meter --- | Degrees, minutes and seconds into radians. +-- | Degrees, minutes and seconds into radians. dms :: Int -> Int -> Double -> Dimensionless Double dms d m s = fromIntegral d *~ degree + fromIntegral m *~ arcminute + s *~ arcsecond @@ -132,7 +133,15 @@ prop_WGS84_and_back p = samePlace p $ toLocal (ellipsoid p) $ toWGS84 p --- | Sample pairs of points with bearings and distances. +-- | Test that for all points p, the ground distance from p to p is zero.+prop_zero_ground :: Geodetic WGS84 -> Bool+prop_zero_ground p =+ case groundDistance p p of+ Nothing -> False+ Just (d, _, _) -> abs d < 1 *~ milli meter+++-- | Sample pairs of points with bearings and distances. -- The Oracle for these values is the @FORWARD@ program from -- <http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html> worldLines :: [(String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double)]@@ -143,25 +152,25 @@ 6695785.820*~meter, 0*~degree, 180*~degree), ("Equator to Pole", Geodetic (0*~degree) (0*~degree) _0 WGS84, Geodetic (90*~degree) (180*~degree) _0 WGS84, 10001965.729*~meter, 0*~degree, 180*~degree)]- - ++ worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double) -> Test worldLineTests (str, g1, g2, d, a, b) = testCase str $ HU.assertBool "" $ ok $ groundDistance g1 g2 where ok Nothing = False- ok (Just (d1, a1, b1)) = - abs (d - d1) < 0.01 *~ meter - && abs (a - a1) < 0.01 *~ arcsecond + ok (Just (d1, a1, b1)) =+ abs (d - d1) < 0.01 *~ meter+ && abs (a - a1) < 0.01 *~ arcsecond && abs (b - b1) < 0.01 *~ arcsecond --- | Sample points for UK tests. The oracle for these values is the script at +-- | Sample points for UK tests. The oracle for these values is the script at -- <http://www.movable-type.co.uk/scripts/latlong-convert-coords.html>, which uses -- the same Helmert transform as this library. Hence the results should match to within 30 cm. ukPoints :: [(String, Geodetic WGS84, Geodetic OSGB36)] ukPoints = [- ("Greenwich", Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) _0 WGS84, + ("Greenwich", Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) _0 WGS84, Geodetic (dms 51 28 39.00) (dms 0 0 0) _0 OSGB36),- ("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) _0 WGS84, + ("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) _0 WGS84, Geodetic (dms 55 56 56.51) (dms (-3) (-11) (-57.61)) _0 OSGB36), ("Lands End", Geodetic (dms 50 03 56.68) (dms (-5) (-42) (-51.20)) _0 WGS84, Geodetic (dms 50 03 54.51) (dms (-5) (-42) (-47.87)) _0 OSGB36),@@ -184,13 +193,13 @@ -- A polar offset multiplied by a scalar is equal to an offset in the same direction with the length multiplied. prop_offset2 :: Distance -> Bearing -> Scalar -> Bool prop_offset2 (Distance d) (Bearing h) (Scalar s) = sameOffset go1 go2- where + where go1 = offsetScale s $ polarOffset d h go2 = polarOffset (d * s) h -- | A polar offset has the offset distance and bearing of its arguments. prop_offset3 :: GridOffset -> Bool-prop_offset3 delta = sameOffset delta0 +prop_offset3 delta = sameOffset delta0 (polarOffset (offsetDistance delta0) (offsetBearing delta)) where delta0 = delta {deltaAltitude = 0 *~ meter} @@ -202,11 +211,11 @@ -- | Converting a UK grid reference to a GridPoint and back is a null operation. prop_ukGrid1 :: GridRef -> Bool-prop_ukGrid1 (GridRef str) = +prop_ukGrid1 (GridRef str) = str == (fromJust $ toUkGridReference ((length str P.- 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str) --- | UK Grid Reference points. The oracle for these points was the +-- | UK Grid Reference points. The oracle for these points was the -- UK Grid Reference Finder (gridreferencefinder.com), retrieved on 26 Jan 2013. ukSampleGrid :: [(String, GridPoint UkNationalGrid, Geodetic WGS84, String)] ukSampleGrid = map convert [@@ -221,7 +230,7 @@ ("ST1922474591", 319224, 174591, 51.464505, -3.1641741, "Torchwood HQ"), ("SK3520736502", 435207, 336502, 52.924784, -1.4777486, "Derby Cathedral")] where- convert (grid, x, y, lat, long, desc) = + convert (grid, x, y, lat, long, desc) = (grid, GridPoint (x *~ meter) (y *~ meter) (0 *~ meter) UkNationalGrid, Geodetic (lat *~ degree) (long *~ degree) (0 *~ meter) WGS84, desc) @@ -229,19 +238,19 @@ -- | Check that grid reference to grid point works for sample points. ukGridTest2 :: GridPointTest-ukGridTest2 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool "" +ukGridTest2 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool "" $ (fst $ fromJust $ fromUkGridReference gridRef) == gp -- | Check that grid point to grid reference works for sample points. ukGridTest3 :: GridPointTest-ukGridTest3 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool "" +ukGridTest3 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool "" $ toUkGridReference 5 gp == Just gridRef --- | Check that grid point to WGS84 works close enough for sample points. +-- | Check that grid point to WGS84 works close enough for sample points. ukGridTest4 :: GridPointTest ukGridTest4 (_, gp, geo, testName) = testCase testName $ HU.assertBool "" $ closeEnough geo $ toWGS84 $ fromGrid gp- + -- | Check that WGS84 to grid point works close enough for sample points. ukGridTest5 :: GridPointTest ukGridTest5 (_, gp, geo, testName) = testCase testName $ HU.assertBool ""@@ -252,7 +261,7 @@ ukTest :: Geodetic OSGB36 ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) (0 *~ meter) OSGB36 -{- +{- v = 6.3885023333E+06 rho = 6.3727564399E+06 eta2 = 2.4708136169E-03@@ -276,11 +285,11 @@ ellipse = LocalEllipsoid "Bessel 1841" (6377397.155 *~ metre) (299.15281 *~ one) mempty tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) (0 *~ meter) ellipse origin = GridOffset (155000 *~ metre) (463000 *~ metre) (0 *~ meter)- - ++ -- | Standard steregraphic grid for point tests in the Southern Hemisphere.--- --- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated. +--+-- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated. stereoGridS :: GridStereo LocalEllipsoid stereoGridS = mkGridStereo tangent origin (0.9999079 *~ one) where@@ -289,12 +298,12 @@ origin = GridOffset ((-155000) *~ metre) (463000 *~ metre) (0 *~ meter) --- | Data for the stereographic tests taken from +-- | Data for the stereographic tests taken from -- <http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf> stereographicToGridN :: Bool stereographicToGridN = sameGrid g1 g1' where- p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN + p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN g1' = toGrid stereoGridN p1 @@ -303,7 +312,7 @@ where p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN- p1' = fromGrid g1 + p1' = fromGrid g1 stereographicToGridS :: Bool@@ -319,7 +328,7 @@ where p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridS g1 = GridPoint ((-196105.283) *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridS- p1' = fromGrid g1 + p1' = fromGrid g1 -- | Check the round trip for a stereographic projection.@@ -328,13 +337,13 @@ let g = fromGrid p r = toGrid (gridBasis p) g in counterexample ("p = " ++ show p ++ "\ng = " ++ show g ++ "\nr = " ++ show r) $- closeGrid p r + closeGrid p r -- | A ray at distance zero returns its original arguments. prop_rayPath1 :: Ray WGS84 -> Bool-prop_rayPath1 r@(Ray pt b e) = +prop_rayPath1 r@(Ray pt b e) = samePlace pt pt1 && sameAngle b b1 && sameAngle e e1 where (pt1,b1,e1) = pathFunc (getRay r) _0 @@ -344,7 +353,7 @@ type ContinuityTest1 e = Geodetic e -> Bearing -> Distance2 -> Distance2 -> Property -- | Many paths can be specified by a start point, bearing and azimuth,--- and have the property that any (point,bearing,azimuth) triple on +-- and have the property that any (point,bearing,azimuth) triple on -- the path will specify the same path with a distance offset. prop_pathContinuity :: (Ellipsoid e) => (Geodetic e -> Angle Double -> Angle Double -> Path e) -> ContinuityTest e@@ -358,9 +367,9 @@ path1 = pf pt1 b1 a1 (pt2, b2, a2) = pathFunc path1 d2 (pt3, b3, a3) = pathFunc path0 (d1 + d2) -- Points 2 and 3 should be the same.- --- | For continuity testing of ground-based paths (azimuth & altitude always zero) ++-- | For continuity testing of ground-based paths (azimuth & altitude always zero) -- where lower accuracy is required. prop_pathContinuity1 :: (Ellipsoid e) => (Geodetic e -> Angle Double -> Path e) -> ContinuityTest1 e prop_pathContinuity1 pf pt0 (Bearing b0) (Distance2 d1) (Distance2 d2) =@@ -377,21 +386,21 @@ -- | A point on a ray will continue along the same ray, and hence give the same points. prop_rayContinuity :: ContinuityTest WGS84-prop_rayContinuity = prop_pathContinuity rayPath +prop_rayContinuity = prop_pathContinuity rayPath -- | A ray bisected to an altitude will give that altitude. -- This is a test of bisection rather than rays. prop_rayBisect :: Ray WGS84 -> Altitude -> Bool-prop_rayBisect r (Altitude height) = +prop_rayBisect r (Altitude height) = case bisect ray0 f (1 *~ centi meter) (0 *~ meter) (1000 *~ kilo meter) of Nothing -> False Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1 *~ centi meter where f g = compare (altitude g) height ray0 = getRay r- + -- | A point on a rhumb line will continue along the same rhumb. prop_rhumbContinuity :: ContinuityTest1 WGS84 prop_rhumbContinuity = prop_pathContinuity1 rhumbPath@@ -399,7 +408,7 @@ -- | Two rhumb paths intersect at the same place. prop_rhumbIntersect :: RhumbPaths2 -> Property-prop_rhumbIntersect rp = +prop_rhumbIntersect rp = case intersect _0 _0 (10.0 *~ centi meter) 100 path1 path2 of Just (d1, d2) -> let (pt1, _, _) = pathFunc path1 d1