-- | -- Module: Data.Geo.Jord.Geodesic -- Copyright: (c) 2020 Cedric Liegeois -- License: BSD3 -- Maintainer: Cedric Liegeois -- Stability: experimental -- Portability: portable -- -- Solutions to the direct and inverse geodesic problems on ellipsoidal models using Vincenty formulaes. -- A geodesic is the shortest path between two points on a curved surface - here an ellispoid. Using these -- functions improves on the accuracy available using "Data.Geo.Jord.GreatCircle" at the expense of higher -- CPU usage. -- -- In order to use this module you should start with the following imports: -- -- @ -- import qualified Data.Geo.Jord.Angle as Angle -- import Data.Geo.Jord.Geodesic (Geodesic) -- import qualified Data.Geo.Jord.Geodesic as Geodesic -- import qualified Data.Geo.Jord.Geodetic as Geodetic -- import qualified Data.Geo.Jord.Length as Length -- @ -- -- module Data.Geo.Jord.Geodesic ( -- * The 'Geodesic' type Geodesic , startPosition , endPosition , initialBearing , finalBearing , length -- * Calculations , direct , inverse ) where import Prelude hiding (length) import Data.Geo.Jord.Angle (Angle) import qualified Data.Geo.Jord.Angle as Angle import Data.Geo.Jord.Ellipsoid import Data.Geo.Jord.Geodetic (HorizontalPosition) import qualified Data.Geo.Jord.Geodetic as Geodetic import Data.Geo.Jord.Length (Length) import qualified Data.Geo.Jord.Length as Length import Data.Geo.Jord.Model (Ellipsoidal, surface) -- | Geodesic line: shortest route between two positions on the surface of a model. -- -- Bearing are in compass angles. -- Compass angles are clockwise angles from true north: 0° = north, 90° = east, 180° = south, 270° = west. -- The final bearing will differ from the initial bearing by varying degrees according to distance and latitude. data Geodesic a = Geodesic { startPosition :: HorizontalPosition a -- ^ geodesic start position. , endPosition :: HorizontalPosition a -- ^ geodesic end position. , initialBearing :: Maybe Angle -- ^ initial bearing from @startPosition@ to @endPosition@, if both are different. , finalBearing :: Maybe Angle -- ^ final bearing from @startPosition@ to @endPosition@, if both are different. , length :: Length -- ^ length of the geodesic: the surface distance between @startPosition@ to @endPosition@. } deriving (Eq, Show) -- | @direct p1 b1 d@ solves the direct geodesic problem using Vicenty formula: position -- along the geodesic, reached from position @p1@ having travelled the __surface__ distance @d@ on -- the initial bearing (compass angle) @b1@ at __constant__ height; it also returns the final bearing -- at the reached position. For example: -- -- >>> Geodesic.direct (Geodetic.northPole WGS84) Angle.zero (Length.kilometres 20003.931458623) -- Just (Geodesic { startPosition = 90°0'0.000"N,0°0'0.000"E (WGS84) -- , endPosition = 90°0'0.000"S,180°0'0.000"E (WGS84) -- , initialBearing = Just 0°0'0.000" -- , finalBearing = Just 180°0'0.000" -- , length = 20003.931458623km}) -- -- The Vincenty formula for the direct problem should always converge, however this function returns -- 'Nothing' if it would ever fail to do so (probably thus indicating a bug in the implementation). direct :: (Ellipsoidal a) => HorizontalPosition a -> Angle -> Length -> Maybe (Geodesic a) direct p1 b1 d | d == Length.zero = Just (Geodesic p1 p1 (Just b1) (Just b1) Length.zero) | otherwise = case rec of Nothing -> Nothing (Just (s, cosS, sinS, cos2S')) -> Just (Geodesic p1 p2 (Just b1) (Just b2) d) where x = sinU1 * sinS - cosU1 * cosS * cosAlpha1 lat2 = atan2 (sinU1 * cosS + cosU1 * sinS * cosAlpha1) ((1.0 - f) * sqrt (sinAlpha * sinAlpha + x * x)) lambda = atan2 (sinS * sinAlpha1) (cosU1 * cosS - sinU1 * sinS * cosAlpha1) _C = f / 16.0 * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)) _L = lambda - (1.0 - _C) * f * sinAlpha * (s + _C * sinS * (cos2S' + _C * cosS * (-1.0 + 2.0 * cos2S' * cos2S'))) lon2 = lon1 + _L b2 = Angle.normalise (Angle.radians (atan2 sinAlpha (-x))) (Angle.decimalDegrees 360.0) p2 = Geodetic.latLongPos' (Angle.radians lat2) (Angle.radians lon2) (Geodetic.model p1) where lat1 = Angle.toRadians . Geodetic.latitude $ p1 lon1 = Angle.toRadians . Geodetic.longitude $ p1 ell = surface . Geodetic.model $ p1 (a, b, f) = abf ell br1 = Angle.toRadians b1 cosAlpha1 = cos br1 sinAlpha1 = sin br1 (tanU1, cosU1, sinU1) = reducedLat lat1 f sigma1 = atan2 tanU1 cosAlpha1 -- angular distance on the sphere from the equator to p1 sinAlpha = cosU1 * sinAlpha1 -- alpha = azimuth of the geodesic at the equator cosSqAlpha = 1.0 - sinAlpha * sinAlpha uSq = cosSqAlpha * (a * a - b * b) / (b * b) _A = 1.0 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) _B = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) dm = Length.toMetres d sigma = dm / (b * _A) rec = directRec sigma1 dm _A _B b sigma 0 -- | @inverse p1 p2@ solves the inverse geodesic problem using Vicenty formula: __surface__ distance, -- and initial/final bearing between the geodesic line between positions @p1@ and @p2@. For example: -- -- >>> Geodesic.inverse (Geodetic.latLongPos 0 0 WGS84) (Geodetic.latLongPos 0.5 179.5 WGS84) -- Just (Geodesic { startPosition = 0°0'0.000"N,0°0'0.000"E 0.0m (WGS84) -- , endPosition = 0°30'0.000"N,179°30'0.000"E 0.0m (WGS84) -- , initialBearing = Just 25°40'18.742" -- , finalBearing = Just 154°19'37.507" -- , length = 19936.288578981km}) -- -- The Vincenty formula for the inverse problem can fail to converge for nearly antipodal points in which -- case this function returns 'Nothing'. For example: -- -- >>> Geodesic.inverse (Geodetic.latLongPos 0 0 WGS84) (Geodetic.latLongPos 0.5 179.7 WGS84) -- Nothing inverse :: (Ellipsoidal a) => HorizontalPosition a -> HorizontalPosition a -> Maybe (Geodesic a) inverse p1 p2 | p1 == p2 = Just (Geodesic p1 p2 Nothing Nothing Length.zero) | otherwise = case rec of Nothing -> Nothing (Just (cosL, sinL, s, cosS, sinS, sinSqS, cos2S', cosSqA)) -> Just (Geodesic p1 p2 (Just b1) (Just b2) d) where uSq = cosSqA * (a * a - b * b) / (b * b) _A = 1 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))) _B = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))) deltaSigma = _B * sinS * (cos2S' + _B / 4.0 * (cosS * (-1.0 + 2.0 * cos2S' * cos2S') - _B / 6.0 * cos2S' * (-3.0 + 4.0 * sinS * sinS) * (-3.0 + 4.0 * cos2S' * cos2S'))) d = Length.metres (b * _A * (s - deltaSigma)) a1R = if abs sinSqS < epsilon then 0.0 else atan2 (cosU2 * sinL) (cosU1 * sinU2 - sinU1 * cosU2 * cosL) a2R = if abs sinSqS < epsilon then pi else atan2 (cosU1 * sinL) (-sinU1 * cosU2 + cosU1 * sinU2 * cosL) b1 = Angle.normalise (Angle.radians a1R) (Angle.decimalDegrees 360.0) b2 = Angle.normalise (Angle.radians a2R) (Angle.decimalDegrees 360.0) where lat1 = Angle.toRadians . Geodetic.latitude $ p1 lon1 = Angle.toRadians . Geodetic.longitude $ p1 lat2 = Angle.toRadians . Geodetic.latitude $ p2 lon2 = Angle.toRadians . Geodetic.longitude $ p2 ell = surface . Geodetic.model $ p1 (a, b, f) = abf ell _L = lon2 - lon1 -- difference in longitude (_, cosU1, sinU1) = reducedLat lat1 f (_, cosU2, sinU2) = reducedLat lat2 f antipodal = abs _L > pi / 2.0 || abs (lat2 - lat1) > pi / 2.0 rec = inverseRec _L cosU1 sinU1 cosU2 sinU2 _L f antipodal 0 directRec :: Double -> Double -> Double -> Double -> Double -> Double -> Int -> Maybe (Double, Double, Double, Double) directRec sigma1 dist _A _B b sigma i | i == 100 = Nothing | abs (sigma - newSigma) <= 1e-12 = Just (newSigma, cosSigma, sinSigma, cos2Sigma') | otherwise = directRec sigma1 dist _A _B b newSigma (i + 1) where cos2Sigma' = cos (2 * sigma1 + sigma) sinSigma = sin sigma cosSigma = cos sigma deltaSigma = _B * sinSigma * (cos2Sigma' + _B / 4.0 * (cosSigma * (-1.0 + 2.0 * cos2Sigma' * cos2Sigma') - _B / 6.0 * cos2Sigma' * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2Sigma' * cos2Sigma'))) newSigma = dist / (b * _A) + deltaSigma inverseRec :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Bool -> Int -> Maybe (Double, Double, Double, Double, Double, Double, Double, Double) inverseRec lambda cosU1 sinU1 cosU2 sinU2 _L f antipodal i | i == 1000 = Nothing -- co-incident/antipodal points (falls back on λ/σ = L) | sinSqSigma < epsilon = Just (inverseFallback cosL sinL sinSqSigma antipodal) | iterationCheck > pi = Nothing | abs (lambda - newLambda) <= 1e-12 = Just (cosL, sinL, sigma, cosSigma, sinSigma, sinSqSigma, cos2Sigma', cosSqAlpha) | otherwise = inverseRec newLambda cosU1 sinU1 cosU2 sinU2 _L f antipodal (i + 1) where sinL = sin lambda cosL = cos lambda sinSqSigma = (cosU2 * sinL) * (cosU2 * sinL) + (cosU1 * sinU2 - sinU1 * cosU2 * cosL) * (cosU1 * sinU2 - sinU1 * cosU2 * cosL) sinSigma = sqrt sinSqSigma cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosL sigma = atan2 sinSigma cosSigma sinAlpha = cosU1 * cosU2 * sinL / sinSigma cosSqAlpha = 1 - sinAlpha * sinAlpha cos2Sigma' = if cosSqAlpha /= 0 then cosSigma - 2.0 * sinU1 * sinU2 / cosSqAlpha else 0 _C = f / 16.0 * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)) newLambda = _L + (1.0 - _C) * f * sinAlpha * (sigma + _C * sinSigma * (cos2Sigma' + _C * cosSigma * (-1.0 + 2.0 * cos2Sigma' * cos2Sigma'))) iterationCheck = if antipodal then abs newLambda - pi else abs newLambda inverseFallback :: Double -> Double -> Double -> Bool -> (Double, Double, Double, Double, Double, Double, Double, Double) inverseFallback cosL sinL sinSqSigma antipodal = (cosL, sinL, sigma, cosSigma, sinSigma, sinSqSigma, cos2Sigma', cosSqAlpha) where sigma = if antipodal then pi else 0 cosSigma = if antipodal then (-1) else 1 sinSigma = 0 cos2Sigma' = 1 cosSqAlpha = 1 -- | see Numeric.Limits epsilon :: Double epsilon = r where r = 1 - encodeFloat (m - 1) e (m, e) = decodeFloat (1 :: Double) reducedLat :: Double -> Double -> (Double, Double, Double) reducedLat lat f = (tanU, cosU, sinU) where tanU = (1.0 - f) * tan lat cosU = 1.0 / sqrt (1 + tanU * tanU) sinU = tanU * cosU abf :: Ellipsoid -> (Double, Double, Double) abf e = (Length.toMetres . equatorialRadius $ e, Length.toMetres . polarRadius $ e, flattening e)