learn-physics 0.6.0.0 → 0.6.0.1
raw patch · 4 files changed
+303/−70 lines, 4 filesdep ~base
Dependency ranges changed: base
Files
- learn-physics.cabal +10/−10
- src/Physics/Learn/BlochSphere.hs +68/−22
- src/Physics/Learn/Ket.hs +150/−36
- src/Physics/Learn/QuantumMat.hs +75/−2
learn-physics.cabal view
@@ -1,5 +1,5 @@ Name: learn-physics-Version: 0.6.0.0+Version: 0.6.0.1 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux,@@ -12,7 +12,7 @@ Category: Physics Build-type: Simple Cabal-version: >=1.8-Tested-with: GHC == 7.10.2+Tested-with: GHC >= 7.10.2 && <= 7.10.3 Library Exposed-modules: Physics.Learn.Charge Physics.Learn.Current@@ -39,7 +39,7 @@ Physics.Learn.Visual.PlotTools Physics.Learn.Visual.VisTools Physics.Learn.Visual.GlossTools- Build-depends: base >= 4.7 && < 4.9,+ Build-depends: base >= 4.7 && < 4.10, vector-space >= 0.8.4 && < 0.11, not-gloss >= 0.5.0.4 && < 0.8, spatial-math >= 0.1.7 && < 0.3,@@ -56,44 +56,44 @@ Executable learn-physics-PlaneWave Main-is: examples/src/PlaneWave.hs Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-eFieldLine3D Main-is: examples/src/eFieldLine3D.hs Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-LorentzForceSimulation Main-is: examples/src/LorentzForceSimulation.hs Build-depends: not-gloss >= 0.7.4 && < 0.8, spatial-math >= 0.2 && < 0.3,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-BCircularLoop Main-is: examples/src/BCircularLoop.hs Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-sunEarth Main-is: examples/src/sunEarthRK4.hs Build-depends: gloss >= 1.8,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-eFieldLine2D Main-is: examples/src/eFieldLine2D.hs Build-depends: gloss >= 1.8,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-Projectile Main-is: examples/src/Projectile.hs Build-depends: gnuplot >= 0.5 && < 0.6,- base >= 4.5 && < 4.9,+ base >= 4.5 && < 4.10, learn-physics Executable learn-physics-NMR
src/Physics/Learn/BlochSphere.hs view
@@ -15,24 +15,27 @@ module Physics.Learn.BlochSphere ( VisObj , toPos+ , ketToPos , staticBlochSphere , displayStaticState , animatedBlochSphere , simulateBlochSphere+ , simulateBlochSphereK , stateProp+ , statePropK , evolutionBlochSphere+ , evolutionBlochSphereK , hamRabi ) where -import Physics.Learn.QuantumMat- ( sy- , sz- , xp- , yp- , ym- , zm- , timeEv+import qualified Physics.Learn.QuantumMat as M+import qualified Physics.Learn.Ket as K+import Physics.Learn.Ket+ ( Ket+ , Operator+ , (<>)+ , dagger ) import Numeric.LinearAlgebra ( Vector@@ -40,7 +43,7 @@ , C , iC -- , (<>) -- matrix multiplication- , (|>) -- vector definition+-- , (|>) -- vector definition , (!) -- vector element access , (><) -- matrix definition , scale@@ -72,6 +75,22 @@ , red ) +{-+3 ways to specify the state of a spin-1/2 particle:+Vector C+Ket+Position (Bloch vector)++2 ways to specify a Hamiltonian:+Matrix C+Operator++3 choices for Vis' world:+(Float, Vector C)+(Float, Ket)+(Float, Position)+-}+ -- | A Vis object. type VisObj = VisObject Double @@ -87,6 +106,18 @@ (2 * imagPart (conjugate z1 * z2)) (realPart (conjugate z1 * z1 - conjugate z2 * z2)) +-- | Convert a qubit ket+-- into Bloch (x,y,z) coordinates.+ketToPos :: Ket -> Position+ketToPos psi+ = if K.dim psi /= 2+ then error "ketToPos only for qubit kets"+ else let z1 = dagger K.zp <> psi+ z2 = dagger K.zm <> psi+ 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@@ -103,15 +134,6 @@ 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)@@ -126,6 +148,14 @@ simulateBlochSphere sampleRate initial statePropFunc = simulate myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc +-- | Given a sample rate, initial qubit state ket, and+-- state propagation function, produce a simulation.+-- The 'Float' in the state propagation function is the time+-- since the beginning of the simulation.+simulateBlochSphereK :: Double -> Ket -> (Float -> (Float,Ket) -> (Float,Ket)) -> IO ()+simulateBlochSphereK sampleRate initial statePropFuncK+ = simulate myOptions sampleRate (0,initial) (staticBlochSphere . ketToPos . snd) statePropFuncK+ {- -- | Given a sample rate, initial qubit state vector, and -- state propagation function, produce a simulation.@@ -139,29 +169,45 @@ -- | 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)+ = (tNew, M.timeEv (realToFrac dt) (ham tMid) v) where dt = tNew - tOld tMid = realToFrac $ (tNew + tOld) / 2 +-- | Produce a state propagation function from a time-dependent Hamiltonian.+statePropK :: (Double -> Operator) -> Float -> (Float,Ket) -> (Float,Ket)+statePropK ham tNew (tOld,psi)+ = (tNew, K.timeEv (realToFrac dt) (ham tMid) psi)+ 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) +-- | Given an initial qubit ket and a time-dependent Hamiltonian,+-- produce a visualization.+evolutionBlochSphereK :: Ket -> (Double -> Operator) -> IO ()+evolutionBlochSphereK psi0 ham+ = simulateBlochSphereK 0.01 psi0 (statePropK ham)+ myOptions :: Options myOptions = defaultOpts {optWindowName = "Bloch Sphere" ,optInitialCamera = Just (Camera0 75 20 4)} +{- staticBz1 :: IO ()-staticBz1 = evolutionBlochSphere xp (const (scale 0.9 sz))+staticBz1 = evolutionBlochSphere M.xp (const (scale 0.9 M.sz)) staticBz2 :: IO ()-staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const sz)+staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const M.sz) staticBy1 :: IO ()-staticBy1 = evolutionBlochSphere xp (const sy)+staticBy1 = evolutionBlochSphere M.xp (const M.sy)+-} -- | Hamiltonian for nuclear magnetic resonance. -- Explain omega0, omegaR, omega.
src/Physics/Learn/Ket.hs view
@@ -17,16 +17,15 @@ -- a Ket layer on top of QuantumMat module Physics.Learn.Ket- ( Ket+ (+ -- * Basic data types+ C+ , i+ , magnitude+ , Ket , Bra , Operator- , Mult(..)- , Dagger(..)- , Representable(..)- , OrthonormalBasis- , makeOB- , listBasis- , size+ -- * Kets for spin-1/2 particles , xp , xm , yp@@ -35,14 +34,41 @@ , zm , np , nm- , xBasis- , yBasis- , zBasis+ -- * Operators for spin-1/2 particles , sx , sy , sz- , prob- , probs+ , sn+ , sn'+ -- * Quantum Dynamics+ , timeEvOp+ , timeEv+ -- * Composition+ , Kron(..)+ -- * Measurement+ , possibleOutcomes+ , outcomesProjectors+ , outcomesProbabilities+-- , prob+-- , probs+ -- * Generic multiplication+ , Mult(..)+ -- * Adjoint operation+ , Dagger(..)+ -- * Normalization+ , HasNorm(..)+ -- * Representation+ , Representable(..)+ -- * Orthonormal bases+ , OrthonormalBasis+ , makeOB+ , listBasis+ , size+ -- * Orthonormal bases for spin-1/2 particles+ , xBasis+ , yBasis+ , zBasis+ , nBasis -- , angularMomentumXMatrix -- , angularMomentumYMatrix -- , angularMomentumZMatrix@@ -60,9 +86,9 @@ -- We try to import only from QuantumMat -- and not from Numeric.LinearAlgebra +import qualified Data.Complex as C import Data.Complex ( Complex(..)- , magnitude , conjugate ) import qualified Physics.Learn.QuantumMat as M@@ -72,7 +98,6 @@ , Matrix , (#>) , (<#)- , couter , conjugateTranspose , scaleV , scaleM@@ -100,7 +125,12 @@ data Operator = Operator (Matrix C) instance Show Operator where- show _ = "<operator>\nTry 'rep zBasis <operator name>'"+ show op =+ let message = "Use 'rep <basis name> <operator name>'."+ in if dim op == 2+ then "Representation in zBasis:\n" +++ show (rep zBasis op) ++ "\n" ++ message+ else message -- | A bra vector describes the state of a quantum system. data Bra = Bra (Vector C)@@ -108,6 +138,12 @@ instance Show Bra where show _ = "<bra>\nTry 'rep zBasis <bra name>'" +magnitude :: C -> Double+magnitude = C.magnitude++i :: C+i = 0 :+ 1+ -- | Generic multiplication including inner product, -- outer product, operator product, and whatever else makes sense. -- No conjugation takes place in this operation.@@ -148,7 +184,9 @@ = Ket (matrixOp #> matrixKet) instance Mult Ket Bra Operator where- Ket k <> Bra b = Operator (couter k b)+ Ket k <> Bra b =+ Operator+ (fromLists [[ x*y | y <- toList b] | x <- toList k]) instance Mult Operator Operator Operator where Operator m1 <> Operator m2 = Operator (m1 M.<> m2)@@ -209,17 +247,6 @@ 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)@@ -256,14 +283,6 @@ 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 -- --------------@@ -314,15 +333,23 @@ = (sin (theta / 2) :+ 0) <> zp - (cos (theta / 2) :+ 0) * (cos phi :+ sin phi) <> zm +-- | The orthonormal basis composed of 'xp' and 'xm'. xBasis :: OrthonormalBasis xBasis = makeOB [xp,xm] +-- | The orthonormal basis composed of 'yp' and 'ym'. yBasis :: OrthonormalBasis yBasis = makeOB [yp,ym] +-- | The orthonormal basis composed of 'zp' and 'zm'. zBasis :: OrthonormalBasis zBasis = makeOB [zp,zm] +-- | Given spherical polar angle theta and azimuthal angle phi,+-- the orthonormal basis composed of 'np' theta phi and 'nm' theta phi.+nBasis :: Double -> Double -> OrthonormalBasis+nBasis theta phi = makeOB [np theta phi,nm theta phi]+ -- | The Pauli X operator. sx :: Operator sx = xp <> dagger xp - xm <> dagger xm@@ -334,6 +361,93 @@ -- | The Pauli Z operator. sz :: Operator sz = zp <> dagger zp - zm <> dagger zm++-- | Pauli operator for an arbitrary direction given+-- by spherical coordinates theta and phi.+sn :: Double -> Double -> Operator+sn theta phi+ = (sin theta * cos phi :+ 0) <> sx ++ (sin theta * sin phi :+ 0) <> sy ++ (cos theta :+ 0) <> sz++-- | Alternative definition+-- of Pauli operator for an arbitrary direction.+sn' :: Double -> Double -> Operator+sn' theta phi+ = np theta phi <> dagger (np theta phi) -+ nm theta phi <> dagger (nm theta phi)++----------------------+-- Quantum Dynamics --+----------------------++-- | Given a time step and a Hamiltonian operator,+-- produce a unitary time evolution operator.+-- Unless you really need the time evolution operator,+-- it is better to use 'timeEv', which gives the+-- same numerical results without doing an explicit+-- matrix inversion. The function assumes hbar = 1.+timeEvOp :: Double -> Operator -> Operator+timeEvOp dt (Operator m) = Operator (M.timeEvMat dt m)++-- | Given a time step and a Hamiltonian operator,+-- advance the state ket using the Schrodinger equation.+-- This method should be faster than using 'timeEvOp'+-- since it solves a linear system rather than calculating+-- an inverse matrix. The function assumes hbar = 1.+timeEv :: Double -> Operator -> Ket -> Ket+timeEv dt (Operator m) (Ket k) = Ket $ M.timeEv dt m k++-----------------+-- Composition --+-----------------++class Kron a where+ kron :: a -> a -> a++instance Kron Ket where+ kron (Ket v1) (Ket v2) = Ket (M.kron v1 v2)++instance Kron Bra where+ kron (Bra v1) (Bra v2) = Bra (M.kron v1 v2)++instance Kron Operator where+ kron (Operator m1) (Operator m2) = Operator (M.kron m1 m2)++-----------------+-- Measurement --+-----------------++-- | The possible outcomes of a measurement+-- of an observable.+-- These are the eigenvalues of the operator+-- of the observable.+possibleOutcomes :: Operator -> [Double]+possibleOutcomes (Operator observable) = M.possibleOutcomes observable++-- | Given an obervable, return a list of pairs+-- of possible outcomes and projectors+-- for each outcome.+outcomesProjectors :: Operator -> [(Double,Operator)]+outcomesProjectors (Operator m)+ = [(val,Operator p) | (val,p) <- M.outcomesProjectors m]++-- | Given an observable and a state ket, return a list of pairs+-- of possible outcomes and probabilites+-- for each outcome.+outcomesProbabilities :: Operator -> Ket -> [(Double,Double)]+outcomesProbabilities (Operator m) (Ket v)+ = M.outcomesProbabilities m v++{-+prob :: Ket -> Ket -> Double+prob k1 k2 = magnitude c ** 2+ where+ c = dagger k1 <> k2++probs :: OrthonormalBasis -> Ket -> [Double]+probs (OB ks) k = map (\bk -> let c = dagger bk <> k in magnitude c ** 2) ks+-} {- ----------------------------------------
src/Physics/Learn/QuantumMat.hs view
@@ -49,6 +49,7 @@ , fromLists , toLists , size+ , matrixFunction -- * Density matrices , couter , dm@@ -56,10 +57,15 @@ , normalizeDM , oneQubitMixed -- * Quantum Dynamics- , timeEv , timeEvMat+ , timeEv+ , timeEvMatSpec+ -- * Composition+ , Kronecker(..) -- * Measurement , possibleOutcomes+ , outcomesProjectors+ , outcomesProbabilities -- * Vector and Matrix , Vector , Matrix@@ -70,6 +76,7 @@ ( C , Vector , Matrix+ , Herm , iC -- square root of negative one , (><) -- matrix definition , ident@@ -79,6 +86,7 @@ , (<\>) , sym , eigenvaluesSH+ , eigSH , cmap , takeDiag , conj@@ -94,6 +102,7 @@ import Data.Complex ( Complex(..) , magnitude+ , realPart ) -- | The state resulting from a measurement of@@ -195,6 +204,10 @@ toList :: Vector C -> [C] toList = H.toList +--------------+-- Matrices --+--------------+ -- | The Pauli X matrix. sx :: Matrix C sx = (2><2) [ 0, 1@@ -242,6 +255,15 @@ size :: Matrix C -> (Int,Int) size = H.size +-- | Apply a function to a matrix.+-- Assumes the matrix is a normal matrix (a matrix+-- with an orthonormal basis of eigenvectors).+matrixFunction :: (C -> C) -> Matrix C -> Matrix C+matrixFunction f m+ = let (valv,vecm) = H.eig m+ fvalv = fromList [f val | val <- toList valv]+ in vecm <> H.diag fvalv <> tr vecm+ ---------------------- -- Density Matrices -- ----------------------@@ -274,7 +296,7 @@ -- 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+-- same numerical results without doing an explicit -- matrix inversion. The function assumes hbar = 1. timeEvMat :: Double -> Matrix C -> Matrix C timeEvMat dt h@@ -297,7 +319,26 @@ identity = ident n in (identity + ah) <\> ((identity - ah) #> v) +-- | Given a Hamiltonian matrix, return a function from time+-- to evolution matrix. Uses spectral decomposition.+-- Assumes hbar = 1.+timeEvMatSpec :: Matrix C -> Double -> Matrix C+timeEvMatSpec m t = matrixFunction (\h -> exp(-iC * h * (t :+ 0))) m+ -----------------+-- Composition --+-----------------++class Kronecker a where+ kron :: a -> a -> a++instance H.Product t => Kronecker (Vector t) where+ kron v1 v2 = H.fromList [c1 * c2 | c1 <- H.toList v1, c2 <- H.toList v2]++instance H.Product t => Kronecker (Matrix t) where+ kron = H.kronecker++----------------- -- Measurement -- ----------------- @@ -308,6 +349,38 @@ possibleOutcomes :: Matrix C -> [Double] possibleOutcomes observable = H.toList $ eigenvaluesSH (sym observable)++-- From a Hermitian matrix, a list of pairs of eigenvalues and eigenvectors.+valsVecs :: Herm C -> [(Double,Vector C)]+valsVecs h = let (valv,m) = eigSH h+ vals = H.toList valv+ vecs = map (conjV . fromList) $ toLists (conjugateTranspose m)+ in zip vals vecs++-- From a Hermitian matrix, a list of pairs of eigenvalues and projectors.+valsPs :: Herm C -> [(Double,Matrix C)]+valsPs h = [(val,couter vec vec) | (val,vec) <- valsVecs h]++combineFst :: (Eq a, Num b) => [(a,b)] -> [(a,b)]+combineFst [] = []+combineFst [p] = [p]+combineFst ((x1,m1):(x2,m2):ps)+ = if x1 == x2+ then combineFst ((x1,m1+m2):ps)+ else (x1,m1):combineFst ((x2,m2):ps)++-- | Given an obervable, return a list of pairs+-- of possible outcomes and projectors+-- for each outcome.+outcomesProjectors :: Matrix C -> [(Double,Matrix C)]+outcomesProjectors m = combineFst (valsPs (sym m))++-- | Given an observable and a state vector, return a list of pairs+-- of possible outcomes and probabilites+-- for each outcome.+outcomesProbabilities :: Matrix C -> Vector C -> [(Double,Double)]+outcomesProbabilities m v+ = [(a,realPart (inner v (p #> v))) | (a,p) <- outcomesProjectors m] ------------------ -- Gram-Schmidt --