    {-# LANGUAGE CPP #-}
    ---------------------------------------------------------------------------
    -- |
    -- Copyright   :  (C) 2015 Edward Kmett
    -- License     :  BSD-style (see the file LICENSE)
    --
    -- Maintainer  :  Edward Kmett <ekmett@gmail.com>
    -- Stability   :  experimental
    -- Portability :  non-portable
    --
    -- Common projection matrices: e.g. perspective/orthographic transformation
    -- matrices.
    --
    -- Analytically derived inverses are also supplied, because they can be
    -- much more accurate in practice than computing them through general
    -- purpose means
    ---------------------------------------------------------------------------
    module Linear.Projection
      ( lookAt
      , perspective, inversePerspective
      , infinitePerspective, inverseInfinitePerspective
      , frustum, inverseFrustum
      , ortho, inverseOrtho
      ) where

    import Control.Lens hiding (index)
    import Linear.V3
    import Linear.V4
    import Linear.Matrix
    import Linear.Epsilon
    import Linear.Metric

    #ifdef HLINT
    {-# ANN module "HLint: ignore Reduce duplication" #-}
    #endif

    -- | Build a look at view matrix
    lookAt
      :: (Epsilon a, Floating a)
      => V3 a -- ^ Eye
      -> V3 a -- ^ Center
      -> V3 a -- ^ Up
      -> M44 a
    lookAt eye center up =
      V4 (V4 (xa^._x)  (xa^._y)  (xa^._z)  xd)
         (V4 (ya^._x)  (ya^._y)  (ya^._z)  yd)
         (V4 (-za^._x) (-za^._y) (-za^._z) zd)
         (V4 0         0         0          1)
      where za = normalize $ center - eye
            xa = normalize $ cross za up
            ya = cross xa za
            xd = -dot xa eye
            yd = -dot ya eye
            zd = dot za eye

    -- | Build a matrix for a symmetric perspective-view frustum
    perspective
      :: Floating a
      => a -- ^ FOV (y direction, in radians)
      -> a -- ^ Aspect ratio
      -> a -- ^ Near plane
      -> a -- ^ Far plane
      -> M44 a
    perspective fovy aspect near far =
      V4 (V4 x 0 0    0)
         (V4 0 y 0    0)
         (V4 0 0 z    w)
         (V4 0 0 (-1) 0)
      where tanHalfFovy = tan $ fovy / 2
            x = 1 / (aspect * tanHalfFovy)
            y = 1 / tanHalfFovy
            fpn = far + near
            fmn = far - near
            oon = 1/near
            oof = 1/far
            -- z = 1 / (near/fpn - far/fpn) -- would be better by .5 bits
            z = -fpn/fmn
            w = 1/(oof-oon) -- 13 bits error reduced to 0.17
            -- w = -(2 * far * near) / fmn

    #ifdef HERBIE
    {-# ANN perspective "NoHerbie" #-}
    #endif

    -- | Build an inverse perspective matrix
    inversePerspective
      :: Floating a
      => a -- ^ FOV (y direction, in radians)
      -> a -- ^ Aspect ratio
      -> a -- ^ Near plane
      -> a -- ^ Far plane
      -> M44 a
    inversePerspective fovy aspect near far =
      V4 (V4 a 0 0 0   )
         (V4 0 b 0 0   )
         (V4 0 0 0 (-1))
         (V4 0 0 c d   )
      where tanHalfFovy = tan $ fovy / 2
            a = aspect * tanHalfFovy
            b = tanHalfFovy
            c = oon - oof
            d = oon + oof
            oon = 1/near
            oof = 1/far

    -- | Build a perspective matrix per the classic @glFrustum@ arguments.
    frustum
      :: Floating a
      => a -- ^ Left
      -> a -- ^ Right
      -> a -- ^ Bottom
      -> a -- ^ Top
      -> a -- ^ Near
      -> a -- ^ Far
      -> M44 a
    frustum l r b t n f =
      V4 (V4 x 0 a    0)
         (V4 0 y e    0)
         (V4 0 0 c    d)
         (V4 0 0 (-1) 0)
      where
        rml = r-l
        tmb = t-b
        fmn = f-n
        x = 2*n/rml
        y = 2*n/tmb
        a = (r+l)/rml
        e = (t+b)/tmb
        c = negate (f+n)/fmn
        d = (-2*f*n)/fmn

    inverseFrustum
      :: Floating a
      => a -- ^ Left
      -> a -- ^ Right
      -> a -- ^ Bottom
      -> a -- ^ Top
      -> a -- ^ Near
      -> a -- ^ Far
      -> M44 a
    inverseFrustum l r b t n f =
      V4 (V4 rx 0 0 ax)
         (V4 0 ry 0 by)
         (V4 0 0 0 (-1))
         (V4 0 0 rd cd)
      where
        hrn  = 0.5/n
        hrnf = 0.5/(n*f)
        rx = (r-l)*hrn
        ry = (t-b)*hrn
        ax = (r+l)*hrn
        by = (t+b)*hrn
        cd = (f+n)*hrnf
        rd = (n-f)*hrnf

    -- | Build a matrix for a symmetric perspective-view frustum with a far plane at infinite
    infinitePerspective
      :: Floating a
      => a -- ^ FOV (y direction, in radians)
      -> a -- ^ Aspect Ratio
      -> a -- ^ Near plane
      -> M44 a
    infinitePerspective fovy a n =
      V4 (V4 x 0 0    0)
         (V4 0 y 0    0)
         (V4 0 0 (-1) w)
         (V4 0 0 (-1) 0)
      where
        t = n*tan(fovy/2)
        b = -t
        l = b*a
        r = t*a
        x = (2*n)/(r-l)
        y = (2*n)/(t-b)
        w = -2*n

    inverseInfinitePerspective
      :: Floating a
      => a -- ^ FOV (y direction, in radians)
      -> a -- ^ Aspect Ratio
      -> a -- ^ Near plane
      -> M44 a
    inverseInfinitePerspective fovy a n =
      V4 (V4 rx 0 0  0)
         (V4 0 ry 0  0)
         (V4 0 0  0  (-1))
         (V4 0 0  rw (-rw))
      where
        t = n*tan(fovy/2)
        b = -t
        l = b*a
        r = t*a
        hrn = 0.5/n
        rx = (r-l)*hrn
        ry = (t-b)*hrn
        rw = -hrn

    -- | Build an orthographic perspective matrix from 6 clipping planes.
    -- This matrix takes the region delimited by these planes and maps it
    -- to normalized device coordinates between [-1,1]
    --
    -- This call is designed to mimic the parameters to the OpenGL @glOrtho@
    -- call, so it has a slightly strange convention: Notably: the near and
    -- far planes are negated.
    --
    -- Consequently:
    --
    -- @
    -- 'ortho' l r b t n f !* 'V4' l b (-n) 1 = 'V4' (-1) (-1) (-1) 1
    -- 'ortho' l r b t n f !* 'V4' r t (-f) 1 = 'V4' 1 1 1 1
    -- @
    --
    -- Examples:
    --
    -- >>> ortho 1 2 3 4 5 6 !* V4 1 3 (-5) 1
    -- V4 (-1.0) (-1.0) (-1.0) 1.0
    --
    -- >>> ortho 1 2 3 4 5 6 !* V4 2 4 (-6) 1
    -- V4 1.0 1.0 1.0 1.0
    ortho
      :: Fractional a
      => a -- ^ Left
      -> a -- ^ Right
      -> a -- ^ Bottom
      -> a -- ^ Top
      -> a -- ^ Near
      -> a -- ^ Far
      -> M44 a
    ortho l r b t n f =
      V4 (V4 (-2*x) 0      0     ((r+l)*x))
         (V4 0      (-2*y) 0     ((t+b)*y))
         (V4 0      0      (2*z) ((f+n)*z))
         (V4 0      0      0     1)
      where x = recip(l-r)
            y = recip(b-t)
            z = recip(n-f)

    -- | Build an inverse orthographic perspective matrix from 6 clipping planes
    inverseOrtho
      :: Fractional a
      => a -- ^ Left
      -> a -- ^ Right
      -> a -- ^ Bottom
      -> a -- ^ Top
      -> a -- ^ Near
      -> a -- ^ Far
      -> M44 a
    inverseOrtho l r b t n f =
      V4 (V4 x 0 0 c)
         (V4 0 y 0 d)
         (V4 0 0 z e)
         (V4 0 0 0 1)
      where x = 0.5*(r-l)
            y = 0.5*(t-b)
            z = 0.5*(n-f)
            c = 0.5*(l+r)
            d = 0.5*(b+t)
            e = -0.5*(n+f)