content: {-# LANGUAGE ViewPatterns #-}
  {-# OPTIONS_GHC -fno-warn-missing-signatures #-}
  {-# LANGUAGE CPP #-}
  
  #ifndef O_LIQUID
  {-# LANGUAGE Safe #-}
  #endif
  
  #if __GLASGOW_HASKELL__ == 810
  -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
  -- This code is meant to reproduce MR 2608 for GHC 8.10
  {-# OPTIONS_GHC -funfolding-keeness-factor=1 -funfolding-use-threshold=80 #-}
  #endif
  
  --------------------------------------------------------------------------------------------
  -- |
  -- Copyright   :  (C) 2018 Nathan Waivio
  -- License     :  BSD3
  -- Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>
  -- Stability   :  Stable
  -- Portability :  unportable
  --
  -- Library implementing standard functions for the Jones Calculus in the Cl3 Library.
  -- This implementation of the Jones Calculus is based on the convensions of SPIE's Field Guide to Polarization (ϕ = ω t − k z).
  -- 
  -- * E. Collett, Field Guide to Polarization, SPIE Field Guides vol. FG05, SPIE (2005). ISBN 0-8194-5868-6.
  -- 
  --  
  -- = Jones Vectors
  -- 
  -- Within the system of the Bloch Sphere, the Jones Vectors in Cl3 are calculated
  -- by generating the left ideal of the rotation of a unit vector to the 'e3' basis.
  -- Standard form for for a versor is 'rot = exp $ (-i/2) * theta * u' for angle 'theta'
  -- and the rotational axis unit vector 'u'.
  --
  --          Bloch Sphere Coordinates:
  -- 
  -- @
  --                 e3
  --                 |
  --                 |____e2
  --                / 
  --               /
  --              e1
  -- @
  --
  --------------------------------------------
  
  module Algebra.Geometric.Cl3.JonesCalculus
  (
   -- * Jones Vectors
   hpv, vpv,
   dpv, apv,
   rpv, lpv,
   jv,
   -- * Jones Matrices
   hpm, vpm,
   dpm, apm,
   rpm, lpm,
   jm,
   hpmRot,
   -- * Wave Plates
   qwp, hwp,
   qwpRot, hwpRot,
   wp,
   wpRot,
   -- * Reflection
   refl,
  #ifndef O_NO_RANDOM
   -- * Random Jones Vectors
   randJonesVec,
   randOrthogonalJonesVec,
  #endif
   -- * Normalization Factorization
   factorize
  ) where
  
  
  import safe Algebra.Geometric.Cl3 (Cl3(..), dag, bar, toR, toV3, toC, project)
  
  #ifndef O_NO_RANDOM
  import safe Algebra.Geometric.Cl3 (randUnitV3)
  import System.Random (RandomGen)
  #endif
  
  e0 = R 1
  e1 = V3 1 0 0
  e2 = V3 0 1 0
  e3 = V3 0 0 1
  
  i = I 1
  
  p1 = 0.5 * (e0 + e1)
  p2 = 0.5 * (e0 + e2)
  p3 = 0.5 * (e0 + e3)
  
  
  
  
  -- | 'hpv' horizontally polarized Jones vector
  hpv = signum $ e0 * p3  -- e0 == exp $ (-i/2) * 0 * e2, any vector orthoganl to e3 could have been selected as the rotational axis because the angle is zero
  
  -- | 'vpv' vertically polarized Jones vector
  vpv = signum $ exp ((-i/2) * pi * e2) * p3  -- e2 is selected to obtain the standard form, e1 or any vector orthoganl to e3 could have been selected
  
  -- | 'dpv' diagonally polarized Jones vector
  dpv = signum $ exp ((-i/2) * (pi/2) * e2) * p3  -- rotate -e1 to e3 around rotational axis e2, an angle of pi/2
  
  -- | 'apv' anti-diagonally polarized Jones vector
  apv = signum $ exp ((-i/2) * (pi/2) * (-e2)) * p3  -- rotate e1 to e3 around rotational axis -e2, an angle of pi/2
  
  -- | 'rpv' right hand circularly polarized Jones vector
  rpv = signum $ exp ((-i/2) * (pi/2) * (-e1)) * p3  -- rotate -e2 to e3 around rotational axis -e1, and angle of pi/2
  
  -- | 'lpv' left hand circularly polarized Jones vector
  lpv = signum $ exp ((-i/2) * (pi/2) * e1) * p3  -- rotate e2 to e3 around rotational axis e1, an angle of pi/2
  
  -- | 'jv' function that returns Jones vector from input vector unit vector
  -- currently converts the input to a unit vector
  jv (signum.toV3 -> v) | v == e3 = hpv
                        | v == -e3 = vpv
                        | otherwise = signum $ sqrt (e3 * v) * p3
  
  
  -- | 'hpm' Horizontal Polarizer Jones Matrix
  hpm = p3
  
  -- | 'vpm' Vertical Polarizer Jones Matrix
  vpm = bar p3
  
  -- | 'dpm' Diagonal Polarizer Jones Matrix
  dpm = p1
  
  -- | 'apm' Anti-diagonal Polarizer Jones Matrix
  apm = bar p1
  
  -- | 'rpm' Right Hand Circular Polarizer Jones Matrix
  rpm = p2
  
  -- | 'lpm' Left Hand Circular Polarizer Jones Matrix
  lpm = bar p2
  
  
  -- | 'jm' funciton that returns a Jones Matrix from an input Bloch Vector
  -- currently converts the input to a unit vector
  jm (signum.toV3 -> v) = project v
  
  -- | 'rot' will produce a versor that rotates a vector by an angle about
  -- a unit vector axis.
  rot (toR -> theta) (signum.toV3 -> axis) = exp $ (-i/2) * theta * axis
  
  -- | 'rotIsh' will produce a versor that rotates by double the input angle
  -- for rotating polarizers and wave plates the axis is e2.
  rotIsh (toR -> theta) = rot (2*theta) e2
  
  -- | 'hpmRot' Jones matrix for a rotated ideal Linear Horizontal Polarizer.
  -- Input value should be a scalar angle in Radians.
  hpmRot (toR -> theta) = 
    let roted = rotIsh theta
    in roted * hpm * dag roted
  
  -- | 'qwp' Quarter Wave Plate Jones Matrix
  qwp = p3 - i * bar p3
  
  -- | 'qwpRot' Rotated Quarter Wave Plate Jones Matrix.
  -- Input value should be a scalar angle in Radians.
  qwpRot (toR -> theta) = 
    let roted = rotIsh theta
    in roted * qwp * dag roted
  
  -- | 'hwp' Half Wave Plate Jones Matrix
  hwp = e3
  
  -- | 'hwpRot' Rotated Half Wave Plate Jones Matrix.
  -- Input value should be a scalar angle in Radians.
  hwpRot (toR -> theta) = 
    let roted = rotIsh theta
    in roted * hwp * dag roted
  
  -- | 'wp' a Wave Plate with phase shift of phi Jones Matrix.
  -- Input value should be a scalar angle in Radians.
  wp (toR -> phi) = exp $ (i/2) * phi * e3
  
  -- | 'wpRot' a Rotated Wave Plate with phase shift of phi and rotation theta Jones Matrix.
  -- The first input value is phi the phase shift as a scalar value in Radians. The second
  -- input value is theta the rotation a scalar angle in Radians.
  wpRot (toR -> phi) (toR -> theta) = 
    let roted = rotIsh theta
    in roted * wp phi * dag roted
  
  -- | 'refl' a Refelection Jones Matrix
  refl = e3
  
  
  -- | 'factorize' is a function that takes an Jones Vector after transformation by an 
  -- optical chain, and returns the amplitude (amp), phase (phi), and normalized Jones 
  -- Vector (vec), by the factorization of the input such that: @__amp * exp (i*phi/2) * vec__@
  factorize :: Cl3 -> (Cl3,Cl3,Cl3)
  factorize jonesVec = 
    let c = toC jonesVec
        jonesVec' = recip c * jonesVec
        ampC = abs c
        ampJonesVec' = abs jonesVec'
        normJonesVec = recip ampJonesVec' * jonesVec'
        amp = ampC * ampJonesVec'
        normC = recip ampC * c
        phi = 2 * (-i) * log normC
    in (amp, phi, normJonesVec)
  
  #ifndef O_NO_RANDOM
  -------------------------------------------------------------------
  --
  --  Random Jones Vectors
  --
  -------------------------------------------------------------------
  
  -- | 'randJonesVec' a Random Jones Vector.
  randJonesVec :: RandomGen g => g -> (Cl3, g)
  randJonesVec g =
    let (v3, g') = randUnitV3 g
    in (jv v3,g')
  
  -- | 'randOrthogonalJonesVec' a Random Orthogonal Complementary pair of Jones
  -- Vectors.
  randOrthogonalJonesVec :: RandomGen g => g -> ((Cl3, Cl3), g)
  randOrthogonalJonesVec g = 
    let (v3, g') = randUnitV3 g
    in ((jv v3, jv (bar v3)),g')
  
  #endif
  