diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011-2015 Scott N. Walck . +Copyright (c) 2011-2016 Scott N. Walck . All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/examples/src/NMR.hs b/examples/src/NMR.hs new file mode 100644 --- /dev/null +++ b/examples/src/NMR.hs @@ -0,0 +1,16 @@ +{-# OPTIONS_GHC -Wall #-} + +-- ^ Nuclear Magnetic Resonance on the Bloch Sphere + +module Main where + +import Physics.Learn.QuantumMat + ( zm + ) +import Physics.Learn.BlochSphere + ( hamRabi + , evolutionBlochSphere + ) + +main :: IO () +main = evolutionBlochSphere zm (hamRabi 10 1 10) diff --git a/learn-physics.cabal b/learn-physics.cabal --- a/learn-physics.cabal +++ b/learn-physics.cabal @@ -1,10 +1,10 @@ Name: learn-physics -Version: 0.5.2 +Version: 0.6.0.0 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux, - magnetic field, and other quantities in mechanics - and electromagnetic theory. + magnetic field, and other quantities in classical mechanics, + electromagnetic theory, and quantum mechanics. License: BSD3 License-file: LICENSE Author: Scott N. Walck @@ -30,16 +30,23 @@ Physics.Learn.CompositeQuadrature Physics.Learn.RootFinding Physics.Learn.Mechanics + Physics.Learn.QuantumMat + Physics.Learn.Ket + Physics.Learn.BlochSphere + Physics.Learn.Schrodinger1D + Physics.Learn.BeamStack Physics.Learn Physics.Learn.Visual.PlotTools Physics.Learn.Visual.VisTools Physics.Learn.Visual.GlossTools - Build-depends: base >= 4.2 && < 4.9, + Build-depends: base >= 4.7 && < 4.9, vector-space >= 0.8.4 && < 0.11, - not-gloss >= 0.7.4 && < 0.8, - spatial-math >= 0.2 && < 0.3, - gloss >= 1.8 && < 1.10, - gnuplot >= 0.5 && < 0.6 + not-gloss >= 0.5.0.4 && < 0.8, + spatial-math >= 0.1.7 && < 0.3, + gloss >= 1.8, + gnuplot >= 0.5 && < 0.6, + linear >= 1.20, + hmatrix >= 0.17, polynomial >= 0.7 Hs-source-dirs: src Source-repository head @@ -73,13 +80,13 @@ Executable learn-physics-sunEarth Main-is: examples/src/sunEarthRK4.hs - Build-depends: gloss >= 1.8 && < 1.10, + Build-depends: gloss >= 1.8, base >= 4.5 && < 4.9, learn-physics Executable learn-physics-eFieldLine2D Main-is: examples/src/eFieldLine2D.hs - Build-depends: gloss >= 1.8 && < 1.10, + Build-depends: gloss >= 1.8, base >= 4.5 && < 4.9, learn-physics @@ -87,4 +94,9 @@ Main-is: examples/src/Projectile.hs Build-depends: gnuplot >= 0.5 && < 0.6, base >= 4.5 && < 4.9, + learn-physics + +Executable learn-physics-NMR + Main-is: examples/src/NMR.hs + Build-depends: base >= 4.5, learn-physics diff --git a/src/Physics/Learn/BeamStack.hs b/src/Physics/Learn/BeamStack.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/BeamStack.hs @@ -0,0 +1,286 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.BeamStack +Copyright : (c) Scott N. Walck 2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +Splitters, recombiners, and detectors for Stern-Gerlach +experiments. +-} + +-- Spin-1/2 mixed states. + +module Physics.Learn.BeamStack + ( + -- * Core laboratory components + BeamStack() + , randomBeam + , split + , recombine + , applyBField + , dropBeam + , flipBeams + , numBeams + , detect + -- * Standard splitters + , splitX + , splitY + , splitZ + -- * Standard magnetic fields + , applyBFieldX + , applyBFieldY + , applyBFieldZ + -- * Standard combiners + , recombineX + , recombineY + , recombineZ + -- * Filters + , xpFilter + , xmFilter + , zpFilter + , zmFilter + ) + where + +import Physics.Learn.QuantumMat + ( zp + , zm + , nm + , np + , couter + , oneQubitMixed + ) +import Numeric.LinearAlgebra + ( C + , Vector + , Matrix + , iC + , (<>) + , kronecker + , fromLists + , toList + , toLists + , scale + , size + , takeDiag + , ident + , tr + ) +import Data.Complex + ( Complex(..) + , realPart + , imagPart + ) +import Data.List + ( intercalate + ) + +data BeamStack = BeamStack (Matrix C) + +showOneBeam :: Double -> String +showOneBeam r = "Beam of intensity " ++ show r + +instance Show BeamStack where + show b = intercalate "\n" $ map showOneBeam (detect b) + +{- +unBeamStack :: BeamStack -> Matrix C +unBeamStack (BeamStack m) = m +-} + +-------------------- +-- Core functions -- +-------------------- + +-- | A beam of randomly oriented spin-1/2 particles. +randomBeam :: BeamStack +randomBeam = BeamStack oneQubitMixed + +extendWithZeros :: Matrix C -> Matrix C +extendWithZeros m + = let (_,q) = size m + ml = toLists m + in fromLists $ map (++ [0,0]) ml + ++ [replicate (q+2) 0, replicate (q+2) 0] + +-- reduce row and column size by 2 +reduceMat :: Matrix C -> Matrix C +reduceMat m + = let (p,q) = size m + ml = toLists m + in fromLists $ take (p-2) $ map (take (q-2)) ml + +checkedRealPart :: C -> Double +checkedRealPart c + = let eps = 1e-14 + in if imagPart c < eps + then realPart c + else error $ "checkRealPart: imagPart = " ++ show (imagPart c) + +-- | Return the intensities of a stack of beams. +detect :: BeamStack -> [Double] +detect (BeamStack m) + = addAlternate $ toList $ takeDiag m + +addAlternate :: [C] -> [Double] +addAlternate [] = [] +addAlternate [_] = error "addAlternate needs even number of elements" +addAlternate (x:y:xs) = checkedRealPart (x+y) : addAlternate xs + +-- | Remove the most recent beam from the stack. +dropBeam :: BeamStack -> BeamStack +dropBeam (BeamStack m) = BeamStack (reduceMat m) + +-- | Return the number of beams in a 'BeamStack'. +numBeams :: BeamStack -> Int +numBeams (BeamStack m) + = let (p,_) = size m + in p `div` 2 + +-- | Interchange the two most recent beams on the stack. +flipBeams :: BeamStack -> BeamStack +flipBeams (BeamStack m) + = let (d,_) = size m + fl = flipMat d + in BeamStack $ fl <> m <> tr fl + +flipMat :: Int -> Matrix C +flipMat d = bigM d (fromLists [[0,0,1,0] + ,[0,0,0,1] + ,[1,0,0,0] + ,[0,1,0,0]]) + +-- Turn a 2x2 into a dxd. +bigM2 :: Int -> Matrix C -> Matrix C +bigM2 d m + | d < 2 = error "bigM2 requires d >= 2" + | odd d = error "bigM2 requires even d" + | otherwise = fromLists $ map (++ [0,0]) (toLists (ident (d-2))) + ++ map (replicate (d-2) 0 ++) (toLists m) + +-- Turn a 4x4 into a dxd. +bigM :: Int -> Matrix C -> Matrix C +bigM d m + | d < 4 = error "bigM requires d >= 4" + | odd d = error "bigM requires even d" + | otherwise = fromLists $ map (++ [0,0,0,0]) (toLists (ident (d-4))) + ++ map (replicate (d-4) 0 ++) (toLists m) + +s :: Double -> Double -> Matrix C +s theta phi = kronecker (u `couter` u) (np theta phi `couter` np theta phi) + + kronecker (l `couter` u) (nm theta phi `couter` nm theta phi) + + kronecker (u `couter` l) (nm theta phi `couter` nm theta phi) + + kronecker (l `couter` l) (np theta phi `couter` np theta phi) + +u :: Vector C +u = zp + +l :: Vector C +l = zm + +-- | Given angles describing the orientation of the splitter, +-- removes an incoming beam from the stack and replaces +-- it with two beams, a spin-up and a spin-down beam. +-- The spin-down beam is the most recent beam on the stack. +split :: Double -> Double -> BeamStack -> BeamStack +split theta phi (BeamStack m) + = let m' = extendWithZeros m + (p,_) = size m' + ss = bigM p (s theta phi) + in BeamStack $ ss <> m' <> tr ss + +-- | Given angles describing the orientation of the recombiner, +-- returns a single beam from an incoming pair of beams. +recombine :: Double -> Double -> BeamStack -> BeamStack +recombine theta phi (BeamStack m) + = let (d,_) = size m + ss = bigM d (s theta phi) + in dropBeam $ BeamStack $ ss <> m <> tr ss + +mag2x2 :: Double -> Double -> Double -> Matrix C +mag2x2 theta phi omegaT + = let z = iC * (omegaT :+ 0) / 2 + np' = np theta phi + nm' = nm theta phi + in scale (exp z ) (np' `couter` np') + + scale (exp (-z)) (nm' `couter` nm') + +-- | Given angles describing the direction of a +-- uniform magnetic field, and given an angle +-- describing the product of the Larmor frequency +-- and the time, return an output beam from an +-- input beam. +applyBField :: Double -> Double -> Double -> BeamStack -> BeamStack +applyBField theta phi omegaT (BeamStack m) + = let (d,_) = size m + uu = bigM2 d (mag2x2 theta phi omegaT) + in BeamStack $ uu <> m <> tr uu + +----------------------- +-- Derived functions -- +----------------------- + +-- | A Stern-Gerlach splitter in the x direction. +splitX :: BeamStack -> BeamStack +splitX = split (pi/2) 0 + +-- | A Stern-Gerlach splitter in the y direction. +splitY :: BeamStack -> BeamStack +splitY = split (pi/2) (pi/2) + +-- | A Stern-Gerlach splitter in the z direction. +splitZ :: BeamStack -> BeamStack +splitZ = split 0 0 + +-- | Given an angle in radians +-- describing the product of the Larmor frequency +-- and the time, apply a magnetic in the x direction +-- to the most recent beam on the stack. +applyBFieldX :: Double -> BeamStack -> BeamStack +applyBFieldX = applyBField (pi/2) 0 + +-- | Given an angle in radians +-- describing the product of the Larmor frequency +-- and the time, apply a magnetic in the y direction +-- to the most recent beam on the stack. +applyBFieldY :: Double -> BeamStack -> BeamStack +applyBFieldY = applyBField (pi/2) (pi/2) + +-- | Given an angle in radians +-- describing the product of the Larmor frequency +-- and the time, apply a magnetic in the z direction +-- to the most recent beam on the stack. +applyBFieldZ :: Double -> BeamStack -> BeamStack +applyBFieldZ = applyBField 0 0 + +-- | A Stern-Gerlach recombiner in the x direction. +recombineX :: BeamStack -> BeamStack +recombineX = recombine (pi/2) 0 + +-- | A Stern-Gerlach recombiner in the y direction. +recombineY :: BeamStack -> BeamStack +recombineY = recombine (pi/2) (pi/2) + +-- | A Stern-Gerlach recombiner in the z direction. +recombineZ :: BeamStack -> BeamStack +recombineZ = recombine 0 0 + +-- | Filter for spin-up particles in the x direction. +xpFilter :: BeamStack -> BeamStack +xpFilter = dropBeam . splitX + +-- | Filter for spin-down particles in the x direction. +xmFilter :: BeamStack -> BeamStack +xmFilter = dropBeam . flipBeams . splitX + +-- | Filter for spin-up particles in the z direction. +zpFilter :: BeamStack -> BeamStack +zpFilter = dropBeam . splitZ + +-- | Filter for spin-down particles in the z direction. +zmFilter :: BeamStack -> BeamStack +zmFilter = dropBeam . flipBeams . splitZ diff --git a/src/Physics/Learn/BlochSphere.hs b/src/Physics/Learn/BlochSphere.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/BlochSphere.hs @@ -0,0 +1,176 @@ +{-# OPTIONS_GHC -Wall #-} + +{- | +Module : Physics.Learn.BlochSphere +Copyright : (c) Scott N. Walck 2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for displaying the +state of a spin-1/2 particle or other quantum two-level +system as a point on the Bloch Sphere. +-} + +module Physics.Learn.BlochSphere + ( VisObj + , toPos + , staticBlochSphere + , displayStaticState + , animatedBlochSphere + , simulateBlochSphere + , stateProp + , evolutionBlochSphere + , hamRabi + ) + where + +import Physics.Learn.QuantumMat + ( sy + , sz + , xp + , yp + , ym + , zm + , timeEv + ) +import Numeric.LinearAlgebra + ( Vector + , Matrix + , C + , iC +-- , (<>) -- matrix multiplication + , (|>) -- vector definition + , (!) -- vector element access + , (><) -- matrix definition + , scale + , size + ) +import Data.Complex + ( Complex(..) + , conjugate + , realPart + , imagPart + ) +import Physics.Learn + ( Position + , v3FromPos + , cart + ) +import SpatialMath + ( Euler(..) + ) +import Vis + ( VisObject(..) + , Flavour(..) + , Options(..) + , Camera0(..) + , defaultOpts + , display + , simulate + , blue + , red + ) + +-- | A Vis object. +type VisObj = VisObject Double + +-- | Convert a 2x1 complex state vector for a qubit +-- into Bloch (x,y,z) coordinates. +toPos :: Vector C -> Position +toPos v + = if size v /= 2 + then error "toPos only for size 2 vectors" + else let z1 = v ! 0 + z2 = v ! 1 + in cart (2 * realPart (conjugate z1 * z2)) + (2 * imagPart (conjugate z1 * z2)) + (realPart (conjugate z1 * z1 - conjugate z2 * z2)) + +-- | A static 'VisObj' for the state of a qubit. +staticBlochSphere :: Position -> VisObj +staticBlochSphere r + = RotEulerDeg (Euler 270 0 0) $ RotEulerDeg (Euler 0 180 0) $ + VisObjects [ Sphere 1 Wireframe blue + , Trans (v3FromPos r) (Sphere 0.05 Solid red) + ] + +displayStaticBlochSphere :: Position -> IO () +displayStaticBlochSphere r + = display myOptions (staticBlochSphere r) + +-- | Display a qubit state vector as a point on the Bloch Sphere. +displayStaticState :: Vector C -> IO () +displayStaticState = displayStaticBlochSphere . toPos + +displayxp :: IO () +displayxp = displayStaticState xp + +displayyp :: IO () +displayyp = displayStaticState yp + +displayym :: IO () +displayym = displayStaticState ym + +-- | Given a Bloch vector as a function of time, +-- return a 'VisObj' as a function of time. +animatedBlochSphere :: (Double -> Position) -> (Float -> VisObj) +animatedBlochSphere f + = staticBlochSphere . f . realToFrac + +-- | Given a sample rate, initial qubit state vector, and +-- state propagation function, produce a simulation. +-- The 'Float' in the state propagation function is the time +-- since the beginning of the simulation. +simulateBlochSphere :: Double -> Vector C -> (Float -> (Float,Vector C) -> (Float,Vector C)) -> IO () +simulateBlochSphere sampleRate initial statePropFunc + = simulate myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc + +{- +-- | Given a sample rate, initial qubit state vector, and +-- state propagation function, produce a simulation. +-- The 'Float' in the state propagation function is the time +-- since the beginning of the simulation. +playBlochSphere :: Double -> Vector C -> (Float -> (Float,Vector C) -> (Float,Vector C)) -> IO () +playBlochSphere sampleRate initial statePropFunc + = play myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc +-} + +-- | Produce a state propagation function from a time-dependent Hamiltonian. +stateProp :: (Double -> Matrix C) -> Float -> (Float,Vector C) -> (Float,Vector C) +stateProp ham tNew (tOld,v) + = (tNew, timeEv (realToFrac dt) (ham tMid) v) + where + dt = tNew - tOld + tMid = realToFrac $ (tNew + tOld) / 2 + +-- | Given an initial qubit state and a time-dependent Hamiltonian, +-- produce a visualization. +evolutionBlochSphere :: Vector C -> (Double -> Matrix C) -> IO () +evolutionBlochSphere psi0 ham + = simulateBlochSphere 0.01 psi0 (stateProp ham) + +myOptions :: Options +myOptions = defaultOpts {optWindowName = "Bloch Sphere" + ,optInitialCamera = Just (Camera0 75 20 4)} + +staticBz1 :: IO () +staticBz1 = evolutionBlochSphere xp (const (scale 0.9 sz)) + +staticBz2 :: IO () +staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const sz) + +staticBy1 :: IO () +staticBy1 = evolutionBlochSphere xp (const sy) + +-- | Hamiltonian for nuclear magnetic resonance. +-- Explain omega0, omegaR, omega. +hamRabi :: Double -> Double -> Double -> Double -> Matrix C +hamRabi omega0 omegaR omega t + = let h11 = omega0 :+ 0 + h12 = (omegaR :+ 0) * exp (-iC * ((omega * t) :+ 0)) + in scale (1/2) $ (2><2) [h11, h12, (conjugate h12), (-h11)] + +-- need to scale time + +-- a pi pulse diff --git a/src/Physics/Learn/Ket.hs b/src/Physics/Learn/Ket.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Ket.hs @@ -0,0 +1,492 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Safe #-} +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} +{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-} + +{- | +Module : Physics.Learn.Ket +Copyright : (c) Scott N. Walck 2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +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 + ( Ket + , Bra + , Operator + , Mult(..) + , Dagger(..) + , Representable(..) + , OrthonormalBasis + , makeOB + , listBasis + , size + , xp + , xm + , yp + , ym + , zp + , zm + , np + , nm + , xBasis + , yBasis + , zBasis + , sx + , sy + , sz + , prob + , probs + -- , 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 Data.Complex + ( Complex(..) + , magnitude + , conjugate + ) +import qualified Physics.Learn.QuantumMat as M +import Physics.Learn.QuantumMat + ( C + , Vector + , Matrix + , (#>) + , (<#) + , couter + , conjugateTranspose + , scaleV + , scaleM + , conjV + , fromList + , toList + , fromLists + ) + +infixl 7 <> + +-- | A ket vector describes the state of a quantum system. +data Ket = Ket (Vector C) + +instance Show Ket where + show k = + let message = "Use 'rep '." + in if dim k == 2 + then "Representation in zBasis:\n" ++ + show (rep zBasis k) ++ "\n" ++ message + else 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 _ = "\nTry 'rep zBasis '" + +-- | A bra vector describes the state of a quantum system. +data Bra = Bra (Vector C) + +instance Show Bra where + show _ = "\nTry 'rep zBasis '" + +-- | 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 + z1 <> z2 = z1 * z2 + +instance Mult C Ket Ket where + c <> Ket matrixKet = Ket (scaleV c matrixKet) + +instance Mult C Bra Bra where + c <> Bra matrixBra = Bra (scaleV c matrixBra) + +instance Mult C Operator Operator where + c <> Operator m = Operator (scaleM c m) + +instance Mult Ket C Ket where + Ket matrixKet <> c = Ket (scaleV c matrixKet) + +instance Mult Bra C Bra where + Bra matrixBra <> c = Bra (scaleV c matrixBra) + +instance Mult Operator C Operator where + Operator m <> c = Operator (scaleM c m) + +instance Mult Bra Ket C where + Bra matrixBra <> Ket matrixKet + = sum $ zipWith (*) (toList matrixBra) (toList matrixKet) + +instance Mult Bra Operator Bra where + Bra matrixBra <> Operator matrixOp + = Bra (matrixBra <# matrixOp) + +instance Mult Operator Ket Ket where + Operator matrixOp <> Ket matrixKet + = Ket (matrixOp #> matrixKet) + +instance Mult Ket Bra Operator where + Ket k <> Bra b = Operator (couter k b) + +instance Mult Operator Operator Operator where + Operator m1 <> Operator m2 = Operator (m1 M.<> m2) + +instance Num Ket where + Ket v1 + Ket v2 = Ket (v1 + v2) + Ket v1 - Ket v2 = Ket (v1 - v2) + (*) = error "Multiplication is not defined for kets" + negate (Ket v) = Ket (negate v) + abs = error "abs is not defined for kets" + signum = error "signum is not defined for kets" + fromInteger = error "fromInteger is not defined for kets" + +instance Num Bra where + Bra v1 + Bra v2 = Bra (v1 + v2) + Bra v1 - Bra v2 = Bra (v1 - v2) + (*) = error "Multiplication is not defined for bra vectors" + negate (Bra v) = Bra (negate v) + abs = error "abs is not defined for bra vectors" + signum = error "signum is not defined for bra vectors" + fromInteger = error "fromInteger is not defined for bra vectors" + +instance Num Operator where + Operator v1 + Operator v2 = Operator (v1 + v2) + Operator v1 - Operator v2 = Operator (v1 - v2) + Operator v1 * Operator v2 = Operator (v1 M.<> v2) + negate (Operator v) = Operator (negate v) + abs = error "abs is not defined for operators" + signum = error "signum is not defined for operators" + fromInteger = error "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 v) = Bra (conjV v) + +instance Dagger Bra Ket where + dagger (Bra v) = Ket (conjV v) + +instance Dagger Operator Operator where + dagger (Operator m) = Operator (conjugateTranspose m) + +instance Dagger C C where + dagger c = conjugate c + +class HasNorm a where + norm :: a -> Double + normalize :: a -> a + +instance HasNorm Ket where + norm (Ket v) = M.norm v + normalize k = (1 / norm k :+ 0) <> k + +instance HasNorm Bra where + norm (Bra v) = M.norm v + normalize b = (1 / norm b :+ 0) <> b + +{- +class HasDim a where + dim :: a -> Int + +instance HasDim Ket where + dim (Ket v) = M.dim v + +instance HasDim Bra where + dim (Bra v) = M.dim v +-} + +-- | An orthonormal basis of kets. +newtype OrthonormalBasis = OB [Ket] + deriving (Show) + +-- | Make an orthonormal basis from a list of linearly independent kets. +makeOB :: [Ket] -> OrthonormalBasis +makeOB = OB . gramSchmidt + +size :: OrthonormalBasis -> Int +size (OB ks) = length ks + +listBasis :: OrthonormalBasis -> [Ket] +listBasis (OB ks) = 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 (OB ks) k = fromList (map (\bk -> dagger bk <> k) ks) + dim (Ket v) = M.dim v + +instance Representable Bra (Vector C) where + rep (OB ks) b = fromList (map (\bk -> b <> bk) ks) + dim (Bra v) = M.dim v + +instance Representable Operator (Matrix C) where + rep (OB ks) op = fromLists [[ dagger k1 <> op <> k2 | k2 <- ks ] | k1 <- ks ] + dim (Operator m) = let (p,q) = M.size m + in if p == q then p else error "dim: non-square operator" + +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 + +-------------- +-- 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 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 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 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 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 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 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 theta phi + = (cos (theta / 2) :+ 0) <> zp + + (sin (theta / 2) :+ 0) * (cos phi :+ sin phi) <> 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 theta phi + = (sin (theta / 2) :+ 0) <> zp + - (cos (theta / 2) :+ 0) * (cos phi :+ sin phi) <> zm + +xBasis :: OrthonormalBasis +xBasis = makeOB [xp,xm] + +yBasis :: OrthonormalBasis +yBasis = makeOB [yp,ym] + +zBasis :: OrthonormalBasis +zBasis = makeOB [zp,zm] + +-- | The Pauli X operator. +sx :: Operator +sx = xp <> dagger xp - xm <> dagger xm + +-- | The Pauli Y operator. +sy :: Operator +sy = yp <> dagger yp - ym <> dagger ym + +-- | The Pauli Z operator. +sz :: Operator +sz = zp <> dagger zp - zm <> dagger zm + +{- +---------------------------------------- +-- 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 [] = [] +gramSchmidt [k] = [normalize k] +gramSchmidt (k:ks) = let nks = gramSchmidt ks + nk = normalize (foldl (-) k [w <> dagger w <> k | w <- nks]) + in nk:nks + +-- todo: Clebsch-Gordon coeffs diff --git a/src/Physics/Learn/QuantumMat.hs b/src/Physics/Learn/QuantumMat.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/QuantumMat.hs @@ -0,0 +1,328 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.QuantumMat +Copyright : (c) Scott N. Walck 2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains state vectors and matrices +for quantum mechanics. +-} + +-- Using only Complex Double here, no cyclotomic + +module Physics.Learn.QuantumMat + ( + -- * Complex numbers + C + -- * State Vectors + , xp + , xm + , yp + , ym + , zp + , zm + , np + , nm + , dim + , scaleV + , inner + , norm + , normalize + , probVector + , gramSchmidt + , conjV + , fromList + , toList + -- * Matrices (operators) + , sx + , sy + , sz + , scaleM + , (<>) + , (#>) + , (<#) + , conjugateTranspose + , fromLists + , toLists + , size + -- * Density matrices + , couter + , dm + , trace + , normalizeDM + , oneQubitMixed + -- * Quantum Dynamics + , timeEv + , timeEvMat + -- * Measurement + , possibleOutcomes + -- * Vector and Matrix + , Vector + , Matrix + ) + where + +import Numeric.LinearAlgebra + ( C + , Vector + , Matrix + , iC -- square root of negative one + , (><) -- matrix definition + , ident + , scale + , norm_2 + , inv + , (<\>) + , sym + , eigenvaluesSH + , cmap + , takeDiag + , conj + , dot + , tr + ) +-- , (<>) -- matrix product (not * !!!!) +-- , (#>) -- matrix-vector product +-- , fromList -- vector definition + +import qualified Numeric.LinearAlgebra as H +-- because H.outer does not conjugate +import Data.Complex + ( Complex(..) + , magnitude + ) + +-- | The state resulting from a measurement of +-- spin angular momentum in the x direction +-- on a spin-1/2 particle +-- when the result of the measurement is hbar/2. +xp :: Vector C +xp = normalize $ fromList [1, 1] + +-- | The state resulting from a measurement of +-- spin angular momentum in the x direction +-- on a spin-1/2 particle +-- when the result of the measurement is -hbar/2. +xm :: Vector C +xm = normalize $ fromList [1, -1] + +-- | The state resulting from a measurement of +-- spin angular momentum in the y direction +-- on a spin-1/2 particle +-- when the result of the measurement is hbar/2. +yp :: Vector C +yp = normalize $ fromList [1, iC] + +-- | The state resulting from a measurement of +-- spin angular momentum in the y direction +-- on a spin-1/2 particle +-- when the result of the measurement is -hbar/2. +ym :: Vector C +ym = normalize $ fromList [1, -iC] + +-- | The state resulting from a measurement of +-- spin angular momentum in the z direction +-- on a spin-1/2 particle +-- when the result of the measurement is hbar/2. +zp :: Vector C +zp = normalize $ fromList [1, 0] + +-- | The state resulting from a measurement of +-- spin angular momentum in the z direction +-- on a spin-1/2 particle +-- when the result of the measurement is -hbar/2. +zm :: Vector C +zm = normalize $ fromList [0, 1] + +-- | The state resulting from a measurement of +-- spin angular momentum in the direction +-- specified by spherical angles theta (polar angle) +-- and phi (azimuthal angle) +-- on a spin-1/2 particle +-- when the result of the measurement is hbar/2. +np :: Double -> Double -> Vector C +np theta phi = fromList [ cos (theta/2) :+ 0 + , exp(0 :+ phi) * (sin (theta/2) :+ 0) ] + +-- | The state resulting from a measurement of +-- spin angular momentum in the direction +-- specified by spherical angles theta (polar angle) +-- and phi (azimuthal angle) +-- on a spin-1/2 particle +-- when the result of the measurement is -hbar/2. +nm :: Double -> Double -> Vector C +nm theta phi = fromList [ sin (theta/2) :+ 0 + , -exp(0 :+ phi) * (cos (theta/2) :+ 0) ] + +-- | Dimension of a vector. +dim :: Vector C -> Int +dim = H.size + +-- | Scale a complex vector by a complex number. +scaleV :: C -> Vector C -> Vector C +scaleV = scale + +-- | Complex inner product. First vector gets conjugated. +inner :: Vector C -> Vector C -> C +inner = dot + +-- | Length of a complex vector. +norm :: Vector C -> Double +norm = norm_2 + +-- | Return a normalized version of a given state vector. +normalize :: Vector C -> Vector C +normalize v = scale (1 / norm_2 v :+ 0) v + +-- | Return a vector of probabilities for a given state vector. +probVector :: Vector C -- ^ state vector + -> Vector Double -- ^ vector of probabilities +probVector = cmap (\c -> magnitude c**2) + +-- | Conjugate the entries of a vector. +conjV :: Vector C -> Vector C +conjV = conj + +-- | Construct a vector from a list of complex numbers. +fromList :: [C] -> Vector C +fromList = H.fromList + +-- | Produce a list of complex numbers from a vector. +toList :: Vector C -> [C] +toList = H.toList + +-- | The Pauli X matrix. +sx :: Matrix C +sx = (2><2) [ 0, 1 + , 1, 0 ] + +-- | The Pauli Y matrix. +sy :: Matrix C +sy = (2><2) [ 0, -iC + , iC, 0 ] + +-- | The Pauli Z matrix. +sz :: Matrix C +sz = (2><2) [ 1, 0 + , 0, -1 ] + +-- | Scale a complex matrix by a complex number. +scaleM :: C -> Matrix C -> Matrix C +scaleM = scale + +-- | Matrix product. +(<>) :: Matrix C -> Matrix C -> Matrix C +(<>) = (H.<>) + +-- | Matrix-vector product. +(#>) :: Matrix C -> Vector C -> Vector C +(#>) = (H.#>) + +-- | Vector-matrix product +(<#) :: Vector C -> Matrix C -> Vector C +(<#) = (H.<#) + +-- | Conjugate transpose of a matrix. +conjugateTranspose :: Matrix C -> Matrix C +conjugateTranspose = tr + +-- | Construct a matrix from a list of lists of complex numbers. +fromLists :: [[C]] -> Matrix C +fromLists = H.fromLists + +-- | Produce a list of lists of complex numbers from a matrix. +toLists :: Matrix C -> [[C]] +toLists = H.toLists + +-- | Size of a matrix. +size :: Matrix C -> (Int,Int) +size = H.size + +---------------------- +-- Density Matrices -- +---------------------- + +-- | Complex outer product +couter :: Vector C -> Vector C -> Matrix C +couter v w = v `H.outer` conj w + +-- | Build a pure-state density matrix from a state vector. +dm :: Vector C -> Matrix C +dm cvec = cvec `couter` cvec + +-- | Trace of a matrix. +trace :: Matrix C -> C +trace = sum . toList . takeDiag + +-- | Normalize a density matrix so that it has trace one. +normalizeDM :: Matrix C -> Matrix C +normalizeDM rho = scale (1 / trace rho) rho + +-- | The one-qubit totally mixed state. +oneQubitMixed :: Matrix C +oneQubitMixed = normalizeDM $ ident 2 + +---------------------- +-- Quantum Dynamics -- +---------------------- + +-- | Given a time step and a Hamiltonian matrix, +-- produce a unitary time evolution matrix. +-- Unless you really need the time evolution matrix, +-- it is better to use 'timeEv', which gives the +-- same numerical results with doing an explicit +-- matrix inversion. The function assumes hbar = 1. +timeEvMat :: Double -> Matrix C -> Matrix C +timeEvMat dt h + = let ah = scale (0 :+ dt / 2) h + (l,m) = size h + n = if l == m then m else error "timeEv needs square Hamiltonian" + identity = ident n + in inv (identity + ah) <> (identity - ah) + +-- | Given a time step and a Hamiltonian matrix, +-- advance the state vector using the Schrodinger equation. +-- This method should be faster than using 'timeEvMat' +-- since it solves a linear system rather than calculating +-- an inverse matrix. The function assumes hbar = 1. +timeEv :: Double -> Matrix C -> Vector C -> Vector C +timeEv dt h v + = let ah = scale (0 :+ dt / 2) h + (l,m) = size h + n = if l == m then m else error "timeEv needs square Hamiltonian" + identity = ident n + in (identity + ah) <\> ((identity - ah) #> v) + +----------------- +-- Measurement -- +----------------- + +-- | The possible outcomes of a measurement +-- of an observable. +-- These are the eigenvalues of the matrix +-- of the observable. +possibleOutcomes :: Matrix C -> [Double] +possibleOutcomes observable + = H.toList $ eigenvaluesSH (sym observable) + +------------------ +-- Gram-Schmidt -- +------------------ + +-- | Form an orthonormal list of complex vectors +-- from a linearly independent list of complex vectors. +gramSchmidt :: [Vector C] -> [Vector C] +gramSchmidt [] = [] +gramSchmidt (v:vs) = let nvs = gramSchmidt vs + nv = normalize (v - sum [scale (inner w v) w | w <- nvs]) + in nv:nvs + +-- To Do +-- Generate higher spin operators and state vectors +-- eigenvectors +-- projection operators + diff --git a/src/Physics/Learn/Schrodinger1D.hs b/src/Physics/Learn/Schrodinger1D.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Schrodinger1D.hs @@ -0,0 +1,389 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE Trustworthy #-} + +{- | +Module : Physics.Learn.Schrodinger1D +Copyright : (c) Scott N. Walck 2015-2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions to +solve the (time dependent) Schrodinger equation +in one spatial dimension for a given potential function. +-} + +module Physics.Learn.Schrodinger1D + ( + -- * Potentials + freeV + , harmonicV + , squareWell + , doubleWell + , stepV + , wall + -- * Initial wavefunctions + , harm + , coherent + , gaussian + , movingGaussian + -- * Utilities + , stateVectorFromWavefunction + , hamiltonianMatrix + , expectX + , picture + ) + where + +import Data.Complex + ( Complex(..) + , magnitude + ) +import Graphics.Gloss + ( Picture(..) + , yellow + ) +import Math.Polynomial.Hermite + ( evalPhysHermite + ) +import Numeric.LinearAlgebra + ( R + , C + , Vector + , Matrix + , (|>) + , (<.>) + , fromLists + , toList + ) +import Physics.Learn.QuantumMat + ( probVector + , timeEv + ) + +hbar :: Double +hbar = 1 + +--i :: Complex Double +--i = 0 :+ 1 + +---------------- +-- Potentials -- +---------------- + +-- | Free potential. +-- The potential energy is zero everywhere. +freeV + :: Double -- ^ position + -> Double -- ^ potential energy +freeV _x = 0 + +-- | Harmonic potential. +-- This is the potential energy of a linear spring. +harmonicV + :: Double -- ^ spring constant + -> Double -- ^ position + -> Double -- ^ potential energy +harmonicV k x = k * x**2 / 2 + +-- | A double well potential. +-- Potential energy is a quartic function of position +-- that gives two wells, each approximately harmonic +-- at the bottom of the well. +doubleWell + :: Double -- ^ width (for both wells and well separation) + -> Double -- ^ energy height of barrier between wells + -> Double -- ^ position + -> Double -- ^ potential energy +doubleWell a v0 x = v0 * ((x**2 - a**2)/a**2)**2 + +-- | Finite square well potential. +-- Potential is zero inside the well, +-- and constant outside the well. +-- Well is centered at the origin. +squareWell + :: Double -- ^ well width + -> Double -- ^ energy height of well + -> Double -- ^ position + -> Double -- ^ potential energy +squareWell l v0 x + | abs x < l/2 = 0 + | otherwise = v0 + +-- | A step barrier potential. +-- Potential is zero to left of origin. +stepV + :: Double -- ^ energy height of barrier (to the right of origin) + -> Double -- ^ position + -> Double -- ^ potential energy +stepV v0 x + | x < 0 = 0 + | otherwise = v0 + +-- | A potential barrier with thickness and height. +wall + :: Double -- ^ thickness of wall + -> Double -- ^ energy height of barrier + -> Double -- ^ position of center of barrier + -> Double -- ^ position + -> Double -- ^ potential energy +wall w v0 x0 x + | abs (x-x0) < w/2 = v0 + | otherwise = 0 + +--------------------------- +-- Initial wavefunctions -- +--------------------------- + +-- | Harmonic oscillator stationary state +harm :: Int -- ^ nonnegative integer n identifying stationary state + -> Double -- ^ x / sqrt(hbar/(m * omega)), i.e. position + -- in units of sqrt(hbar/(m * omega)) + -> C -- ^ complex amplitude +harm n u + = exp (-u**2/2) * evalPhysHermite n u / sqrt (2^n * fact n * sqrt pi) :+ 0 + +coherent + :: Double -- ^ mass of particle + -> Double -- ^ angular frequency + -> Complex Double -- ^ parameter z + -> Double -> Complex Double -- ^ wavefunction +coherent m omega z x + = ((m*omega/(pi*hbar))**0.25 * exp(-m*omega*x**2/(2*hbar)) :+ 0) + * exp(-z**2/2 + (sqrt(2*m*omega/hbar) * x :+ 0) * z) + +gaussian + :: Double -- ^ width parameter + -> Double -- ^ center of wave packet + -> Double -> Complex Double -- ^ wavefunction +gaussian a x0 x = exp(-(x-x0)**2/(2*a**2)) / sqrt(a * sqrt pi) :+ 0 + +movingGaussian + :: Double -- ^ width parameter + -> Double -- ^ center of wave packet + -> Double -- ^ momentum + -> Double -> Complex Double -- ^ wavefunction +movingGaussian a x0 p0 x = exp((0 :+ p0*x/hbar) - ((x-x0)**2/(2*a**2) :+ 0)) / (sqrt(a * sqrt pi) :+ 0) + +--------------- +-- Utilities -- +--------------- + +fact :: Int -> Double +fact 0 = 1 +fact n = fromIntegral n * fact (n-1) + +linspace :: Double -> Double -> Int -> [Double] +linspace left right num + = let dx = (right - left) / fromIntegral (num - 1) + in [ left + dx * fromIntegral n | n <- [0..num-1]] + +-- | Transform a wavefunction into a state vector. +stateVectorFromWavefunction :: R -- ^ lowest x + -> R -- ^ highest x + -> Int -- ^ dimension of state vector + -> (R -> C) -- ^ wavefunction + -> Vector C -- ^ state vector +stateVectorFromWavefunction left right num psi + = (num |>) [psi x | x <- linspace left right num] + +hamiltonianMatrix :: R -- ^ lowest x + -> R -- ^ highest x + -> Int -- ^ dimension of state vector + -> R -- ^ hbar + -> R -- ^ mass + -> (R -> R) -- ^ potential energy function + -> Matrix C -- ^ Hamiltonian Matrix +hamiltonianMatrix xmin xmax num hbar m pe + = let coeff = -hbar**2/(2*m) + dx = (xmax - xmin) / fromIntegral (num - 1) + diagKEterm = -2 * coeff / dx**2 + offdiagKEterm = coeff / dx**2 + xs = linspace xmin xmax num + in fromLists [[case abs(i-j) of + 0 -> (diagKEterm + pe x) :+ 0 + 1 -> offdiagKEterm :+ 0 + _ -> 0 + | j <- [1..num] ] | (i,x) <- zip [1..num] xs] + +expectX :: Vector C -- ^ state vector + -> Vector R -- ^ vector of x values + -> R -- ^ , expectation value of X +expectX psi xs = probVector psi <.> xs + + +glossScaleX :: Int -> (Double,Double) -> Double -> Float +glossScaleX screenWidth (xmin,xmax) x + = let w = fromIntegral screenWidth :: Double + in realToFrac $ (x - xmin) / (xmax - xmin) * w - w / 2 + +glossScaleY :: Int -> (Double,Double) -> Double -> Float +glossScaleY screenHeight (ymin,ymax) y + = let h = fromIntegral screenHeight :: Double + in realToFrac $ (y - ymin) / (ymax - ymin) * h - h / 2 + +glossScalePoint :: (Int,Int) -- ^ (screenWidth,screenHeight) + -> (Double,Double) -- ^ (xmin,xmax) + -> (Double,Double) -- ^ (ymin,ymax) + -> (Double,Double) -- ^ (x,y) + -> (Float,Float) +glossScalePoint (screenWidth,screenHeight) xMinMax yMinMax (x,y) + = (glossScaleX screenWidth xMinMax x + ,glossScaleY screenHeight yMinMax y) + + +-- | Produce a gloss 'Picture' of state vector +-- for 1D wavefunction. +picture :: (Double, Double) -- ^ y range + -> [Double] -- ^ xs + -> Vector C -- ^ state vector + -> Picture +picture (ymin,ymax) xs psi + = Color + yellow + (Line + [glossScalePoint + (screenWidth,screenHeight) + (head xs, last xs) + (ymin,ymax) + p | p <- zip xs (map magSq $ toList psi)]) + where + magSq = \z -> magnitude z ** 2 + screenWidth = 1000 + screenHeight = 750 + +{- +-- | Given an initial state vector and +-- state propagation function, produce a simulation. +-- The 'Float' in the state propagation function is the time +-- interval for one timestep. +simulate1D :: [Double] -> Vector C -> (Float -> (Float,[Double],Vector C) -> (Float,[Double],Vector C)) -> IO () +simulate1D xs initial statePropFunc + = simulate display black 10 (0,initial) displayFunc (const statePropFunc) + where + display = InWindow "Animation" (screenWidth,screenHeight) (10,10) + displayFunc (_t,v) = Color yellow (Line [( + + white (\tFloat -> Pictures [Color blue (Line (points (realToFrac tFloat))) + ,axes (screenWidth,screenHeight) (xmin,xmax) (ymin,ymax)]) + +-- | Produce a state propagation function from a time-dependent Hamiltonian. +-- The float is dt. +statePropGloss :: (Double -> Matrix C) -> Float -> (Float,Vector C) -> (Float,Vector C) +statePropGloss ham dt (tOld,v) + = (tNew, timeEv (realToFrac dt) (ham tMid) v) + where + tNew = tOld + dt + tMid = realToFrac $ (tNew + tOld) / 2 + +-- | Given an initial state vector and a time-dependent Hamiltonian, +-- produce a visualization of a 1D wavefunction. +evolutionBlochSphere :: Vector C -> (Double -> Matrix C) -> IO () +evolutionBlochSphere psi0 ham + = simulateBlochSphere 0.01 psi0 (stateProp ham) + +-} + + +{- +def triDiagMatrixMult(square_arr,arr): + num = len(arr) + result = array([0 for n in range(num)],dtype=complex128) + result[0] = square_arr[0][0] * arr[0] + square_arr[0][1] * arr[1] + for n in range(1,num-1): + result[n] = square_arr[n][n-1] * arr[n-1] + square_arr[n][n] * arr[n] \ + + square_arr[n][n+1] * arr[n+1] + result[num-1] = square_arr[num-1][num-2] * arr[num-2] \ + + square_arr[num-1][num-1] * arr[num-1] + return result + +################ +# Main program # +################ + +if __name__ == '__main__': + m = 1 + omega = 10 + xmin = -2.0 + xmax = 2.0 +# num = 256 + num = 128 + dt = 0.0002 +# dt = 0.01 + xs = linspace(xmin,xmax,num) + dx = xs[1] - xs[0] + + super = lambda x: (harm0(m,omega)(x) + harm1(m,omega)(x))/sqrt(2) + shiftedHarm = lambda x: harm0(m,omega)(x-1) + coh = coherent(m,omega,1) + +# print sum(conj(psi)*psi)*dx + + harmV = harmonicV(m * omega**2) + +# V = doubleWell(1,0.1*hbar*omega) + V = squareWell(1.0,hbar*omega) +# V = harmonicV(m*omega**2) +# V = stepV(10*hbar*omega) +# V = wall(0.1,14.0*hbar*omega,0) +# V = freeV + + H = matrixH(m,xmin,xmax,num,V) + I = matrixI(num) + + (vals,vecs) = eigh(H) + + E0 = vals[0] + E1 = vals[1] + psi0 = normalize(transpose(vecs)[0],dx) + psi1 = normalize(transpose(vecs)[1],dx) + +# psi = func2psi(gaussian(0.3,1),xmin,xmax,num) +# psi = func2psi(coh,xmin,xmax,num) +# psi = func2psi(movingGaussian(0.3,10,-1),xmin,xmax,num) + + psi = psi0 +# psi = psi1 +# psi = (psi0 + psi1)/sqrt(2) + + E = sum(conj(psi)*triDiagMatrixMult(H,psi)).real*dx + + Escale = hbar*omega + + print E + print Escale + + leftM = I + 0.5 * i * H / hbar * dt + rightM = I - 0.5 * i * H / hbar * dt + + box = display(title='Schrodinger Equation',width=1000,height=1000) + + c = curve(pos = psi2rho(psi,xs)) + c.color = color.blue + c.radius = 0.02 + + ball = sphere(radius=0.05,color=color.red,pos=(expectX(psi,xs),0,0)) + + pot_curve = [(x,V(x)/Escale,0) for x in xs if V(x)/Escale < xmax] + pot = curve(color=color.green,pos=pot_curve,radius=0.01) + + Eline = curve(color=(1,1,0),pos=[(x,E/Escale) for x in xs]) +# axis = curve(color=color.white,pos=[(x,0) for x in xs]) + + while 1: + psi = solve(leftM,triDiagMatrixMult(rightM,psi)) + c.pos = psi2rho(psi,xs) + ball.x = expectX(psi,xs) + +To Do: +add combinators for potentials +to shift horizontally and vertically, +and to add potentials + +-} + +-- Are we committed to SI units for hbar? +-- harmonic oscillator functions depend only on sqrt(hbar/m omega) +-- which is a length parameter +-- for moving gaussian, could give hbar/p0 instead of p0 +-- (is that debrogie wavelength? I think it's h/p0)