{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Safe #-}
{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
{-# LANGUAGE CPP #-}

{- | 
Module      :  Physics.Learn.Ket
Copyright   :  (c) Scott N. Walck 2016-2018
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <[email protected]>
Stability   :  experimental

This module contains ket vectors, bra vectors,
and operators for quantum mechanics.
-}

-- a Ket layer on top of QuantumMat

module Physics.Learn.Ket
    (
    -- * Basic data types
      C
    , i
    , magnitude
    , Ket
    , Bra
    , Operator
    -- * Kets for spin-1/2 particles
    , xp
    , xm
    , yp
    , ym
    , zp
    , zm
    , np
    , nm
    -- * Operators for spin-1/2 particles
    , sx
    , sy
    , sz
    , sn
    , sn'
    -- * Quantum Dynamics
    , timeEvOp
    , timeEv
    -- * Composition
    , Kron(..)
    -- * Measurement
    , possibleOutcomes
    , outcomesProjectors
    , outcomesProbabilities
--    , prob
--    , probs
    -- * Generic multiplication
    , Mult(..)
    -- * Adjoint operation
    , Dagger(..)
    -- * Normalization
    , HasNorm(..)
    -- * Representation
    , Representable(..)
    -- * Orthonormal bases
    , OrthonormalBasis
    , makeOB
    , listBasis
    , size
    -- * Orthonormal bases for spin-1/2 particles
    , xBasis
    , yBasis
    , zBasis
    , nBasis
    -- , angularMomentumXMatrix
    -- , angularMomentumYMatrix
    -- , angularMomentumZMatrix
    -- , angularMomentumPlusMatrix
    -- , angularMomentumMinusMatrix
    -- , jXMatrix
    -- , jYMatrix
    -- , jZMatrix
    -- , matrixCommutator
    -- , rotationMatrix
    -- , jmColumn
    )
    where

-- We try to import only from QuantumMat
-- and not from Numeric.LinearAlgebra

import qualified Data.Complex as C
import Data.Complex
    ( Complex(..)
    , conjugate
    )
import qualified Physics.Learn.QuantumMat as M
import Physics.Learn.QuantumMat
    ( C
    , Vector
    , Matrix
    , (#>)
    , (<#)
    , conjugateTranspose
    , scaleV
    , scaleM
    , conjV
    , fromList
    , toList
    , fromLists
    )
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

infixl 7 <>

-- | A ket vector describes the state of a quantum system.
data Ket = Ket (Vector C)

instance Show Ket where
    show :: Ket -> String
show Ket
k =
        let message :: String
message = String
"Use 'rep <basis name> <ket name>'."
        in if Ket -> Int
forall a b. Representable a b => a -> Int
dim Ket
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
           then String
"Representation in zBasis:\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++
                Vector (Complex Double) -> String
forall a. Show a => a -> String
show (OrthonormalBasis -> Ket -> Vector (Complex Double)
forall a b. Representable a b => OrthonormalBasis -> a -> b
rep OrthonormalBasis
zBasis Ket
k) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
message
           else String
message

-- | An operator describes an observable (a Hermitian operator)
--   or an action (a unitary operator).
data Operator = Operator (Matrix C)

instance Show Operator where
    show :: Operator -> String
show Operator
op =
        let message :: String
message = String
"Use 'rep <basis name> <operator name>'."
        in if Operator -> Int
forall a b. Representable a b => a -> Int
dim Operator
op Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
           then String
"Representation in zBasis:\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++
                Matrix (Complex Double) -> String
forall a. Show a => a -> String
show (OrthonormalBasis -> Operator -> Matrix (Complex Double)
forall a b. Representable a b => OrthonormalBasis -> a -> b
rep OrthonormalBasis
zBasis Operator
op) String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
message
           else String
message

-- | A bra vector describes the state of a quantum system.
data Bra = Bra (Vector C)

instance Show Bra where
    show :: Bra -> String
show Bra
_ = String
"<bra>\nTry 'rep zBasis <bra name>'"

magnitude :: C -> Double
magnitude :: Complex Double -> Double
magnitude = Complex Double -> Double
forall a. RealFloat a => Complex a -> a
C.magnitude

i :: C
i :: Complex Double
i = Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
1

-- | Generic multiplication including inner product,
--   outer product, operator product, and whatever else makes sense.
--   No conjugation takes place in this operation.
class Mult a b c | a b -> c where
    (<>) :: a -> b -> c

instance Mult C C C where
    Complex Double
z1 <> :: Complex Double -> Complex Double -> Complex Double
<> Complex Double
z2 = Complex Double
z1 Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
z2

instance Mult C Ket Ket where
    Complex Double
c <> :: Complex Double -> Ket -> Ket
<> Ket Vector (Complex Double)
matrixKet = Vector (Complex Double) -> Ket
Ket (Complex Double
-> Vector (Complex Double) -> Vector (Complex Double)
scaleV Complex Double
c Vector (Complex Double)
matrixKet)

instance Mult C Bra Bra where
    Complex Double
c <> :: Complex Double -> Bra -> Bra
<> Bra Vector (Complex Double)
matrixBra = Vector (Complex Double) -> Bra
Bra (Complex Double
-> Vector (Complex Double) -> Vector (Complex Double)
scaleV Complex Double
c Vector (Complex Double)
matrixBra)

instance Mult C Operator Operator where
    Complex Double
c <> :: Complex Double -> Operator -> Operator
<> Operator Matrix (Complex Double)
m = Matrix (Complex Double) -> Operator
Operator (Complex Double
-> Matrix (Complex Double) -> Matrix (Complex Double)
scaleM Complex Double
c Matrix (Complex Double)
m)

instance Mult Ket C Ket where
    Ket Vector (Complex Double)
matrixKet <> :: Ket -> Complex Double -> Ket
<> Complex Double
c = Vector (Complex Double) -> Ket
Ket (Complex Double
-> Vector (Complex Double) -> Vector (Complex Double)
scaleV Complex Double
c Vector (Complex Double)
matrixKet)

instance Mult Bra C Bra where
    Bra Vector (Complex Double)
matrixBra <> :: Bra -> Complex Double -> Bra
<> Complex Double
c = Vector (Complex Double) -> Bra
Bra (Complex Double
-> Vector (Complex Double) -> Vector (Complex Double)
scaleV Complex Double
c Vector (Complex Double)
matrixBra)

instance Mult Operator C Operator where
    Operator Matrix (Complex Double)
m <> :: Operator -> Complex Double -> Operator
<> Complex Double
c = Matrix (Complex Double) -> Operator
Operator (Complex Double
-> Matrix (Complex Double) -> Matrix (Complex Double)
scaleM Complex Double
c Matrix (Complex Double)
m)

instance Mult Bra Ket C where
    Bra Vector (Complex Double)
matrixBra <> :: Bra -> Ket -> Complex Double
<> Ket Vector (Complex Double)
matrixKet
        = [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Complex Double] -> Complex Double)
-> [Complex Double] -> Complex Double
forall a b. (a -> b) -> a -> b
$ (Complex Double -> Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double] -> [Complex Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
(*) (Vector (Complex Double) -> [Complex Double]
toList Vector (Complex Double)
matrixBra) (Vector (Complex Double) -> [Complex Double]
toList Vector (Complex Double)
matrixKet)

instance Mult Bra Operator Bra where
    Bra Vector (Complex Double)
matrixBra <> :: Bra -> Operator -> Bra
<> Operator Matrix (Complex Double)
matrixOp
        = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double)
matrixBra Vector (Complex Double)
-> Matrix (Complex Double) -> Vector (Complex Double)
<# Matrix (Complex Double)
matrixOp)

instance Mult Operator Ket Ket where
    Operator Matrix (Complex Double)
matrixOp <> :: Operator -> Ket -> Ket
<> Ket Vector (Complex Double)
matrixKet
        = Vector (Complex Double) -> Ket
Ket (Matrix (Complex Double)
matrixOp Matrix (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
#> Vector (Complex Double)
matrixKet)

instance Mult Ket Bra Operator where
    Ket Vector (Complex Double)
k <> :: Ket -> Bra -> Operator
<> Bra Vector (Complex Double)
b =
        Matrix (Complex Double) -> Operator
Operator
        ([[Complex Double]] -> Matrix (Complex Double)
fromLists [[ Complex Double
xComplex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
*Complex Double
y | Complex Double
y <- Vector (Complex Double) -> [Complex Double]
toList Vector (Complex Double)
b] | Complex Double
x <- Vector (Complex Double) -> [Complex Double]
toList Vector (Complex Double)
k])

instance Mult Operator Operator Operator where
    Operator Matrix (Complex Double)
m1 <> :: Operator -> Operator -> Operator
<> Operator Matrix (Complex Double)
m2 = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double)
m1 Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
M.<> Matrix (Complex Double)
m2)

instance Num Ket where
    Ket Vector (Complex Double)
v1 + :: Ket -> Ket -> Ket
+ Ket Vector (Complex Double)
v2 = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double)
v1 Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a -> a
+ Vector (Complex Double)
v2)
    Ket Vector (Complex Double)
v1 - :: Ket -> Ket -> Ket
- Ket Vector (Complex Double)
v2 = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double)
v1 Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a -> a
- Vector (Complex Double)
v2)
    * :: Ket -> Ket -> Ket
(*)         = String -> Ket -> Ket -> Ket
forall a. HasCallStack => String -> a
error String
"Multiplication is not defined for kets"
    negate :: Ket -> Ket
negate (Ket Vector (Complex Double)
v) = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a
negate Vector (Complex Double)
v)
    abs :: Ket -> Ket
abs         = String -> Ket -> Ket
forall a. HasCallStack => String -> a
error String
"abs is not defined for kets"
    signum :: Ket -> Ket
signum      = String -> Ket -> Ket
forall a. HasCallStack => String -> a
error String
"signum is not defined for kets"
    fromInteger :: Integer -> Ket
fromInteger = String -> Integer -> Ket
forall a. HasCallStack => String -> a
error String
"fromInteger is not defined for kets"

instance Num Bra where
    Bra Vector (Complex Double)
v1 + :: Bra -> Bra -> Bra
+ Bra Vector (Complex Double)
v2 = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double)
v1 Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a -> a
+ Vector (Complex Double)
v2)
    Bra Vector (Complex Double)
v1 - :: Bra -> Bra -> Bra
- Bra Vector (Complex Double)
v2 = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double)
v1 Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a -> a
- Vector (Complex Double)
v2)
    * :: Bra -> Bra -> Bra
(*)         = String -> Bra -> Bra -> Bra
forall a. HasCallStack => String -> a
error String
"Multiplication is not defined for bra vectors"
    negate :: Bra -> Bra
negate (Bra Vector (Complex Double)
v) = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double) -> Vector (Complex Double)
forall a. Num a => a -> a
negate Vector (Complex Double)
v)
    abs :: Bra -> Bra
abs         = String -> Bra -> Bra
forall a. HasCallStack => String -> a
error String
"abs is not defined for bra vectors"
    signum :: Bra -> Bra
signum      = String -> Bra -> Bra
forall a. HasCallStack => String -> a
error String
"signum is not defined for bra vectors"
    fromInteger :: Integer -> Bra
fromInteger = String -> Integer -> Bra
forall a. HasCallStack => String -> a
error String
"fromInteger is not defined for bra vectors"

instance Num Operator where
    Operator Matrix (Complex Double)
v1 + :: Operator -> Operator -> Operator
+ Operator Matrix (Complex Double)
v2 = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double)
v1 Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
forall a. Num a => a -> a -> a
+ Matrix (Complex Double)
v2)
    Operator Matrix (Complex Double)
v1 - :: Operator -> Operator -> Operator
- Operator Matrix (Complex Double)
v2 = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double)
v1 Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
forall a. Num a => a -> a -> a
- Matrix (Complex Double)
v2)
    Operator Matrix (Complex Double)
v1 * :: Operator -> Operator -> Operator
* Operator Matrix (Complex Double)
v2 = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double)
v1 Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
M.<> Matrix (Complex Double)
v2)
    negate :: Operator -> Operator
negate (Operator Matrix (Complex Double)
v) = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double) -> Matrix (Complex Double)
forall a. Num a => a -> a
negate Matrix (Complex Double)
v)
    abs :: Operator -> Operator
abs         = String -> Operator -> Operator
forall a. HasCallStack => String -> a
error String
"abs is not defined for operators"
    signum :: Operator -> Operator
signum      = String -> Operator -> Operator
forall a. HasCallStack => String -> a
error String
"signum is not defined for operators"
    fromInteger :: Integer -> Operator
fromInteger = String -> Integer -> Operator
forall a. HasCallStack => String -> a
error String
"fromInteger is not defined for operators"

-- | The adjoint operation on complex numbers, kets,
--   bras, and operators.
class Dagger a b | a -> b where
    dagger :: a -> b

instance Dagger Ket Bra where
    dagger :: Ket -> Bra
dagger (Ket Vector (Complex Double)
v) = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double) -> Vector (Complex Double)
conjV Vector (Complex Double)
v)

instance Dagger Bra Ket where
    dagger :: Bra -> Ket
dagger (Bra Vector (Complex Double)
v) = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double) -> Vector (Complex Double)
conjV Vector (Complex Double)
v)

instance Dagger Operator Operator where
    dagger :: Operator -> Operator
dagger (Operator Matrix (Complex Double)
m) = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double) -> Matrix (Complex Double)
conjugateTranspose Matrix (Complex Double)
m)

instance Dagger C C where
    dagger :: Complex Double -> Complex Double
dagger Complex Double
c = Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate Complex Double
c

class HasNorm a where
    norm      :: a -> Double
    normalize :: a -> a

instance HasNorm Ket where
    norm :: Ket -> Double
norm (Ket Vector (Complex Double)
v) = Vector (Complex Double) -> Double
M.norm Vector (Complex Double)
v
    normalize :: Ket -> Ket
normalize Ket
k  = (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Ket -> Double
forall a. HasNorm a => a -> Double
norm Ket
k Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
k

instance HasNorm Bra where
    norm :: Bra -> Double
norm (Bra Vector (Complex Double)
v) = Vector (Complex Double) -> Double
M.norm Vector (Complex Double)
v
    normalize :: Bra -> Bra
normalize Bra
b  = (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Bra -> Double
forall a. HasNorm a => a -> Double
norm Bra
b Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Bra -> Bra
forall a b c. Mult a b c => a -> b -> c
<> Bra
b

-- | An orthonormal basis of kets.
newtype OrthonormalBasis = OB [Ket]
    deriving (Int -> OrthonormalBasis -> ShowS
[OrthonormalBasis] -> ShowS
OrthonormalBasis -> String
(Int -> OrthonormalBasis -> ShowS)
-> (OrthonormalBasis -> String)
-> ([OrthonormalBasis] -> ShowS)
-> Show OrthonormalBasis
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> OrthonormalBasis -> ShowS
showsPrec :: Int -> OrthonormalBasis -> ShowS
$cshow :: OrthonormalBasis -> String
show :: OrthonormalBasis -> String
$cshowList :: [OrthonormalBasis] -> ShowS
showList :: [OrthonormalBasis] -> ShowS
Show)

-- | Make an orthonormal basis from a list of linearly independent kets.
makeOB :: [Ket] -> OrthonormalBasis
makeOB :: [Ket] -> OrthonormalBasis
makeOB = [Ket] -> OrthonormalBasis
OB ([Ket] -> OrthonormalBasis)
-> ([Ket] -> [Ket]) -> [Ket] -> OrthonormalBasis
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Ket] -> [Ket]
gramSchmidt

size :: OrthonormalBasis -> Int
size :: OrthonormalBasis -> Int
size (OB [Ket]
ks) = [Ket] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Ket]
ks

listBasis :: OrthonormalBasis -> [Ket]
listBasis :: OrthonormalBasis -> [Ket]
listBasis (OB [Ket]
ks) = [Ket]
ks

{-
newOrthonormalBasis :: Int -> OrthonormalBasis
newOrthonormalBasis = undefined
-}

class Representable a b | a -> b where
    rep :: OrthonormalBasis -> a -> b
    dim :: a -> Int

instance Representable Ket (Vector C) where
    rep :: OrthonormalBasis -> Ket -> Vector (Complex Double)
rep (OB [Ket]
ks) Ket
k = [Complex Double] -> Vector (Complex Double)
fromList ((Ket -> Complex Double) -> [Ket] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Ket
bk -> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
bk Bra -> Ket -> Complex Double
forall a b c. Mult a b c => a -> b -> c
<> Ket
k) [Ket]
ks)
    dim :: Ket -> Int
dim (Ket Vector (Complex Double)
v) = Vector (Complex Double) -> Int
M.dim Vector (Complex Double)
v

instance Representable Bra (Vector C) where
    rep :: OrthonormalBasis -> Bra -> Vector (Complex Double)
rep (OB [Ket]
ks) Bra
b = [Complex Double] -> Vector (Complex Double)
fromList ((Ket -> Complex Double) -> [Ket] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Ket
bk -> Bra
b Bra -> Ket -> Complex Double
forall a b c. Mult a b c => a -> b -> c
<> Ket
bk) [Ket]
ks)
    dim :: Bra -> Int
dim (Bra Vector (Complex Double)
v) = Vector (Complex Double) -> Int
M.dim Vector (Complex Double)
v

instance Representable Operator (Matrix C) where
    rep :: OrthonormalBasis -> Operator -> Matrix (Complex Double)
rep (OB [Ket]
ks) Operator
op = [[Complex Double]] -> Matrix (Complex Double)
fromLists [[ Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
k1 Bra -> Operator -> Bra
forall a b c. Mult a b c => a -> b -> c
<> Operator
op Bra -> Ket -> Complex Double
forall a b c. Mult a b c => a -> b -> c
<> Ket
k2 | Ket
k2 <- [Ket]
ks ] | Ket
k1 <- [Ket]
ks ]
    dim :: Operator -> Int
dim (Operator Matrix (Complex Double)
m) = let (Int
p,Int
q) = Matrix (Complex Double) -> (Int, Int)
M.size Matrix (Complex Double)
m
                       in if Int
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
q then Int
p else String -> Int
forall a. HasCallStack => String -> a
error String
"dim: non-square operator"

--------------
-- Spin 1/2 --
--------------

-- | State of a spin-1/2 particle if measurement
--   in the x-direction would give angular momentum +hbar/2.
xp :: Ket
xp :: Ket
xp = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.xp

-- | State of a spin-1/2 particle if measurement
--   in the x-direction would give angular momentum -hbar/2.
xm :: Ket
xm :: Ket
xm = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.xm

-- | State of a spin-1/2 particle if measurement
--   in the y-direction would give angular momentum +hbar/2.
yp :: Ket
yp :: Ket
yp = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.yp

-- | State of a spin-1/2 particle if measurement
--   in the y-direction would give angular momentum -hbar/2.
ym :: Ket
ym :: Ket
ym = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.ym

-- | State of a spin-1/2 particle if measurement
--   in the z-direction would give angular momentum +hbar/2.
zp :: Ket
zp :: Ket
zp = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.zp

-- | State of a spin-1/2 particle if measurement
--   in the z-direction would give angular momentum -hbar/2.
zm :: Ket
zm :: Ket
zm = Vector (Complex Double) -> Ket
Ket Vector (Complex Double)
M.zm

-- | State of a spin-1/2 particle if measurement
--   in the n-direction, described by spherical polar angle theta
--   and azimuthal angle phi, would give angular momentum +hbar/2.
np :: Double -> Double -> Ket
np :: Double -> Double -> Ket
np Double
theta Double
phi
    = (Double -> Double
forall a. Floating a => a -> a
cos (Double
theta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
zp
      Ket -> Ket -> Ket
forall a. Num a => a -> a -> a
+ (Double -> Double
forall a. Floating a => a -> a
sin (Double
theta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double -> Double
forall a. Floating a => a -> a
sin Double
phi) Complex Double -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
zm

-- | State of a spin-1/2 particle if measurement
--   in the n-direction, described by spherical polar angle theta
--   and azimuthal angle phi, would give angular momentum -hbar/2.
nm :: Double -> Double -> Ket
nm :: Double -> Double -> Ket
nm Double
theta Double
phi
    = (Double -> Double
forall a. Floating a => a -> a
sin (Double
theta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
zp
      Ket -> Ket -> Ket
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
cos (Double
theta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double -> Double
forall a. Floating a => a -> a
sin Double
phi) Complex Double -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
zm

-- | The orthonormal basis composed of 'xp' and 'xm'.
xBasis :: OrthonormalBasis
xBasis :: OrthonormalBasis
xBasis = [Ket] -> OrthonormalBasis
makeOB [Ket
xp,Ket
xm]

-- | The orthonormal basis composed of 'yp' and 'ym'.
yBasis :: OrthonormalBasis
yBasis :: OrthonormalBasis
yBasis = [Ket] -> OrthonormalBasis
makeOB [Ket
yp,Ket
ym]

-- | The orthonormal basis composed of 'zp' and 'zm'.
zBasis :: OrthonormalBasis
zBasis :: OrthonormalBasis
zBasis = [Ket] -> OrthonormalBasis
makeOB [Ket
zp,Ket
zm]

-- | Given spherical polar angle theta and azimuthal angle phi,
--   the orthonormal basis composed of 'np' theta phi and 'nm' theta phi.
nBasis :: Double -> Double -> OrthonormalBasis
nBasis :: Double -> Double -> OrthonormalBasis
nBasis Double
theta Double
phi = [Ket] -> OrthonormalBasis
makeOB [Double -> Double -> Ket
np Double
theta Double
phi,Double -> Double -> Ket
nm Double
theta Double
phi]

-- | The Pauli X operator.
sx :: Operator
sx :: Operator
sx = Ket
xp Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
xp Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
- Ket
xm Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
xm

-- | The Pauli Y operator.
sy :: Operator
sy :: Operator
sy = Ket
yp Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
yp Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
- Ket
ym Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
ym

-- | The Pauli Z operator.
sz :: Operator
sz :: Operator
sz = Ket
zp Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
zp Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
- Ket
zm Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
zm

-- | Pauli operator for an arbitrary direction given
--   by spherical coordinates theta and phi.
sn :: Double -> Double -> Operator
sn :: Double -> Double -> Operator
sn Double
theta Double
phi
    = (Double -> Double
forall a. Floating a => a -> a
sin Double
theta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Operator -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Operator
sx Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
+
      (Double -> Double
forall a. Floating a => a -> a
sin Double
theta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Operator -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Operator
sy Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
+
      (Double -> Double
forall a. Floating a => a -> a
cos Double
theta           Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Operator -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Operator
sz

-- | Alternative definition
--   of Pauli operator for an arbitrary direction.
sn' :: Double -> Double -> Operator
sn' :: Double -> Double -> Operator
sn' Double
theta Double
phi
    = Double -> Double -> Ket
np Double
theta Double
phi Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger (Double -> Double -> Ket
np Double
theta Double
phi) Operator -> Operator -> Operator
forall a. Num a => a -> a -> a
-
      Double -> Double -> Ket
nm Double
theta Double
phi Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger (Double -> Double -> Ket
nm Double
theta Double
phi)

----------------------
-- Quantum Dynamics --
----------------------

-- | Given a time step and a Hamiltonian operator,
--   produce a unitary time evolution operator.
--   Unless you really need the time evolution operator,
--   it is better to use 'timeEv', which gives the
--   same numerical results without doing an explicit
--   matrix inversion.  The function assumes hbar = 1.
timeEvOp :: Double -> Operator -> Operator
timeEvOp :: Double -> Operator -> Operator
timeEvOp Double
dt (Operator Matrix (Complex Double)
m) = Matrix (Complex Double) -> Operator
Operator (Double -> Matrix (Complex Double) -> Matrix (Complex Double)
M.timeEvMat Double
dt Matrix (Complex Double)
m)

-- | Given a time step and a Hamiltonian operator,
--   advance the state ket using the Schrodinger equation.
--   This method should be faster than using 'timeEvOp'
--   since it solves a linear system rather than calculating
--   an inverse matrix.  The function assumes hbar = 1.
timeEv :: Double -> Operator -> Ket -> Ket
timeEv :: Double -> Operator -> Ket -> Ket
timeEv Double
dt (Operator Matrix (Complex Double)
m) (Ket Vector (Complex Double)
k) = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double) -> Ket) -> Vector (Complex Double) -> Ket
forall a b. (a -> b) -> a -> b
$ Double
-> Matrix (Complex Double)
-> Vector (Complex Double)
-> Vector (Complex Double)
M.timeEv Double
dt Matrix (Complex Double)
m Vector (Complex Double)
k

-----------------
-- Composition --
-----------------

class Kron a where
    kron :: a -> a -> a

instance Kron Ket where
    kron :: Ket -> Ket -> Ket
kron (Ket Vector (Complex Double)
v1) (Ket Vector (Complex Double)
v2) = Vector (Complex Double) -> Ket
Ket (Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Kronecker a => a -> a -> a
M.kron Vector (Complex Double)
v1 Vector (Complex Double)
v2)

instance Kron Bra where
    kron :: Bra -> Bra -> Bra
kron (Bra Vector (Complex Double)
v1) (Bra Vector (Complex Double)
v2) = Vector (Complex Double) -> Bra
Bra (Vector (Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a. Kronecker a => a -> a -> a
M.kron Vector (Complex Double)
v1 Vector (Complex Double)
v2)

instance Kron Operator where
    kron :: Operator -> Operator -> Operator
kron (Operator Matrix (Complex Double)
m1) (Operator Matrix (Complex Double)
m2) = Matrix (Complex Double) -> Operator
Operator (Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
forall a. Kronecker a => a -> a -> a
M.kron Matrix (Complex Double)
m1 Matrix (Complex Double)
m2)

-----------------
-- Measurement --
-----------------

-- | The possible outcomes of a measurement
--   of an observable.
--   These are the eigenvalues of the operator
--   of the observable.
possibleOutcomes :: Operator -> [Double]
possibleOutcomes :: Operator -> [Double]
possibleOutcomes (Operator Matrix (Complex Double)
observable) = Matrix (Complex Double) -> [Double]
M.possibleOutcomes Matrix (Complex Double)
observable

-- | Given an obervable, return a list of pairs
--   of possible outcomes and projectors
--   for each outcome.
outcomesProjectors :: Operator -> [(Double,Operator)]
outcomesProjectors :: Operator -> [(Double, Operator)]
outcomesProjectors (Operator Matrix (Complex Double)
m)
    = [(Double
val,Matrix (Complex Double) -> Operator
Operator Matrix (Complex Double)
p) | (Double
val,Matrix (Complex Double)
p) <- Matrix (Complex Double) -> [(Double, Matrix (Complex Double))]
M.outcomesProjectors Matrix (Complex Double)
m]

-- | Given an observable and a state ket, return a list of pairs
--   of possible outcomes and probabilites
--   for each outcome.
outcomesProbabilities :: Operator -> Ket -> [(Double,Double)]
outcomesProbabilities :: Operator -> Ket -> [(Double, Double)]
outcomesProbabilities (Operator Matrix (Complex Double)
m) (Ket Vector (Complex Double)
v)
    = Matrix (Complex Double)
-> Vector (Complex Double) -> [(Double, Double)]
M.outcomesProbabilities Matrix (Complex Double)
m Vector (Complex Double)
v

{-
prob :: Ket -> Ket -> Double
prob k1 k2 = magnitude c ** 2
    where
      c = dagger k1 <> k2

probs :: OrthonormalBasis -> Ket -> [Double]
probs (OB ks) k = map (\bk -> let c = dagger bk <> k in magnitude c ** 2) ks
-}

{-
----------------------------------------
-- Angular Momentum of arbitrary size --
----------------------------------------

angularMomentumZMatrix :: Rational -> Matrix Cyclotomic
angularMomentumZMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> if r == c then fromRational (j + 1 - fromIntegral r) else 0)

twoJPlusOne :: Rational -> Maybe Int
twoJPlusOne j
    | j >= 0 && (denominator j == 1 || denominator j == 2)  = Just $ fromIntegral $ numerator (2 * j + 1)
    | otherwise                                             = Nothing

angularMomentumPlusMatrix :: Rational -> Matrix Cyclotomic
angularMomentumPlusMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let mr = j + 1 - fromIntegral r
                                             mc = j + 1 - fromIntegral c
                                         in if mr == mc + 1
                                            then sqrtRat (j*(j+1) - mc*mr)
                                            else 0
                              )

angularMomentumMinusMatrix :: Rational -> Matrix Cyclotomic
angularMomentumMinusMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let mr = j + 1 - fromIntegral r
                                             mc = j + 1 - fromIntegral c
                                         in if mr + 1 == mc
                                            then sqrtRat (j*(j+1) - mc*mr)
                                            else 0
                              )

angularMomentumXMatrix :: Rational -> Matrix Cyclotomic
angularMomentumXMatrix j
    = scaleMatrix (1/2) (angularMomentumPlusMatrix j + angularMomentumMinusMatrix j)

angularMomentumYMatrix :: Rational -> Matrix Cyclotomic
angularMomentumYMatrix j
    = scaleMatrix (1/(2*i)) (angularMomentumPlusMatrix j - angularMomentumMinusMatrix j)

jXMatrix :: Rational -> Matrix Cyclotomic
jXMatrix = angularMomentumXMatrix

jYMatrix :: Rational -> Matrix Cyclotomic
jYMatrix = angularMomentumYMatrix

jZMatrix :: Rational -> Matrix Cyclotomic
jZMatrix = angularMomentumZMatrix

matrixCommutator :: Matrix Cyclotomic -> Matrix Cyclotomic -> Matrix Cyclotomic
matrixCommutator m1 m2 = m1 * m2 - m2 * m1

-----------------------
-- Rotation matrices --
-----------------------

r2i :: Rational -> Integer
r2i r
    | denominator r == 1  = numerator r
    | otherwise           = error "r2i:  not an integer"

-- from Sakurai, revised, (3.8.33)
-- beta in degrees
smallDRotElement :: Rational -> Rational -> Rational -> Rational -> Cyclotomic
smallDRotElement j m' m beta
    = sum [parity(k-m+m') * sqrtRat (fact(j+m) * fact(j-m) * fact(j+m') * fact(j-m'))
                   / fromRational (fact(j+m-k) * fact(k) * fact(j-k-m') * fact(k-m+m'))
                         * cosDeg (beta/2) ^ r2i(2*j-2*k+m-m')
                         * sinDeg (beta/2) ^ r2i(2*k-m+m') | k <- [max 0 (m-m') .. min (j+m) (j-m')]]

parity :: Rational -> Cyclotomic
parity = fromIntegral . parityInteger . r2i

-- | (-1)^n, where n might be negative
parityInteger :: Integer -> Integer
parityInteger n
    | odd n      = -1
    | otherwise  =  1

factInteger :: Integer -> Integer
factInteger n
    | n == 0     = 1
    | n >  0     = n * factInteger (n-1)
    | otherwise  = error "factInteger called on negative argument"

fact :: Rational -> Rational
fact = fromIntegral . factInteger . r2i

-- | Rotation matrix elements.
--   From Sakurai, Revised Edition, (3.5.50).
--   The matrix desribes a rotation by gamma about the z axis,
--   followed by a rotation by beta about the y axis,
--   followed by a rotation by alpha about the z axis.
bigDRotElement :: Rational  -- ^ j, a nonnegative integer or half-integer
               -> Rational  -- ^ m', an integer or half-integer index for the row
               -> Rational  -- ^ m, an integer or half-integer index for the column
               -> Rational  -- ^ alpha, in degrees
               -> Rational  -- ^ beta, in degrees
               -> Rational  -- ^ gamma, in degrees
               -> Cyclotomic  -- ^ rotation matrix element
bigDRotElement j m' m alpha beta gamma
    = polarRat 1 (-(m' * alpha + m * gamma) / 360) * smallDRotElement j m' m beta

-- | Rotation matrix for a spin-j particle.
--   The matrix desribes a rotation by gamma about the z axis,
--   followed by a rotation by beta about the y axis,
--   followed by a rotation by alpha about the z axis.
rotationMatrix :: Rational  -- ^ j, a nonnegative integer or half-integer
               -> Rational  -- ^ alpha, in degrees
               -> Rational  -- ^ beta, in degrees
               -> Rational  -- ^ gamma, in degrees
               -> Matrix Cyclotomic  -- ^ rotation matrix
rotationMatrix j alpha beta gamma
    = case twoJPlusOne j of
        Nothing -> error "bigDRotMatrix: j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let m' = j + 1 - fromIntegral r
                                             m  = j + 1 - fromIntegral c
                                         in bigDRotElement j m' m alpha beta gamma
                              )

----------------------------------
-- Angular Momentum eigenstates --
----------------------------------

jmColumn :: Rational -> Rational -> Matrix Cyclotomic
jmColumn j m
    = case twoJPlusOne j of
        Nothing -> error "bigDRotMatrix: j must be a nonnegative integer or half-integer"
        Just d  -> matrix d 1 (\(r,_) -> let m' = j + 1 - fromIntegral r
                                         in if m == m'
                                            then 1
                                            else 0
                              )
-}

------------------
-- Gram-Schmidt --
------------------

-- | Form an orthonormal list of kets from
--   a list of linearly independent kets.
gramSchmidt :: [Ket] -> [Ket]
gramSchmidt :: [Ket] -> [Ket]
gramSchmidt [] = []
gramSchmidt [Ket
k] = [Ket -> Ket
forall a. HasNorm a => a -> a
normalize Ket
k]
gramSchmidt (Ket
k:[Ket]
ks) = let nks :: [Ket]
nks = [Ket] -> [Ket]
gramSchmidt [Ket]
ks
                         nk :: Ket
nk = Ket -> Ket
forall a. HasNorm a => a -> a
normalize ((Ket -> Ket -> Ket) -> Ket -> [Ket] -> Ket
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (-) Ket
k [Ket
w Ket -> Bra -> Operator
forall a b c. Mult a b c => a -> b -> c
<> Ket -> Bra
forall a b. Dagger a b => a -> b
dagger Ket
w Operator -> Ket -> Ket
forall a b c. Mult a b c => a -> b -> c
<> Ket
k | Ket
w <- [Ket]
nks])
                     in Ket
nkKet -> [Ket] -> [Ket]
forall a. a -> [a] -> [a]
:[Ket]
nks

-- todo:  Clebsch-Gordon coeffs