module Geodetics.Geodetic (
   -- * Geodetic Coordinates
   Geodetic (..),
   readGroundPosition,
   toLocal,
   toWGS84,
   antipode,
   geometricalDistance,
   geometricalDistanceSq,
   groundDistance,
   properAngle,
   showAngle,
   showGeodeticLatLong,
   showGeodeticSignedDecimal,
   showGeodeticNSEWDecimal,
   showGeodeticDDDMMSS,
   -- * Earth Centred Earth Fixed Coordinates
   ECEF,
   geoToEarth,
   earthToGeo,
   -- * Re-exported for convenience
   WGS84 (..)
) where

import Data.Bool (bool)
import Data.Maybe
import Geodetics.Altitude
import Geodetics.Ellipsoids
import Geodetics.LatLongParser
import Text.ParserCombinators.ReadP
import Text.Printf

-- | 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
-- 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
-- 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
-- in the same place. However if different ellipsoids are used then
-- two co-ordinates with different numbers can still refer to the same
-- 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 Geodetic e = Geodetic {
   forall e. Geodetic e -> Double
latitude, forall e. Geodetic e -> Double
longitude :: Double,  -- ^ In radians.
   forall e. Geodetic e -> Double
geoAlt :: Double,  -- ^ In meters.
   forall e. Geodetic e -> e
ellipsoid :: e
}

instance (Ellipsoid e) => Show (Geodetic e) where
   show :: Geodetic e -> String
show Geodetic e
g =
     Geodetic e -> String
forall e. Geodetic e -> String
showGeodeticLatLong Geodetic e
g String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
", " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show (Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
g) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" " String -> ShowS
forall a. [a] -> [a] -> [a]
++ e -> String
forall a. Show a => a -> String
show (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)

-- | 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.
-- 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' W
--
-- * 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 :: forall e. Ellipsoid e => e -> String -> Maybe (Geodetic e)
readGroundPosition e
e String
str =
   case (((Double, Double), String) -> (Double, Double))
-> [((Double, Double), String)] -> [(Double, Double)]
forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double), String) -> (Double, Double)
forall a b. (a, b) -> a
fst ([((Double, Double), String)] -> [(Double, Double)])
-> [((Double, Double), String)] -> [(Double, Double)]
forall a b. (a -> b) -> a -> b
$ (((Double, Double), String) -> Bool)
-> [((Double, Double), String)] -> [((Double, Double), String)]
forall a. (a -> Bool) -> [a] -> [a]
filter (String -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null (String -> Bool)
-> (((Double, Double), String) -> String)
-> ((Double, Double), String)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Double, Double), String) -> String
forall a b. (a, b) -> b
snd) ([((Double, Double), String)] -> [((Double, Double), String)])
-> [((Double, Double), String)] -> [((Double, Double), String)]
forall a b. (a -> b) -> a -> b
$ ReadP (Double, Double) -> ReadS (Double, Double)
forall a. ReadP a -> ReadS a
readP_to_S ReadP (Double, Double)
latLong String
str of
      [(Double
lat,Double
long)] -> Geodetic e -> Maybe (Geodetic e)
forall a. a -> Maybe a
Just (Geodetic e -> Maybe (Geodetic e))
-> Geodetic e -> Maybe (Geodetic e)
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Geodetic e
forall a. HasAltitude a => a -> a
groundPosition (Geodetic e -> Geodetic e) -> Geodetic e -> Geodetic e
forall a b. (a -> b) -> a -> b
$ Geodetic
        { latitude :: Double
latitude = Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree,
          longitude :: Double
longitude = Double
long Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree,
          geoAlt :: Double
geoAlt = Double
0.0,
          ellipsoid :: e
ellipsoid = e
e
        }
      -- It appears incorrect to accept more than 1 interpretation since
      -- that would mean a pathological ambiguity in the parser. As stated,
      -- the accepted grammars do not overlap. Furthermore, should there
      -- exist several different ways to parse the input, it appears
      -- dangerous to arbitrarily choose one of them while disregarding the
      -- others, hence that case should be treated as a failure.
      [(Double, Double)]
_ -> Maybe (Geodetic e)
forall a. Maybe a
Nothing

-- | Show an angle as degrees, minutes and seconds to two decimal places.
showAngle :: Double -> String
showAngle :: Double -> String
showAngle Double
a
   | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
a        = String
"NaN"  -- Not an angle
   | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
a   = String
sign String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"Infinity"
   | Bool
otherwise      = String -> String -> Integer -> Integer -> Double -> String
forall r. PrintfType r => String -> r
printf String
"%s%d° %d′ %.2f″" String
sign  Integer
d Integer
m Double
s
   where
     sign :: String
sign = if Bool
isPositive then String
"" else String
"-"
     (Integer
d, Integer
m, Double
s, Bool
isPositive) = Double -> (Integer, Integer, Double, Bool)
radianToDegrees Double
a

-- | Show 'Geodetic' as a pair of degrees, minutes, and seconds. This is
-- similar to the 'Show' instance for 'Geodetic', except it does not include
-- the altitude and the ellipsoid.
--
-- @since 1.1.0
showGeodeticLatLong :: Geodetic e -> String
showGeodeticLatLong :: forall e. Geodetic e -> String
showGeodeticLatLong Geodetic e
x =
  String
-> Integer
-> Integer
-> Double
-> Char
-> Integer
-> Integer
-> Double
-> Char
-> String
forall r. PrintfType r => String -> r
printf
    String
"%d° %d′ %.2f″ %c, %d° %d′ %.2f″ %c"
    Integer
latD Integer
latM Double
latS (Char -> Char -> Bool -> Char
forall a. a -> a -> Bool -> a
bool Char
'S' Char
'N' Bool
isNorth)
    Integer
longD Integer
longM Double
longS (Char -> Char -> Bool -> Char
forall a. a -> a -> Bool -> a
bool Char
'W' Char
'E' Bool
isEast)
  where
    (Integer
latD, Integer
latM, Double
latS, Bool
isNorth) = Double -> (Integer, Integer, Double, Bool)
radianToDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
x)
    (Integer
longD, Integer
longM, Double
longS, Bool
isEast) = Double -> (Integer, Integer, Double, Bool)
radianToDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
x)

-- | Show 'Geodetic' as a pair of signed decimal degrees (5 decimal places
-- of precision), e.g. 34.52327, -46.23234.
--
-- @since 1.1.0
showGeodeticSignedDecimal :: Geodetic e -> String
showGeodeticSignedDecimal :: forall e. Geodetic e -> String
showGeodeticSignedDecimal Geodetic e
x =
  String -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf
    String
"%.5f, %.5f"
    (Double -> Double
forall {a}. Floating a => a -> a
toDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
x))
    (Double -> Double
forall {a}. Floating a => a -> a
toDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
x))
  where
    toDegrees :: a -> a
toDegrees a
r = a
r a -> a -> a
forall a. Num a => a -> a -> a
* a
180 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
forall a. Floating a => a
pi

-- | Show 'Geodetic' as a pair of decimal degrees NSEW, e.g. 34.52327N,
-- 46.23234W.
--
-- @since 1.1.0
showGeodeticNSEWDecimal :: Geodetic e -> String
showGeodeticNSEWDecimal :: forall e. Geodetic e -> String
showGeodeticNSEWDecimal Geodetic e
x =
  String -> Double -> Char -> Double -> Char -> String
forall r. PrintfType r => String -> r
printf
    String
"%.5f%c, %.5f%c"
    (Double -> Double
forall {a}. Floating a => a -> a
toAbsDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
x))
    (Char -> Char -> Double -> Char
forall {a} {p}. (Ord a, Num a) => p -> p -> a -> p
c Char
'N' Char
'S' (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
x))
    (Double -> Double
forall {a}. Floating a => a -> a
toAbsDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
x))
    (Char -> Char -> Double -> Char
forall {a} {p}. (Ord a, Num a) => p -> p -> a -> p
c Char
'E' Char
'W' (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
x))
  where
    toAbsDegrees :: a -> a
toAbsDegrees a
r = a -> a
forall a. Num a => a -> a
abs (a
r a -> a -> a
forall a. Num a => a -> a -> a
* a
180 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
forall a. Floating a => a
pi)
    c :: p -> p -> a -> p
c p
forPositive p
forNegative a
r =
      if a
r a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0
        then p
forNegative
        else p
forPositive

-- | Show 'Geodetic' as a pair of angles in the DDDMMSS format, e.g.
-- 343123.52N, 461356.43W.
--
-- @since 1.1.0
showGeodeticDDDMMSS ::
  -- | Use leading zeros
  Bool ->
  Geodetic e ->
  String
showGeodeticDDDMMSS :: forall e. Bool -> Geodetic e -> String
showGeodeticDDDMMSS Bool
useLeadingZeros Geodetic e
x =
  String
-> Integer
-> Integer
-> Double
-> Char
-> Integer
-> Integer
-> Double
-> Char
-> String
forall r. PrintfType r => String -> r
printf
    String
fmt
    Integer
latD Integer
latM Double
latS (Char -> Char -> Bool -> Char
forall a. a -> a -> Bool -> a
bool Char
'S' Char
'N' Bool
isNorth)
    Integer
longD Integer
longM Double
longS (Char -> Char -> Bool -> Char
forall a. a -> a -> Bool -> a
bool Char
'W' Char
'E' Bool
isEast)
  where
    fmt :: String
fmt =
      if Bool
useLeadingZeros
        then String
"%03d%02d%05.2f%c, %03d%02d%05.2f%c"
        else String
"%d%02d%05.2f%c, %d%02d%05.2f%c"
    (Integer
latD, Integer
latM, Double
latS, Bool
isNorth) = Double -> (Integer, Integer, Double, Bool)
radianToDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
x)
    (Integer
longD, Integer
longM, Double
longS, Bool
isEast) = Double -> (Integer, Integer, Double, Bool)
radianToDegrees (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
x)

-- | An internal helper function for 'showAngle' and 'showGeodeticDDDMMSS'.
radianToDegrees :: Double -> (Integer, Integer, Double, Bool)
radianToDegrees :: Double -> (Integer, Integer, Double, Bool)
radianToDegrees Double
a = (Integer
d, Integer
m, Double
s, Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0)
  where
    centisecs :: Integer
    centisecs :: Integer
centisecs = Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
arcsecond Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
100))
    (Integer
d, Integer
m1) = Integer
centisecs Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
360000
    (Integer
m, Integer
_) = Integer
m1 Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
6000   -- hundredths of arcsec per arcmin
    s :: Double
s = Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$
      (Double -> Double
forall a. Num a => a -> a
abs Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree Double -> Double -> Double
forall a. Num a => a -> a -> a
- Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
arcminute) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
arcsecond

instance (Ellipsoid e) => HasAltitude (Geodetic e) where
   altitude :: Geodetic e -> Double
altitude = Geodetic e -> Double
forall e. Geodetic e -> Double
geoAlt
   setAltitude :: Double -> Geodetic e -> Geodetic e
setAltitude Double
h Geodetic e
g = Geodetic e
g{geoAlt = h}



-- | The point on the Earth diametrically opposite the argument, with
-- the same altitude.
antipode :: Geodetic e -> Geodetic e
antipode :: forall e. Geodetic e -> Geodetic e
antipode Geodetic e
g = Double -> Double -> Double -> e -> Geodetic e
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
long (Geodetic e -> Double
forall e. Geodetic e -> Double
geoAlt Geodetic e
g) (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)
   where
      lat :: Double
lat = Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
g
      long' :: Double
long' = Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
g Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
180 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree
      long :: Double
long | Double
long' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0  = Double
long' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
360 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
degree
           | Bool
otherwise  = Double
long'



-- | Convert a geodetic coordinate into earth centered, relative to the
-- ellipsoid in use.
geoToEarth :: (Ellipsoid e) => Geodetic e -> ECEF
geoToEarth :: forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
geo = (
      (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslong,
      (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coslat Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinlong,
      (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinlat)
   where
      n :: Double
n = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
normal e
e (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      e :: e
e = Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
geo
      coslat :: Double
coslat = Double -> Double
forall {a}. Floating a => a -> a
cos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      coslong :: Double
coslong = Double -> Double
forall {a}. Floating a => a -> a
cos (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo
      sinlat :: Double
sinlat = Double -> Double
forall {a}. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
      sinlong :: Double
sinlong = Double -> Double
forall {a}. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo
      h :: Double
h = Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
geo


-- | Convert an earth centred coordinate into a geodetic coordinate on
-- the specified geoid.
--
-- Uses the closed form solution of H. Vermeille: Direct
-- transformation from geocentric coordinates to geodetic coordinates.
-- Journal of Geodesy Volume 76, Number 8 (2002), 451-454. Result is in the form
-- @(latitude, longitude, altitude)@.
earthToGeo :: (Ellipsoid e) => e -> ECEF -> (Double, Double, Double)
earthToGeo :: forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo e
e (Double
x,Double
y,Double
z) = (Double
phi, Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x, Double -> Double
forall {a}. Floating a => a -> a
sqrt (Double
l Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p2) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
norm)
   where
      -- Naming: numeric suffix inicates power. Hence x2 = x * x, x3 = x2 * x, etc.
      p2 :: Double
p2 = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y
      a :: Double
a = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius e
e
      a2 :: Double
a2 = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a
      e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
e
      e4 :: Double
e4 = Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e2
      zeta :: Double
zeta = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
e2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a2)
      rho :: Double
rho = (Double
p2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
zeta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e4) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6
      rho2 :: Double
rho2 = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho
      rho3 :: Double
rho3 = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho2
      s :: Double
s = Double
e4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zeta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a2)
      t :: Double
t = (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rho3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall {a}. Floating a => a -> a
sqrt (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho3))) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3) -- Cube root
      u :: Double
u = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rho2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t
      v :: Double
v = Double -> Double
forall {a}. Floating a => a -> a
sqrt (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zeta)
      w :: Double
w = Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
zeta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
      kappa :: Double
kappa = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall {a}. Floating a => a -> a
sqrt (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v)
      phi :: Double
phi = Double -> Double
forall {a}. Floating a => a -> a
atan (Double
kappa Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall {a}. Floating a => a -> a
sqrt Double
p2)
      norm :: Double
norm = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
normal e
e Double
phi
      l :: Double
l = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
norm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall {a}. Floating a => a -> a
sin Double
phi


-- | Convert a position from any geodetic to another one, assuming local altitude stays constant.
toLocal :: (Ellipsoid e1, Ellipsoid e2) => e2 -> Geodetic e1 -> Geodetic e2
toLocal :: forall e1 e2.
(Ellipsoid e1, Ellipsoid e2) =>
e2 -> Geodetic e1 -> Geodetic e2
toLocal e2
e2 Geodetic e1
g = Double -> Double -> Double -> e2 -> Geodetic e2
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
lon Double
alt e2
e2
   where
      alt :: Double
alt = Geodetic e1 -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e1
g
      (Double
lat, Double
lon, Double
_) = e2 -> ECEF -> ECEF
forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo e2
e2 (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Helmert -> ECEF -> ECEF
applyHelmert Helmert
h (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Geodetic e1 -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e1
g
      h :: Helmert
h = e1 -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert (Geodetic e1 -> e1
forall e. Geodetic e -> e
ellipsoid Geodetic e1
g) Helmert -> Helmert -> Helmert
forall a. Monoid a => a -> a -> a
`mappend` Helmert -> Helmert
inverseHelmert (e2 -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert e2
e2)

-- | Convert a position from any geodetic to WGS84, assuming local
-- altitude stays constant.
toWGS84 :: (Ellipsoid e) => Geodetic e -> Geodetic WGS84
toWGS84 :: forall e. Ellipsoid e => Geodetic e -> Geodetic WGS84
toWGS84 Geodetic e
g = Double -> Double -> Double -> WGS84 -> Geodetic WGS84
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic Double
lat Double
lon Double
alt WGS84
WGS84
   where
      alt :: Double
alt = Geodetic e -> Double
forall a. HasAltitude a => a -> Double
altitude Geodetic e
g
      (Double
lat, Double
lon, Double
_) = WGS84 -> ECEF -> ECEF
forall e. Ellipsoid e => e -> ECEF -> ECEF
earthToGeo WGS84
WGS84 (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Helmert -> ECEF -> ECEF
applyHelmert Helmert
h (ECEF -> ECEF) -> ECEF -> ECEF
forall a b. (a -> b) -> a -> b
$ Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g
      h :: Helmert
h = e -> Helmert
forall a. Ellipsoid a => a -> Helmert
helmert (Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
g)


-- | 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
-- the curvature of the earth.
geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
geometricalDistance :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistance Geodetic e
g1 Geodetic e
g2 = Double -> Double
forall {a}. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Geodetic e -> Double
forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq Geodetic e
g1 Geodetic e
g2

-- | The square of the absolute distance.
geometricalDistanceSq :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Double
geometricalDistanceSq Geodetic e
g1 Geodetic e
g2 = (Double
x1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
y1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
z1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
z2) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
   where
      (Double
x1,Double
y1,Double
z1) = Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g1
      (Double
x2,Double
y2,Double
z2) = Geodetic e -> ECEF
forall e. Ellipsoid e => Geodetic e -> ECEF
geoToEarth Geodetic e
g2


-- | The shortest ellipsoidal distance between two points on the
-- ground with reference to the same ellipsoid. Altitude is ignored.
--
-- The results are the distance between the points, the bearing of
-- the second point from the first, and (180 degrees - the bearing
-- of the first point from the second).
--
-- The algorithm can fail to converge where the arguments are near to
-- antipodal. In this case it returns @Nothing@.
--
-- Uses Vincenty's formula. \"Direct and inverse solutions of
-- geodesics on the ellipsoid with application of nested
-- equations\". T. Vincenty. Survey Review XXII 176, April
-- 1975. <http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf>
groundDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e ->
                  Maybe (Double, Double, Double)
groundDistance :: forall e. Ellipsoid e => Geodetic e -> Geodetic e -> Maybe ECEF
groundDistance Geodetic e
p1 Geodetic e
p2 = do
     ((Double, (Double, Double, Double, Double, Double))
_, (Double
lambda, (Double
cos2Alpha, Double
delta, Double
sinDelta, Double
cosDelta, Double
cos2DeltaM))) <-
       [((Double, (Double, Double, Double, Double, Double)),
  (Double, (Double, Double, Double, Double, Double)))]
-> Maybe
     ((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))
forall a. [a] -> Maybe a
listToMaybe ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> Maybe
      ((Double, (Double, Double, Double, Double, Double)),
       (Double, (Double, Double, Double, Double, Double))))
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> Maybe
     ((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))
forall a b. (a -> b) -> a -> b
$ (((Double, (Double, Double, Double, Double, Double)),
  (Double, (Double, Double, Double, Double, Double)))
 -> Bool)
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((Double, (Double, Double, Double, Double, Double)),
 (Double, (Double, Double, Double, Double, Double)))
-> Bool
forall {a} {b} {b}.
(Ord a, Fractional a) =>
((a, b), (a, b)) -> Bool
converging ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ Int
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a. Int -> [a] -> [a]
take Int
100 ([((Double, (Double, Double, Double, Double, Double)),
   (Double, (Double, Double, Double, Double, Double)))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ [(Double, (Double, Double, Double, Double, Double))]
-> [(Double, (Double, Double, Double, Double, Double))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Double, (Double, Double, Double, Double, Double))]
lambdas ([(Double, (Double, Double, Double, Double, Double))]
 -> [((Double, (Double, Double, Double, Double, Double)),
      (Double, (Double, Double, Double, Double, Double)))])
-> [(Double, (Double, Double, Double, Double, Double))]
-> [((Double, (Double, Double, Double, Double, Double)),
     (Double, (Double, Double, Double, Double, Double)))]
forall a b. (a -> b) -> a -> b
$ Int
-> [(Double, (Double, Double, Double, Double, Double))]
-> [(Double, (Double, Double, Double, Double, Double))]
forall a. Int -> [a] -> [a]
drop Int
1 [(Double, (Double, Double, Double, Double, Double))]
lambdas
     let
       uSq :: Double
uSq = Double
cos2Alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
aDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
bDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
       bigA :: Double
bigA = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSqDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
16384 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4096 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
768) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
320 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
175Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
uSq)))
       bigB :: Double
bigB =     Double
uSqDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
1024  Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
256  Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
128) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
uSq Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
74 Double -> Double -> Double
forall a. Num a => a -> a -> a
-  Double
47Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
uSq)))
       deltaDelta :: Double
deltaDelta =
         Double
bigB Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
+
                             Double
bigBDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaMDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
                                       Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bigBDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDeltaDouble -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3)
                                          Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3)))
       s :: Double
s = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
bigA Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
deltaDelta)
       alpha1 :: Double
alpha1 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall {a}. Floating a => a -> a
sin Double
lambda) (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 -> Double
forall {a}. Floating a => a -> a
cos Double
lambda)
       alpha2 :: Double
alpha2 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
cosU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall {a}. Floating a => a -> a
sin Double
lambda) (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 -> Double
forall {a}. Floating a => a -> a
cos Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinU1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosU2)
     ECEF -> Maybe ECEF
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double
s, Double
alpha1, Double
alpha2)
  where
    f :: Double
f = e -> Double
forall e. Ellipsoid e => e -> Double
flattening (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    a :: Double
a = e -> Double
forall e. Ellipsoid e => e -> Double
majorRadius (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    b :: Double
b = e -> Double
forall e. Ellipsoid e => e -> Double
minorRadius (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
p1
    l :: Double
l = Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
p1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
p2
    u1 :: Double
u1 = Double -> Double
forall {a}. Floating a => a -> a
atan ((Double
1Double -> 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 (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
p1))
    u2 :: Double
u2 = Double -> Double
forall {a}. Floating a => a -> a
atan ((Double
1Double -> 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 (Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
p2))
    sinU1 :: Double
sinU1 = Double -> Double
forall {a}. Floating a => a -> a
sin Double
u1
    cosU1 :: Double
cosU1 = Double -> Double
forall {a}. Floating a => a -> a
cos Double
u1
    sinU2 :: Double
sinU2 = Double -> Double
forall {a}. Floating a => a -> a
sin Double
u2
    cosU2 :: Double
cosU2 = Double -> Double
forall {a}. Floating a => a -> a
cos Double
u2

    nextLambda :: Double -> (Double, (Double, Double, Double, Double, Double))
nextLambda Double
lambda = (Double
lambda1, (Double
cos2Alpha, Double
delta, Double
sinDelta, Double
cosDelta, Double
cos2DeltaM))
      where
        sinLambda :: Double
sinLambda = Double -> Double
forall {a}. Floating a => a -> a
sin Double
lambda
        cosLambda :: Double
cosLambda = Double -> Double
forall {a}. Floating a => a -> a
cos Double
lambda
        sinDelta :: Double
sinDelta = Double -> Double
forall {a}. Floating a => a -> a
sqrt ((Double
cosU2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLambda) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 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
cosLambda) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2)
        cosDelta :: Double
cosDelta = 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
cosLambda
        delta :: Double
delta = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
sinDelta Double
cosDelta
        sinAlpha :: Double
sinAlpha = if Double
sinDelta Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else 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
sinLambda Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sinDelta
        cos2Alpha :: Double
cos2Alpha = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinAlpha Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2
        cos2DeltaM :: Double
cos2DeltaM = if Double
cos2Alpha Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                     then Double
0
                     else Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 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
cos2Alpha
        c :: Double
c = (Double
fDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
16) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2Alpha))
        lambda1 :: Double
lambda1 = Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1Double -> 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
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinDelta
                     Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
cos2DeltaM Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosDelta Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cos2DeltaM Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)))
    lambdas :: [(Double, (Double, Double, Double, Double, Double))]
lambdas = ((Double, (Double, Double, Double, Double, Double))
 -> (Double, (Double, Double, Double, Double, Double)))
-> (Double, (Double, Double, Double, Double, Double))
-> [(Double, (Double, Double, Double, Double, Double))]
forall a. (a -> a) -> a -> [a]
iterate (Double -> (Double, (Double, Double, Double, Double, Double))
nextLambda (Double -> (Double, (Double, Double, Double, Double, Double)))
-> ((Double, (Double, Double, Double, Double, Double)) -> Double)
-> (Double, (Double, Double, Double, Double, Double))
-> (Double, (Double, Double, Double, Double, Double))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, (Double, Double, Double, Double, Double)) -> Double
forall a b. (a, b) -> a
fst) (Double
l, (Double, Double, Double, Double, Double)
forall a. HasCallStack => a
undefined)
    converging :: ((a, b), (a, b)) -> Bool
converging ((a
l1,b
_),(a
l2,b
_)) = a -> a
forall a. Num a => a -> a
abs (a
l1 a -> a -> a
forall a. Num a => a -> a -> a
- a
l2) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
1e-14


-- | Add or subtract multiples of 2*pi so that for all @t@, @-pi < properAngle t < pi@.
properAngle :: Double -> Double
properAngle :: Double -> Double
properAngle Double
t
   | Double
r1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double -> Double
forall a. Num a => a -> a
negate Double
forall a. Floating a => a
pi    = Double
r1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
pi2
   | Double
r1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
forall a. Floating a => a
pi            = Double
r1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
pi2
   | Bool
otherwise          = Double
r1
   where
      pf :: Double -> (Int, Double)
      pf :: Double -> (Int, Double)
pf = Double -> (Int, Double)
forall b. Integral b => Double -> (b, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction  -- Shut up GHC warning about defaulting to Integer.
      (Int
_,Double
r) = Double -> (Int, Double)
pf (Double
tDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
pi2)
      r1 :: Double
r1 = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
pi2
      pi2 :: Double
pi2 = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2