learn-physics 0.5.2 → 0.6.0.0

raw patch · 8 files changed

+1710/−11 lines, 8 filesdep +hmatrixdep +lineardep +polynomialdep ~basedep ~glossdep ~not-glossnew-component:exe:learn-physics-NMRPVP ok

version bump matches the API change (PVP)

Dependencies added: hmatrix, linear, polynomial

Dependency ranges changed: base, gloss, not-gloss, spatial-math

API changes (from Hackage documentation)

Files

• LICENSE +1/−1
• examples/src/NMR.hs +16/−0
• learn-physics.cabal +22/−10
• src/Physics/Learn/BeamStack.hs +286/−0
• src/Physics/Learn/BlochSphere.hs +176/−0
• src/Physics/Learn/Ket.hs +492/−0
• src/Physics/Learn/QuantumMat.hs +328/−0
• src/Physics/Learn/Schrodinger1D.hs +389/−0

@@ -1,4 +1,4 @@-Copyright (c) 2011-2015 Scott N. Walck <walck@lvc.edu>.+Copyright (c) 2011-2016 Scott N. Walck <walck@lvc.edu>. All rights reserved.  Redistribution and use in source and binary forms, with or without

@@ -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)

@@ -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

@@ -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 <walck@lvc.edu>+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

@@ -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 <walck@lvc.edu>+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

@@ -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 <walck@lvc.edu>+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 <basis name> <ket name>'."+        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 _ = "<operator>\nTry 'rep zBasis <operator name>'"++-- | A bra vector describes the state of a quantum system.+data Bra = Bra (Vector C)++instance Show Bra where+    show _ = "<bra>\nTry 'rep zBasis <bra name>'"++-- | 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

@@ -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 <walck@lvc.edu>+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+

@@ -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 <walck@lvc.edu>+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         -- ^ <X>, 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)