-- | Module containing several useful projections for generating maps.
module GIS.Math.Projections where

import           Control.Lens
import           GIS.Graphics.Types
import           GIS.Math.Utils
import           GIS.Types

-- | For use as a reference point in certain projections.
washingtonDC :: Point
washingtonDC = over _2 radians $ over _1 radians (38.9072, -77.0369)

-- | For use as a reference point with the Craig retroazimuthal projection
mecca :: Point
mecca = over _2 radians $ over _1 radians (21.3891, 39.8579)

-- | Littow retroazimuthal + conformal projection
littow :: Projection
littow (long, lat) = (sin(long - referenceLong)/cos lat, cos (long - referenceLong) * tan lat)
    where referenceLong = radians 0

-- | Craig retroazimuthal projection
-- (works on a subset of the world)
craig :: Point -> Projection
craig referencePoint (long, lat) = (long - referenceLong, y)
    where (referenceLong, referenceLat) = referencePoint
          y = if long - referenceLong == 0 then expr else (long - referenceLong)/sin(long - referenceLong) * expr
          expr = sin lat * cos (long - referenceLong) - tan referenceLat * cos lat

-- | Winkel Tripel projection (standard for the National Geographic Society since
winkel3 :: Projection
winkel3 (long, lat) = ((lambda * cos phi1 + (2 * cos lat * sin (lambda/2)/sinc alpha))/2, (lat + sin lat/sinc alpha)/2)
    where lambda = long - lambda0
          phi1 = acos $ 2 / pi
          alpha = acos $ cos lat * cos (lambda/2)
          lambda0 = radians (-77.0369)

-- | Mercator projection.
mercator :: Projection
mercator (long, lat) = (long - meridian, asinh (tan lat))
    where meridian = radians (-98.5795)

-- | Bonne projection with standard parallel at 45 N and central meridian
-- centered at Washington DC
bonne :: Projection
bonne (long, lat) = (rho * sin e, cot phi1 - rho * cos e)
    where rho = cot phi1 + phi1 - lat
          e = (long - meridian) * cos lat / rho
          phi1 = radians 45 -- standard parallel @ 45 N
          meridian = radians (-77.0369) -- central meridian @ dc
          cot = (1/) . tan

-- | Albers projection for a given reference point. To make it usable you can
-- use
-- > ablers washingtonDC
albers :: Point -> Projection
albers referencePoint (long, lat) = (rho * sin theta, rho' - rho * cos theta)
    where n = (sin phi1 + sin phi2)/2
          theta = n * (long - referenceLong)
          c = cos phi1^(2 :: Int) + 2 * n * sin phi1
          rho = sqrt (c - 2 * n * sin lat) / n
          rho' = sqrt (c - 2 * n * sin referenceLat) / n
          phi1 = radians 20 -- standard parallels @ 20, 50 degrees
          phi2 = radians 50
          (referenceLong, referenceLat) = referencePoint

-- | Helper to project given a `Polygon`.
project :: Projection -> Polygon -> Polygon
project f = fmap (f . toRadians)

-- | Helper to apply a projection given a `Map`.
projectMap :: Projection -> Map -> Map
projectMap p = over labelledDistricts (fmap (over _1  (fmap (project p))))