{- |
The following is based on equations in Section 1.4.7.1 in 
OGP Surveying and Positioning Guidance Note number 7, part 2 – August 2006
<http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf>
-}

{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
module Geodetics.Stereographic (
   GridStereo (gridTangent, gridOrigin, gridScale),
   mkGridStereo
) where

import Geodetics.Ellipsoids
import Geodetics.Geodetic
import Geodetics.Grid

import qualified Data.Stream as Stream

-- | A stereographic projection with its origin at an arbitrary point on Earth, other than the poles.
data GridStereo e = GridStereo {
      forall e. GridStereo e -> Geodetic e
gridTangent :: Geodetic e, -- ^ Point where the plane of projection touches the ellipsoid. Often known as the Natural Origin.
      forall e. GridStereo e -> GridOffset
gridOrigin :: GridOffset,  -- ^ Grid position of the tangent point. Often known as the False Origin.
      forall e. GridStereo e -> Double
gridScale :: Double, -- ^ Scaling factor that balances the distortion between the center and the edges.
                                         -- Should be slightly less than unity.
      
      -- Memoised parameters derived from the tangent point.
      forall e. GridStereo e -> Double
gridR :: Double,
      forall e. GridStereo e -> Double
gridN, forall e. GridStereo e -> Double
gridC, forall e. GridStereo e -> Double
gridSin, forall e. GridStereo e -> Double
gridCos :: Double,
      forall e. GridStereo e -> Double
gridLatC :: Double,
      forall e. GridStereo e -> Double
gridG, forall e. GridStereo e -> Double
gridH :: Double
   } deriving (Int -> GridStereo e -> ShowS
[GridStereo e] -> ShowS
GridStereo e -> String
(Int -> GridStereo e -> ShowS)
-> (GridStereo e -> String)
-> ([GridStereo e] -> ShowS)
-> Show (GridStereo e)
forall e. Ellipsoid e => Int -> GridStereo e -> ShowS
forall e. Ellipsoid e => [GridStereo e] -> ShowS
forall e. Ellipsoid e => GridStereo e -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall e. Ellipsoid e => Int -> GridStereo e -> ShowS
showsPrec :: Int -> GridStereo e -> ShowS
$cshow :: forall e. Ellipsoid e => GridStereo e -> String
show :: GridStereo e -> String
$cshowList :: forall e. Ellipsoid e => [GridStereo e] -> ShowS
showList :: [GridStereo e] -> ShowS
Show)
   
-- | Create a stereographic projection. The tangency point must not be one of the poles.  
mkGridStereo :: (Ellipsoid e) => Geodetic e -> GridOffset -> Double -> GridStereo e
mkGridStereo :: forall e.
Ellipsoid e =>
Geodetic e -> GridOffset -> Double -> GridStereo e
mkGridStereo Geodetic e
tangent GridOffset
origin Double
scale = GridStereo {
      gridTangent :: Geodetic e
gridTangent = Geodetic e
tangent,
      gridOrigin :: GridOffset
gridOrigin = GridOffset
origin,
      gridScale :: Double
gridScale = Double
scale,
      gridR :: Double
gridR = Double
r,
      gridN :: Double
gridN = Double
n,
      gridC :: Double
gridC = Double
c,
      gridSin :: Double
gridSin = Double
sinLatC1,
      gridCos :: Double
gridCos = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLatC1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLatC1,
      gridLatC :: Double
gridLatC = Double -> Double
forall a. Floating a => a -> a
asin Double
sinLatC1,
      gridG :: Double
gridG = Double
g,
      gridH :: Double
gridH = Double
h
   }
   where 
      -- The reference seems to use χO to refer to two slightly different values. 
      -- Here these will be called LatC0 and LatC1.
      ellipse :: e
ellipse = Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
tangent
      op :: Num a => a -> a    -- Values of longitude, tangent longitude, E and N
      op :: forall a. Num a => a -> a
op = if Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
tangent Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then a -> a
forall a. Num a => a -> a
negate else a -> a
forall a. a -> a
id  -- must be negated in the southern hemisphere.
      lat0 :: Double
lat0 = Double -> Double
forall a. Num a => a -> a
op (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
tangent
      sinLat0 :: Double
sinLat0 = Double -> Double
forall a. Floating a => a -> a
sin Double
lat0
      e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 e
ellipse
      e :: Double
e = Double -> Double
forall a. Floating a => a -> a
sqrt Double
e2
      r :: Double
r = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
meridianRadius e
ellipse Double
lat0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
primeVerticalRadius e
ellipse Double
lat0
      n :: Double
n = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ 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
cos Double
lat0 Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_4)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2))
      s1 :: Double
s1 = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sinLat0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLat0)
      s2 :: Double
s2 = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLat0)
      w1 :: Double
w1 = (Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
e) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
n
      sinLatC0 :: Double
sinLatC0 = (Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
      c :: Double
c = ((Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sin Double
lat0) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLatC0)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sin Double
lat0) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sinLatC0))
      w2 :: Double
w2 = Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1
      sinLatC1 :: Double
sinLatC1 = (Double
w2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
w2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
      g :: Double
g = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
scale Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
latC1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
      h :: Double
h = Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
scale Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan Double
latC1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g
      latC1 :: Double
latC1 = Double -> Double
forall a. Floating a => a -> a
asin Double
sinLatC1
      

instance (Ellipsoid e) => GridClass (GridStereo e) e where
   toGrid :: GridStereo e -> Geodetic e -> GridPoint (GridStereo e)
toGrid GridStereo e
grid Geodetic e
geo = GridOffset -> GridPoint (GridStereo e) -> GridPoint (GridStereo e)
forall g. GridOffset -> GridPoint g -> GridPoint g
applyOffset (GridStereo e -> GridOffset
forall e. GridStereo e -> GridOffset
gridOrigin GridStereo e
grid) (GridPoint (GridStereo e) -> GridPoint (GridStereo e))
-> GridPoint (GridStereo e) -> GridPoint (GridStereo e)
forall a b. (a -> b) -> a -> b
$ Double
-> Double -> Double -> GridStereo e -> GridPoint (GridStereo e)
forall r. Double -> Double -> Double -> r -> GridPoint r
GridPoint Double
east Double
north (Geodetic e -> Double
forall e. Geodetic e -> Double
geoAlt Geodetic e
geo) GridStereo e
grid
      where
         op :: Num a => a -> a    -- Values of longitude, tangent longitude, E and N
         op :: forall a. Num a => a -> a
op = if Geodetic e -> Double
forall e. Geodetic e -> Double
latitude (GridStereo e -> Geodetic e
forall e. GridStereo e -> Geodetic e
gridTangent GridStereo e
grid) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then a -> a
forall a. Num a => a -> a
negate else a -> a
forall a. a -> a
id  -- must be negated in the southern hemisphere.
         sinLatC :: Double
sinLatC = (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
         cosLatC :: Double
cosLatC = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sinLatC
         longC :: Double
longC = GridStereo e -> Double
forall e. GridStereo e -> Double
gridN GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
op (Geodetic e -> Double
forall e. Geodetic e -> Double
longitude Geodetic e
geo) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
long0) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
long0
         w :: Double
w = GridStereo e -> Double
forall e. GridStereo e -> Double
gridC GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
sA Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sB Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
e) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** GridStereo e -> Double
forall e. GridStereo e -> Double
gridN GridStereo e
grid
         sA :: Double
sA = (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sinLat) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLat)
         sB :: Double
sB = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sinLat) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sinLat)
         sinLat :: Double
sinLat = Double -> Double
forall a. Floating a => a -> a
sin (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Num a => a -> a
op (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
latitude Geodetic e
geo
         e :: Double
e = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid Geodetic e
geo
         long0 :: Double
long0 = Double -> Double
forall a. Num a => a -> a
op (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude (Geodetic e -> Double) -> Geodetic e -> Double
forall a b. (a -> b) -> a -> b
$ GridStereo e -> Geodetic e
forall e. GridStereo e -> Geodetic e
gridTangent GridStereo e
grid
         b :: Double
b = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sinLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridSin GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cosLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridCos GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
longC Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
long0)
         east :: Double
east = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridR GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridScale GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
longC Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
long0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b
         north :: Double
north = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridR GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridScale GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
sinLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridCos GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cosLatC Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridSin GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
longC Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
long0)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b
   
   fromGrid :: GridPoint (GridStereo e) -> Geodetic e
fromGrid GridPoint (GridStereo e)
gp = 
      {- trace (    -- Remove comment brackets for debugging.
         "fromGrid values:\n   i = " ++ show i ++ "\n   j = " ++ show j ++
         "\n   longC = " ++ show longC ++ "\n   long = " ++ show long ++
         "\n   latC = " ++ show latC ++
         "\n   lat1 = " ++ show lat1 ++ "\n   latN = " ++ show latN ) $ -}
         Double -> Double -> Double -> e -> Geodetic e
forall e. Double -> Double -> Double -> e -> Geodetic e
Geodetic (Double -> Double
forall a. Num a => a -> a
op Double
latN) (Double -> Double
forall a. Num a => a -> a
op Double
long) Double
height (e -> Geodetic e) -> e -> Geodetic e
forall a b. (a -> b) -> a -> b
$ GridStereo e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridStereo e
grid
      where
         op :: Num a => a -> a                   -- Values of longitude, tangent longitude, E and N
         op :: forall a. Num a => a -> a
op = if Geodetic e -> Double
forall e. Geodetic e -> Double
latitude (GridStereo e -> Geodetic e
forall e. GridStereo e -> Geodetic e
gridTangent GridStereo e
grid) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then a -> a
forall a. Num a => a -> a
negate else a -> a
forall a. a -> a
id  -- must be negated in the southern hemisphere.
         GridPoint Double
east Double
north Double
height GridStereo e
_ = GridOffset -> GridPoint (GridStereo e) -> GridPoint (GridStereo e)
forall g. GridOffset -> GridPoint g -> GridPoint g
applyOffset (GridOffset -> GridOffset
offsetNegate (GridOffset -> GridOffset) -> GridOffset -> GridOffset
forall a b. (a -> b) -> a -> b
$ GridStereo e -> GridOffset
forall e. GridStereo e -> GridOffset
gridOrigin GridStereo e
grid) GridPoint (GridStereo e)
gp
         east' :: Double
east' = Double
east
         north' :: Double
north' = Double
north
         grid :: GridStereo e
grid = GridPoint (GridStereo e) -> GridStereo e
forall r. GridPoint r -> r
gridBasis GridPoint (GridStereo e)
gp
         long0 :: Double
long0 = Double -> Double
forall a. Num a => a -> a
op (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Geodetic e -> Double
forall e. Geodetic e -> Double
longitude (Geodetic e -> Double) -> Geodetic e -> Double
forall a b. (a -> b) -> a -> b
$ GridStereo e -> Geodetic e
forall e. GridStereo e -> Geodetic e
gridTangent GridStereo e
grid
         i :: Double
i = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
east' (GridStereo e -> Double
forall e. GridStereo e -> Double
gridH GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
north')
         j :: Double
j = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
east' (GridStereo e -> Double
forall e. GridStereo e -> Double
gridG GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
north') Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
i
         latC :: Double
latC = GridStereo e -> Double
forall e. GridStereo e -> Double
gridLatC GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
north' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
east' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Double
jDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)) (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridR GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridScale GridStereo e
grid)
         longC :: Double
longC = Double
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
long0
         sinLatC :: Double
sinLatC = Double -> Double
forall a. Floating a => a -> a
sin Double
latC
         long :: Double
long = (Double
longC Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
long0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ GridStereo e -> Double
forall e. GridStereo e -> Double
gridN GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
long0
         isoLat :: Double
isoLat = Double -> Double
forall a. Floating a => a -> a
log ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sinLatC) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (GridStereo e -> Double
forall e. GridStereo e -> Double
gridC GridStereo e
grid Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinLatC))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* GridStereo e -> Double
forall e. GridStereo e -> Double
gridN GridStereo e
grid)
         lat1 :: Double
lat1 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
atan (Double -> Double
forall a. Floating a => a -> a
exp Double
isoLat) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2
         next :: Double -> Double
next Double
lat = Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
isoN Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
isoLat) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
lat Double -> Double -> Double
forall a. Num a => a -> a -> a
* (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
sin Double
lat Double -> 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
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
e2)
            where isoN :: Double
isoN = e -> Double -> Double
forall e. Ellipsoid e => e -> Double -> Double
isometricLatitude (GridStereo e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridStereo e
grid) Double
lat
                  e2 :: Double
e2 = e -> Double
forall e. Ellipsoid e => e -> Double
eccentricity2 (e -> Double) -> e -> Double
forall a b. (a -> b) -> a -> b
$ GridStereo e -> e
forall r e. GridClass r e => r -> e
gridEllipsoid GridStereo e
grid
         lats :: Stream Double
lats = (Double -> Double) -> Double -> Stream Double
forall a. (a -> a) -> a -> Stream a
Stream.iterate Double -> Double
next Double
lat1
         latN :: Double
latN = (Double, Double) -> Double
forall a b. (a, b) -> b
snd ((Double, Double) -> Double) -> (Double, Double) -> Double
forall a b. (a -> b) -> a -> b
$ Stream (Double, Double) -> (Double, Double)
forall a. Stream a -> a
Stream.head (Stream (Double, Double) -> (Double, Double))
-> Stream (Double, Double) -> (Double, Double)
forall a b. (a -> b) -> a -> b
$ ((Double, Double) -> Bool)
-> Stream (Double, Double) -> Stream (Double, Double)
forall a. (a -> Bool) -> Stream a -> Stream a
Stream.dropWhile (\(Double
v1, Double
v2) -> Double -> Double
forall a. Num a => a -> a
abs (Double
v1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
v2) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.01 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
arcsecond) (Stream (Double, Double) -> Stream (Double, Double))
-> Stream (Double, Double) -> Stream (Double, Double)
forall a b. (a -> b) -> a -> b
$ Stream Double -> Stream Double -> Stream (Double, Double)
forall a b. Stream a -> Stream b -> Stream (a, b)
Stream.zip Stream Double
lats (Stream Double -> Stream (Double, Double))
-> Stream Double -> Stream (Double, Double)
forall a b. (a -> b) -> a -> b
$ Int -> Stream Double -> Stream Double
forall a. Int -> Stream a -> Stream a
Stream.drop Int
1 Stream Double
lats
   gridEllipsoid :: GridStereo e -> e
gridEllipsoid = Geodetic e -> e
forall e. Geodetic e -> e
ellipsoid (Geodetic e -> e)
-> (GridStereo e -> Geodetic e) -> GridStereo e -> e
forall b c a. (b -> c) -> (a -> b) -> a -> c
. GridStereo e -> Geodetic e
forall e. GridStereo e -> Geodetic e
gridTangent