{-# 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 :: Cl3
e0 = Double -> Cl3
R Double
1
e1 :: Cl3
e1 = Double -> Double -> Double -> Cl3
V3 Double
1 Double
0 Double
0
e2 :: Cl3
e2 = Double -> Double -> Double -> Cl3
V3 Double
0 Double
1 Double
0
e3 :: Cl3
e3 = Double -> Double -> Double -> Cl3
V3 Double
0 Double
0 Double
1

i :: Cl3
i = Double -> Cl3
I Double
1

p1 :: Cl3
p1 = Cl3
0.5 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
e0 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
+ Cl3
e1)
p2 :: Cl3
p2 = Cl3
0.5 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
e0 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
+ Cl3
e2)
p3 :: Cl3
p3 = Cl3
0.5 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
e0 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
+ Cl3
e3)




-- | 'hpv' horizontally polarized Jones vector
hpv :: Cl3
hpv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3
e0 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
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 :: Cl3
vpv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
exp ((-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
forall a. Floating a => a
pi Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
e2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
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 :: Cl3
dpv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
exp ((-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
forall a. Floating a => a
piCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
e2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
p3  -- rotate -e1 to e3 around rotational axis e2, an angle of pi/2

-- | 'apv' anti-diagonally polarized Jones vector
apv :: Cl3
apv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
exp ((-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
forall a. Floating a => a
piCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (-Cl3
e2)) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
p3  -- rotate e1 to e3 around rotational axis -e2, an angle of pi/2

-- | 'rpv' right hand circularly polarized Jones vector
rpv :: Cl3
rpv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
exp ((-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
forall a. Floating a => a
piCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (-Cl3
e1)) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
p3  -- rotate -e2 to e3 around rotational axis -e1, and angle of pi/2

-- | 'lpv' left hand circularly polarized Jones vector
lpv :: Cl3
lpv = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
exp ((-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (Cl3
forall a. Floating a => a
piCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
e1) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
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 :: Cl3 -> Cl3
jv (Cl3 -> Cl3
forall a. Num a => a -> a
signum(Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Cl3 -> Cl3
toV3 -> Cl3
v) | Cl3
v Cl3 -> Cl3 -> Bool
forall a. Eq a => a -> a -> Bool
== Cl3
e3 = Cl3
hpv
                      | Cl3
v Cl3 -> Cl3 -> Bool
forall a. Eq a => a -> a -> Bool
== -Cl3
e3 = Cl3
vpv
                      | Bool
otherwise = Cl3 -> Cl3
forall a. Num a => a -> a
signum (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ Cl3 -> Cl3
forall a. Floating a => a -> a
sqrt (Cl3
e3 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
v) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
p3


-- | 'hpm' Horizontal Polarizer Jones Matrix
hpm :: Cl3
hpm = Cl3
p3

-- | 'vpm' Vertical Polarizer Jones Matrix
vpm :: Cl3
vpm = Cl3 -> Cl3
bar Cl3
p3

-- | 'dpm' Diagonal Polarizer Jones Matrix
dpm :: Cl3
dpm = Cl3
p1

-- | 'apm' Anti-diagonal Polarizer Jones Matrix
apm :: Cl3
apm = Cl3 -> Cl3
bar Cl3
p1

-- | 'rpm' Right Hand Circular Polarizer Jones Matrix
rpm :: Cl3
rpm = Cl3
p2

-- | 'lpm' Left Hand Circular Polarizer Jones Matrix
lpm :: Cl3
lpm = Cl3 -> Cl3
bar Cl3
p2


-- | 'jm' funciton that returns a Jones Matrix from an input Bloch Vector
-- currently converts the input to a unit vector
jm :: Cl3 -> Cl3
jm (Cl3 -> Cl3
forall a. Num a => a -> a
signum(Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Cl3 -> Cl3
toV3 -> Cl3
v) = Cl3 -> Cl3
project Cl3
v

-- | 'rot' will produce a versor that rotates a vector by an angle about
-- a unit vector axis.
rot :: Cl3 -> Cl3 -> Cl3
rot (Cl3 -> Cl3
toR -> Cl3
theta) (Cl3 -> Cl3
forall a. Num a => a -> a
signum(Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Cl3 -> Cl3
toV3 -> Cl3
axis) = Cl3 -> Cl3
forall a. Floating a => a -> a
exp (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ (-Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
theta Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
axis

-- | 'rotIsh' will produce a versor that rotates by double the input angle
-- for rotating polarizers and wave plates the axis is e2.
rotIsh :: Cl3 -> Cl3
rotIsh (Cl3 -> Cl3
toR -> Cl3
theta) = Cl3 -> Cl3 -> Cl3
rot (Cl3
2Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
*Cl3
theta) Cl3
e2

-- | 'hpmRot' Jones matrix for a rotated ideal Linear Horizontal Polarizer.
-- Input value should be a scalar angle in Radians.
hpmRot :: Cl3 -> Cl3
hpmRot (Cl3 -> Cl3
toR -> Cl3
theta) = 
  let roted :: Cl3
roted = Cl3 -> Cl3
rotIsh Cl3
theta
  in Cl3
roted Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
hpm Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
dag Cl3
roted

-- | 'qwp' Quarter Wave Plate Jones Matrix
qwp :: Cl3
qwp = Cl3
p3 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
- Cl3
i Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
bar Cl3
p3

-- | 'qwpRot' Rotated Quarter Wave Plate Jones Matrix.
-- Input value should be a scalar angle in Radians.
qwpRot :: Cl3 -> Cl3
qwpRot (Cl3 -> Cl3
toR -> Cl3
theta) = 
  let roted :: Cl3
roted = Cl3 -> Cl3
rotIsh Cl3
theta
  in Cl3
roted Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
qwp Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
dag Cl3
roted

-- | 'hwp' Half Wave Plate Jones Matrix
hwp :: Cl3
hwp = Cl3
e3

-- | 'hwpRot' Rotated Half Wave Plate Jones Matrix.
-- Input value should be a scalar angle in Radians.
hwpRot :: Cl3 -> Cl3
hwpRot (Cl3 -> Cl3
toR -> Cl3
theta) = 
  let roted :: Cl3
roted = Cl3 -> Cl3
rotIsh Cl3
theta
  in Cl3
roted Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
hwp Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
dag Cl3
roted

-- | 'wp' a Wave Plate with phase shift of phi Jones Matrix.
-- Input value should be a scalar angle in Radians.
wp :: Cl3 -> Cl3
wp (Cl3 -> Cl3
toR -> Cl3
phi) = Cl3 -> Cl3
forall a. Floating a => a -> a
exp (Cl3 -> Cl3) -> Cl3 -> Cl3
forall a b. (a -> b) -> a -> b
$ (Cl3
iCl3 -> Cl3 -> Cl3
forall a. Fractional a => a -> a -> a
/Cl3
2) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
phi Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
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 :: Cl3 -> Cl3 -> Cl3
wpRot (Cl3 -> Cl3
toR -> Cl3
phi) (Cl3 -> Cl3
toR -> Cl3
theta) = 
  let roted :: Cl3
roted = Cl3 -> Cl3
rotIsh Cl3
theta
  in Cl3
roted Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
wp Cl3
phi Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
dag Cl3
roted

-- | 'refl' a Refelection Jones Matrix
refl :: Cl3
refl = Cl3
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 :: Cl3 -> (Cl3, Cl3, Cl3)
factorize Cl3
jonesVec = 
  let c :: Cl3
c = Cl3 -> Cl3
toC Cl3
jonesVec
      jonesVec' :: Cl3
jonesVec' = Cl3 -> Cl3
forall a. Fractional a => a -> a
recip Cl3
c Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
jonesVec
      ampC :: Cl3
ampC = Cl3 -> Cl3
forall a. Num a => a -> a
abs Cl3
c
      ampJonesVec' :: Cl3
ampJonesVec' = Cl3 -> Cl3
forall a. Num a => a -> a
abs Cl3
jonesVec'
      normJonesVec :: Cl3
normJonesVec = Cl3 -> Cl3
forall a. Fractional a => a -> a
recip Cl3
ampJonesVec' Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
jonesVec'
      amp :: Cl3
amp = Cl3
ampC Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
ampJonesVec'
      normC :: Cl3
normC = Cl3 -> Cl3
forall a. Fractional a => a -> a
recip Cl3
ampC Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3
c
      phi :: Cl3
phi = Cl3
2 Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* (-Cl3
i) Cl3 -> Cl3 -> Cl3
forall a. Num a => a -> a -> a
* Cl3 -> Cl3
forall a. Floating a => a -> a
log Cl3
normC
  in (Cl3
amp, Cl3
phi, Cl3
normJonesVec)

#ifndef O_NO_RANDOM
-------------------------------------------------------------------
--
--  Random Jones Vectors
--
-------------------------------------------------------------------

-- | 'randJonesVec' a Random Jones Vector.
randJonesVec :: RandomGen g => g -> (Cl3, g)
randJonesVec :: g -> (Cl3, g)
randJonesVec g
g =
  let (Cl3
v3, g
g') = g -> (Cl3, g)
forall g. RandomGen g => g -> (Cl3, g)
randUnitV3 g
g
  in (Cl3 -> Cl3
jv Cl3
v3,g
g')

-- | 'randOrthogonalJonesVec' a Random Orthogonal Complementary pair of Jones
-- Vectors.
randOrthogonalJonesVec :: RandomGen g => g -> ((Cl3, Cl3), g)
randOrthogonalJonesVec :: g -> ((Cl3, Cl3), g)
randOrthogonalJonesVec g
g = 
  let (Cl3
v3, g
g') = g -> (Cl3, g)
forall g. RandomGen g => g -> (Cl3, g)
randUnitV3 g
g
  in ((Cl3 -> Cl3
jv Cl3
v3, Cl3 -> Cl3
jv (Cl3 -> Cl3
bar Cl3
v3)),g
g')

#endif