diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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. diff --git a/geodetics.cabal b/geodetics.cabal --- a/geodetics.cabal +++ b/geodetics.cabal @@ -1,5 +1,5 @@ name: geodetics -version: 0.1.0 +version: 0.1.2 cabal-version: >= 1.10 build-type: Simple author: Paul Johnson diff --git a/src/Geodetics/Geodetic.hs b/src/Geodetics/Geodetic.hs --- a/src/Geodetics/Geodetic.hs +++ b/src/Geodetics/Geodetic.hs @@ -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 - diff --git a/src/Geodetics/LatLongParser.hs b/src/Geodetics/LatLongParser.hs --- a/src/Geodetics/LatLongParser.hs +++ b/src/Geodetics/LatLongParser.hs @@ -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) - diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 -- 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 -- , 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 -- 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