module Data.Geo.Jord.Geodesic
(
Geodesic
, startPosition
, endPosition
, initialBearing
, finalBearing
, length
, 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)
data Geodesic a =
Geodesic
{ Geodesic a -> HorizontalPosition a
startPosition :: HorizontalPosition a
, Geodesic a -> HorizontalPosition a
endPosition :: HorizontalPosition a
, Geodesic a -> Maybe Angle
initialBearing :: Maybe Angle
, Geodesic a -> Maybe Angle
finalBearing :: Maybe Angle
, Geodesic a -> Length
length :: Length
}
deriving (Geodesic a -> Geodesic a -> Bool
(Geodesic a -> Geodesic a -> Bool)
-> (Geodesic a -> Geodesic a -> Bool) -> Eq (Geodesic a)
forall a. Model a => Geodesic a -> Geodesic a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Geodesic a -> Geodesic a -> Bool
$c/= :: forall a. Model a => Geodesic a -> Geodesic a -> Bool
== :: Geodesic a -> Geodesic a -> Bool
$c== :: forall a. Model a => Geodesic a -> Geodesic a -> Bool
Eq, Int -> Geodesic a -> ShowS
[Geodesic a] -> ShowS
Geodesic a -> String
(Int -> Geodesic a -> ShowS)
-> (Geodesic a -> String)
-> ([Geodesic a] -> ShowS)
-> Show (Geodesic a)
forall a. Model a => Int -> Geodesic a -> ShowS
forall a. Model a => [Geodesic a] -> ShowS
forall a. Model a => Geodesic a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Geodesic a] -> ShowS
$cshowList :: forall a. Model a => [Geodesic a] -> ShowS
show :: Geodesic a -> String
$cshow :: forall a. Model a => Geodesic a -> String
showsPrec :: Int -> Geodesic a -> ShowS
$cshowsPrec :: forall a. Model a => Int -> Geodesic a -> ShowS
Show)
direct :: (Ellipsoidal a) => HorizontalPosition a -> Angle -> Length -> Maybe (Geodesic a)
direct :: HorizontalPosition a -> Angle -> Length -> Maybe (Geodesic a)
direct HorizontalPosition a
p1 Angle
b1 Length
d
| Length
d Length -> Length -> Bool
forall a. Eq a => a -> a -> Bool
== Length
Length.zero = Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p1 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) Length
Length.zero)
| Bool
otherwise =
case Maybe (Double, Double, Double, Double)
rec of
Maybe (Double, Double, Double, Double)
Nothing -> Maybe (Geodesic a)
forall a. Maybe a
Nothing
(Just (Double
s, Double
cosS, Double
sinS, Double
cos2S')) -> Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b2) Length
d)
where x :: Double
x = Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1
lat2 :: Double
lat2 =
Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2
(Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1)
((Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
lambda :: Double
lambda = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha1) (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosAlpha1)
_C :: Double
_C = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha))
_L :: Double
_L =
Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
-
(Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
_C) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S')))
lon2 :: Double
lon2 = Double
lon1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_L
b2 :: Angle
b2 =
Angle -> Angle -> Angle
Angle.normalise
(Double -> Angle
Angle.radians (Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinAlpha (-Double
x)))
(Double -> Angle
Angle.decimalDegrees Double
360.0)
p2 :: HorizontalPosition a
p2 =
Angle -> Angle -> a -> HorizontalPosition a
forall a. Model a => Angle -> Angle -> a -> HorizontalPosition a
Geodetic.latLongPos'
(Double -> Angle
Angle.radians Double
lat2)
(Double -> Angle
Angle.radians Double
lon2)
(HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
p1)
where
lat1 :: Double
lat1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
lon1 :: Double
lon1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
ell :: Ellipsoid
ell = a -> Ellipsoid
forall a. Model a => a -> Ellipsoid
surface (a -> Ellipsoid)
-> (HorizontalPosition a -> a) -> HorizontalPosition a -> Ellipsoid
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model (HorizontalPosition a -> Ellipsoid)
-> HorizontalPosition a -> Ellipsoid
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
(Double
a, Double
b, Double
f) = Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
ell
br1 :: Double
br1 = Angle -> Double
Angle.toRadians Angle
b1
cosAlpha1 :: Double
cosAlpha1 = Double -> Double
forall a. Floating a => a -> a
cos Double
br1
sinAlpha1 :: Double
sinAlpha1 = Double -> Double
forall a. Floating a => a -> a
sin Double
br1
(Double
tanU1, Double
cosU1, Double
sinU1) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat1 Double
f
sigma1 :: Double
sigma1 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
tanU1 Double
cosAlpha1
sinAlpha :: Double
sinAlpha = Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha1
cosSqAlpha :: Double
cosSqAlpha = Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha
uSq :: Double
uSq = Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
_A :: Double
_A = Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16384.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
768.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
320.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
_B :: Double
_B = Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1024.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
128.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
74.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
47.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
dm :: Double
dm = Length -> Double
Length.toMetres Length
d
sigma :: Double
sigma = Double
dm Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A)
rec :: Maybe (Double, Double, Double, Double)
rec = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dm Double
_A Double
_B Double
b Double
sigma Int
0
inverse :: (Ellipsoidal a) => HorizontalPosition a -> HorizontalPosition a -> Maybe (Geodesic a)
inverse :: HorizontalPosition a -> HorizontalPosition a -> Maybe (Geodesic a)
inverse HorizontalPosition a
p1 HorizontalPosition a
p2
| HorizontalPosition a
p1 HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== HorizontalPosition a
p2 = Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 Maybe Angle
forall a. Maybe a
Nothing Maybe Angle
forall a. Maybe a
Nothing Length
Length.zero)
| Bool
otherwise =
case Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
rec of
Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
Nothing -> Maybe (Geodesic a)
forall a. Maybe a
Nothing
(Just (Double
cosL, Double
sinL, Double
s, Double
cosS, Double
sinS, Double
sinSqS, Double
cos2S', Double
cosSqA)) ->
Geodesic a -> Maybe (Geodesic a)
forall a. a -> Maybe a
Just (HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
forall a.
HorizontalPosition a
-> HorizontalPosition a
-> Maybe Angle
-> Maybe Angle
-> Length
-> Geodesic a
Geodesic HorizontalPosition a
p1 HorizontalPosition a
p2 (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b1) (Angle -> Maybe Angle
forall a. a -> Maybe a
Just Angle
b2) Length
d)
where uSq :: Double
uSq = Double
cosSqA Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
_A :: Double
_A =
Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+
Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16384.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
768.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
320.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
_B :: Double
_B = Double
uSq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1024.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
128.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
74.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
47.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
deltaSigma :: Double
deltaSigma =
Double
_B Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
+
Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
cosS Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S') Double -> Double -> Double
forall a. Num a => a -> a -> a
-
Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinS) Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2S')))
d :: Length
d = Double -> Length
Length.metres (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
deltaSigma))
a1R :: Double
a1R =
if Double -> Double
forall a. Num a => a -> a
abs Double
sinSqS Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
then Double
0.0
else Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL)
a2R :: Double
a2R =
if Double -> Double
forall a. Num a => a -> a
abs Double
sinSqS Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
then Double
forall a. Floating a => a
pi
else Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) (-Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL)
b1 :: Angle
b1 = Angle -> Angle -> Angle
Angle.normalise (Double -> Angle
Angle.radians Double
a1R) (Double -> Angle
Angle.decimalDegrees Double
360.0)
b2 :: Angle
b2 = Angle -> Angle -> Angle
Angle.normalise (Double -> Angle
Angle.radians Double
a2R) (Double -> Angle
Angle.decimalDegrees Double
360.0)
where
lat1 :: Double
lat1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
lon1 :: Double
lon1 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
lat2 :: Double
lat2 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p2
lon2 :: Double
lon2 = Angle -> Double
Angle.toRadians (Angle -> Double)
-> (HorizontalPosition a -> Angle)
-> HorizontalPosition a
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude (HorizontalPosition a -> Double) -> HorizontalPosition a -> Double
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p2
ell :: Ellipsoid
ell = a -> Ellipsoid
forall a. Model a => a -> Ellipsoid
surface (a -> Ellipsoid)
-> (HorizontalPosition a -> a) -> HorizontalPosition a -> Ellipsoid
forall b c a. (b -> c) -> (a -> b) -> a -> c
. HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model (HorizontalPosition a -> Ellipsoid)
-> HorizontalPosition a -> Ellipsoid
forall a b. (a -> b) -> a -> b
$ HorizontalPosition a
p1
(Double
a, Double
b, Double
f) = Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
ell
_L :: Double
_L = Double
lon2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lon1
(Double
_, Double
cosU1, Double
sinU1) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat1 Double
f
(Double
_, Double
cosU2, Double
sinU2) = Double -> Double -> (Double, Double, Double)
reducedLat Double
lat2 Double
f
antipodal :: Bool
antipodal = Double -> Double
forall a. Num a => a -> a
abs Double
_L Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs (Double
lat2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lat1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0
rec :: Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
rec = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
_L Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal Int
0
directRec ::
Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec :: Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dist Double
_A Double
_B Double
b Double
sigma Int
i
| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
100 = Maybe (Double, Double, Double, Double)
forall a. Maybe a
Nothing
| Double -> Double
forall a. Num a => a -> a
abs (Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
newSigma) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-12 = (Double, Double, Double, Double)
-> Maybe (Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
newSigma, Double
cosSigma, Double
sinSigma, Double
cos2Sigma')
| Bool
otherwise = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Int
-> Maybe (Double, Double, Double, Double)
directRec Double
sigma1 Double
dist Double
_A Double
_B Double
b Double
newSigma (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
where
cos2Sigma' :: Double
cos2Sigma' = Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sigma1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sigma)
sinSigma :: Double
sinSigma = Double -> Double
forall a. Floating a => a -> a
sin Double
sigma
cosSigma :: Double
cosSigma = Double -> Double
forall a. Floating a => a -> a
cos Double
sigma
deltaSigma :: Double
deltaSigma =
Double
_B Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
+
Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma') Double -> Double -> Double
forall a. Num a => a -> a -> a
-
Double
_B Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma) Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(-Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma')))
newSigma :: Double
newSigma = Double
dist Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
_A) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
deltaSigma
inverseRec ::
Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe (Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec :: Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
lambda Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal Int
i
| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1000 = Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
forall a. Maybe a
Nothing
| Double
sinSqSigma Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon = (Double, Double, Double, Double, Double, Double, Double, Double)
-> Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
-> Double
-> Double
-> Bool
-> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback Double
cosL Double
sinL Double
sinSqSigma Bool
antipodal)
| Double
iterationCheck Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi = Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
forall a. Maybe a
Nothing
| Double -> Double
forall a. Num a => a -> a
abs (Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
newLambda) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-12 =
(Double, Double, Double, Double, Double, Double, Double, Double)
-> Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
forall a. a -> Maybe a
Just (Double
cosL, Double
sinL, Double
sigma, Double
cosSigma, Double
sinSigma, Double
sinSqSigma, Double
cos2Sigma', Double
cosSqAlpha)
| Bool
otherwise = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Bool
-> Int
-> Maybe
(Double, Double, Double, Double, Double, Double, Double, Double)
inverseRec Double
newLambda Double
cosU1 Double
sinU1 Double
cosU2 Double
sinU2 Double
_L Double
f Bool
antipodal (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
where
sinL :: Double
sinL = Double -> Double
forall a. Floating a => a -> a
sin Double
lambda
cosL :: Double
cosL = Double -> Double
forall a. Floating a => a -> a
cos Double
lambda
sinSqSigma :: Double
sinSqSigma =
(Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL) Double -> Double -> Double
forall a. Num a => a -> a -> a
+
(Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL)
sinSigma :: Double
sinSigma = Double -> Double
forall a. Floating a => a -> a
sqrt Double
sinSqSigma
cosSigma :: Double
cosSigma = Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosL
sigma :: Double
sigma = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinSigma Double
cosSigma
sinAlpha :: Double
sinAlpha = Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinL Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sinSigma
cosSqAlpha :: Double
cosSqAlpha = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha
cos2Sigma' :: Double
cos2Sigma' =
if Double
cosSqAlpha Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0
then Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinU2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
cosSqAlpha
else Double
0
_C :: Double
_C = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
16.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSqAlpha))
newLambda :: Double
newLambda =
Double
_L Double -> Double -> Double
forall a. Num a => a -> a -> a
+
(Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
_C) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinAlpha Double -> Double -> Double
forall a. Num a => a -> a -> a
*
(Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
+
Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
_C Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosSigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Sigma')))
iterationCheck :: Double
iterationCheck =
if Bool
antipodal
then Double -> Double
forall a. Num a => a -> a
abs Double
newLambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
pi
else Double -> Double
forall a. Num a => a -> a
abs Double
newLambda
inverseFallback ::
Double
-> Double
-> Double
-> Bool
-> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback :: Double
-> Double
-> Double
-> Bool
-> (Double, Double, Double, Double, Double, Double, Double, Double)
inverseFallback Double
cosL Double
sinL Double
sinSqSigma Bool
antipodal =
(Double
cosL, Double
sinL, Double
sigma, Double
cosSigma, Double
sinSigma, Double
sinSqSigma, Double
cos2Sigma', Double
cosSqAlpha)
where
sigma :: Double
sigma =
if Bool
antipodal
then Double
forall a. Floating a => a
pi
else Double
0
cosSigma :: Double
cosSigma =
if Bool
antipodal
then (-Double
1)
else Double
1
sinSigma :: Double
sinSigma = Double
0
cos2Sigma' :: Double
cos2Sigma' = Double
1
cosSqAlpha :: Double
cosSqAlpha = Double
1
epsilon :: Double
epsilon :: Double
epsilon = Double
r
where
r :: Double
r = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat (Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) Int
e
(Integer
m, Int
e) = Double -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat (Double
1 :: Double)
reducedLat :: Double -> Double -> (Double, Double, Double)
reducedLat :: Double -> Double -> (Double, Double, Double)
reducedLat Double
lat Double
f = (Double
tanU, Double
cosU, Double
sinU)
where
tanU :: Double
tanU = (Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan Double
lat
cosU :: Double
cosU = Double
1.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
tanU Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanU)
sinU :: Double
sinU = Double
tanU Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU
abf :: Ellipsoid -> (Double, Double, Double)
abf :: Ellipsoid -> (Double, Double, Double)
abf Ellipsoid
e = (Length -> Double
Length.toMetres (Length -> Double) -> (Ellipsoid -> Length) -> Ellipsoid -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ellipsoid -> Length
equatorialRadius (Ellipsoid -> Double) -> Ellipsoid -> Double
forall a b. (a -> b) -> a -> b
$ Ellipsoid
e, Length -> Double
Length.toMetres (Length -> Double) -> (Ellipsoid -> Length) -> Ellipsoid -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ellipsoid -> Length
polarRadius (Ellipsoid -> Double) -> Ellipsoid -> Double
forall a b. (a -> b) -> a -> b
$ Ellipsoid
e, Ellipsoid -> Double
flattening Ellipsoid
e)