learn-physics 0.5.2 → 0.6.7
raw patch · 27 files changed
Files
- LICENSE +1/−1
- examples/src/LorentzForceSimulation.hs +0/−3
- examples/src/NMR.hs +16/−0
- learn-physics.cabal +56/−38
- src/Physics/Learn.hs +1/−1
- src/Physics/Learn/BeamStack.hs +290/−0
- src/Physics/Learn/BlochSphere.hs +224/−0
- src/Physics/Learn/CarrotVec.hs +5/−3
- src/Physics/Learn/Charge.hs +2/−2
- src/Physics/Learn/CommonVec.hs +10/−7
- src/Physics/Learn/CompositeQuadrature.hs +2/−2
- src/Physics/Learn/CoordinateFields.hs +2/−2
- src/Physics/Learn/CoordinateSystem.hs +2/−2
- src/Physics/Learn/Current.hs +2/−2
- src/Physics/Learn/Curve.hs +2/−2
- src/Physics/Learn/Ket.hs +610/−0
- src/Physics/Learn/Mechanics.hs +2/−2
- src/Physics/Learn/Position.hs +2/−2
- src/Physics/Learn/QuantumMat.hs +405/−0
- src/Physics/Learn/RootFinding.hs +2/−2
- src/Physics/Learn/RungeKutta.hs +2/−2
- src/Physics/Learn/Schrodinger1D.hs +415/−0
- src/Physics/Learn/SimpleVec.hs +8/−22
- src/Physics/Learn/StateSpace.hs +2/−2
- src/Physics/Learn/Surface.hs +2/−2
- src/Physics/Learn/Visual/VisTools.hs +2/−2
- src/Physics/Learn/Volume.hs +2/−2
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2011-2015 Scott N. Walck <walck@lvc.edu>.+Copyright (c) 2011-2020 Scott N. Walck <walck@lvc.edu>. All rights reserved. Redistribution and use in source and binary forms, with or without
examples/src/LorentzForceSimulation.hs view
@@ -2,9 +2,6 @@ import Physics.Learn import Vis-import SpatialMath- ( Euler(..)- ) drawFunction :: SimpleState -> VisObject Double drawFunction (_t,r,_v)
+ examples/src/NMR.hs view
@@ -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)
learn-physics.cabal view
@@ -1,46 +1,52 @@ Name: learn-physics-Version: 0.5.2+Version: 0.6.7 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 Maintainer: Scott N. Walck <walck@lvc.edu> Category: Physics Build-type: Simple-Cabal-version: >=1.8-Tested-with: GHC == 7.10.2+Cabal-version: >=1.10 Library- Exposed-modules: Physics.Learn.Charge- Physics.Learn.Current- Physics.Learn.Position- Physics.Learn.Curve- Physics.Learn.Surface- Physics.Learn.Volume+ Exposed-modules: Physics.Learn+ Physics.Learn.BeamStack+ Physics.Learn.BlochSphere Physics.Learn.CarrotVec- Physics.Learn.SimpleVec+ Physics.Learn.Charge Physics.Learn.CommonVec+ Physics.Learn.CompositeQuadrature Physics.Learn.CoordinateFields Physics.Learn.CoordinateSystem- Physics.Learn.StateSpace- Physics.Learn.RungeKutta- Physics.Learn.CompositeQuadrature- Physics.Learn.RootFinding+ Physics.Learn.Current+ Physics.Learn.Curve+ Physics.Learn.Ket Physics.Learn.Mechanics- Physics.Learn+ Physics.Learn.Position+ Physics.Learn.QuantumMat+ Physics.Learn.RootFinding+ Physics.Learn.RungeKutta+ Physics.Learn.Schrodinger1D+ Physics.Learn.SimpleVec+ Physics.Learn.StateSpace+ Physics.Learn.Surface+ Physics.Learn.Visual.GlossTools Physics.Learn.Visual.PlotTools Physics.Learn.Visual.VisTools- Physics.Learn.Visual.GlossTools- Build-depends: base >= 4.2 && < 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+ Physics.Learn.Volume+ Build-depends: base >= 4.7 && < 5,+ gloss >= 1.8 && < 1.14,+ gnuplot >= 0.5 && < 0.6,+ hmatrix >= 0.17 && < 0.21,+ Vis >= 1.0.0 && < 1.1,+ linear >= 1.23 && < 1.24,+ vector-space >= 0.8.4 && < 0.17 Hs-source-dirs: src+ default-language: Haskell2010 Source-repository head type: git@@ -48,43 +54,55 @@ Executable learn-physics-PlaneWave Main-is: examples/src/PlaneWave.hs- Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ Build-depends: Vis >= 1.0.0 && < 1.1,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 Executable learn-physics-eFieldLine3D Main-is: examples/src/eFieldLine3D.hs- Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ Build-depends: Vis >= 1.0.0 && < 1.1,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 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,+ Build-depends: Vis >= 1.0.0 && < 1.1,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 Executable learn-physics-BCircularLoop Main-is: examples/src/BCircularLoop.hs- Build-depends: not-gloss >= 0.7.4 && < 0.8,- base >= 4.5 && < 4.9,+ Build-depends: Vis >= 1.0.0 && < 1.1,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 Executable learn-physics-sunEarth Main-is: examples/src/sunEarthRK4.hs- Build-depends: gloss >= 1.8 && < 1.10,- base >= 4.5 && < 4.9,+ Build-depends: gloss >= 1.8 && < 1.14,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 Executable learn-physics-eFieldLine2D Main-is: examples/src/eFieldLine2D.hs- Build-depends: gloss >= 1.8 && < 1.10,- base >= 4.5 && < 4.9,+ Build-depends: gloss >= 1.8 && < 1.14,+ base >= 4.5 && < 5, learn-physics+ default-language: Haskell2010 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 && < 5, learn-physics+ default-language: Haskell2010++Executable learn-physics-NMR+ Main-is: examples/src/NMR.hs+ Build-depends: base >= 4.5,+ learn-physics+ default-language: Haskell2010
src/Physics/Learn.hs view
@@ -3,7 +3,7 @@ {- | Module : Physics.Learn-Copyright : (c) Scott N. Walck 2014+Copyright : (c) Scott N. Walck 2014-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
+ src/Physics/Learn/BeamStack.hs view
@@ -0,0 +1,290 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE CPP #-}++{- | +Module : Physics.Learn.BeamStack+Copyright : (c) Scott N. Walck 2016-2018+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+ )+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif++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
+ src/Physics/Learn/BlochSphere.hs view
@@ -0,0 +1,224 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE CPP #-}++{- |+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+ , ketToPos+ , staticBlochSphere+ , displayStaticState+ , animatedBlochSphere+ , simulateBlochSphere+ , simulateBlochSphereK+ , stateProp+ , statePropK+ , evolutionBlochSphere+ , evolutionBlochSphereK+ , hamRabi+ )+ where++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+ , 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 Vis+ ( VisObject(..)+ , Flavour(..)+ , Options(..)+ , Camera0(..)+ , Euler(..)+ , defaultOpts+ , display+ , simulate+ , blue+ , red+ )+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif++{-+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++-- | 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))++-- | 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+ = 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++-- | 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 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.+-- 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, 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 M.xp (const (scale 0.9 M.sz))++staticBz2 :: IO ()+staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const M.sz)++staticBy1 :: IO ()+staticBy1 = evolutionBlochSphere M.xp (const M.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
src/Physics/Learn/CarrotVec.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_GHC -Wall -fno-warn-orphans #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.CarrotVec-Copyright : (c) Scott N. Walck 2011-2014+Copyright : (c) Scott N. Walck 2011-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental@@ -26,6 +26,7 @@ module Physics.Learn.CarrotVec ( Vec+ , R , xComp , yComp , zComp@@ -65,6 +66,7 @@ ) import Physics.Learn.CommonVec ( Vec(..)+ , R , xComp , yComp , zComp@@ -81,7 +83,7 @@ Vec ax ay az ^+^ Vec bx by bz = Vec (ax+bx) (ay+by) (az+bz) instance VectorSpace Vec where- type Scalar Vec = Double+ type Scalar Vec = R c *^ Vec ax ay az = Vec (c*ax) (c*ay) (c*az) instance InnerSpace Vec where
src/Physics/Learn/Charge.hs view
@@ -1,9 +1,9 @@ {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Charge-Copyright : (c) Scott N. Walck 2011-2014+Copyright : (c) Scott N. Walck 2011-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/CommonVec.hs view
@@ -3,7 +3,7 @@ {- | Module : Physics.Learn.CommonVec-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental@@ -23,6 +23,7 @@ module Physics.Learn.CommonVec ( Vec(..)+ , R , vec , (><) , iHat@@ -33,10 +34,12 @@ infixl 7 >< +type R = Double+ -- | A type for vectors.-data Vec = Vec { xComp :: Double -- ^ x component- , yComp :: Double -- ^ y component- , zComp :: Double -- ^ z component+data Vec = Vec { xComp :: R -- ^ x component+ , yComp :: R -- ^ y component+ , zComp :: R -- ^ z component } deriving (Eq) instance Show Vec where@@ -50,9 +53,9 @@ | otherwise = show x -- | Form a vector by giving its x, y, and z components.-vec :: Double -- ^ x component- -> Double -- ^ y component- -> Double -- ^ z component+vec :: R -- ^ x component+ -> R -- ^ y component+ -> R -- ^ z component -> Vec vec = Vec
src/Physics/Learn/CompositeQuadrature.hs view
@@ -1,10 +1,10 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.CompositeQuadrature-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/CoordinateFields.hs view
@@ -1,9 +1,9 @@ {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.CoordinateFields-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/CoordinateSystem.hs view
@@ -1,9 +1,9 @@ {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.CoordinateSystem-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/Current.hs view
@@ -1,9 +1,9 @@ {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Current-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/Curve.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE TypeFamilies, FlexibleContexts #-} {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Curve-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
+ src/Physics/Learn/Ket.hs view
@@ -0,0 +1,610 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE Safe #-}+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}+{-# LANGUAGE CPP #-}++{- | +Module : Physics.Learn.Ket+Copyright : (c) Scott N. Walck 2016-2018+License : BSD3 (see LICENSE)+Maintainer : Scott N. Walck <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+ (+ -- * Basic data types+ C+ , i+ , magnitude+ , Ket+ , Bra+ , Operator+ -- * Kets for spin-1/2 particles+ , xp+ , xm+ , yp+ , ym+ , zp+ , zm+ , np+ , nm+ -- * Operators for spin-1/2 particles+ , sx+ , sy+ , sz+ , sn+ , sn'+ -- * Quantum Dynamics+ , timeEvOp+ , timeEv+ -- * Composition+ , Kron(..)+ -- * Measurement+ , possibleOutcomes+ , outcomesProjectors+ , outcomesProbabilities+-- , prob+-- , probs+ -- * Generic multiplication+ , Mult(..)+ -- * Adjoint operation+ , Dagger(..)+ -- * Normalization+ , HasNorm(..)+ -- * Representation+ , Representable(..)+ -- * Orthonormal bases+ , OrthonormalBasis+ , makeOB+ , listBasis+ , size+ -- * Orthonormal bases for spin-1/2 particles+ , xBasis+ , yBasis+ , zBasis+ , nBasis+ -- , angularMomentumXMatrix+ -- , angularMomentumYMatrix+ -- , angularMomentumZMatrix+ -- , angularMomentumPlusMatrix+ -- , angularMomentumMinusMatrix+ -- , jXMatrix+ -- , jYMatrix+ -- , jZMatrix+ -- , matrixCommutator+ -- , rotationMatrix+ -- , jmColumn+ )+ where++-- We try to import only from QuantumMat+-- and not from Numeric.LinearAlgebra++import qualified Data.Complex as C+import Data.Complex+ ( Complex(..)+ , conjugate+ )+import qualified Physics.Learn.QuantumMat as M+import Physics.Learn.QuantumMat+ ( C+ , Vector+ , Matrix+ , (#>)+ , (<#)+ , conjugateTranspose+ , scaleV+ , scaleM+ , conjV+ , fromList+ , toList+ , fromLists+ )+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif++infixl 7 <>++-- | A ket vector describes the state of a quantum system.+data Ket = Ket (Vector C)++instance Show Ket where+ show 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 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)++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.+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+ (fromLists [[ x*y | y <- toList b] | x <- toList k])++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++-- | 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"++--------------+-- 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++-- | 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++-- | 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++-- | 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+-}++{-+----------------------------------------+-- 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
src/Physics/Learn/Mechanics.hs view
@@ -1,10 +1,10 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Mechanics-Copyright : (c) Scott N. Walck 2014+Copyright : (c) Scott N. Walck 2014-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/Position.hs view
@@ -1,10 +1,10 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Position-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
+ src/Physics/Learn/QuantumMat.hs view
@@ -0,0 +1,405 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE CPP #-}++{- | +Module : Physics.Learn.QuantumMat+Copyright : (c) Scott N. Walck 2016-2018+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+ , matrixFunction+ -- * Density matrices+ , couter+ , dm+ , trace+ , normalizeDM+ , oneQubitMixed+ -- * Quantum Dynamics+ , timeEvMat+ , timeEv+ , timeEvMatSpec+ -- * Composition+ , Kronecker(..)+ -- * Measurement+ , possibleOutcomes+ , outcomesProjectors+ , outcomesProbabilities+ -- * Vector and Matrix+ , Vector+ , Matrix+ )+ where++import Numeric.LinearAlgebra+ ( C+ , Vector+ , Matrix+ , Herm+ , iC -- square root of negative one+ , (><) -- matrix definition+ , ident+ , scale+ , norm_2+ , inv+ , (<\>)+ , sym+ , eigenvaluesSH+ , eigSH+ , 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+ , realPart+ )+#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif++-- | 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++--------------+-- Matrices --+--------------++-- | 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++-- | 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 --+----------------------++-- | 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 without 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)++-- | 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 --+-----------------++-- | 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)++-- 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 --+------------------++-- | 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+
src/Physics/Learn/RootFinding.hs view
@@ -3,7 +3,7 @@ {- | Module : Physics.Learn.RootFinding-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2017 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental@@ -26,7 +26,7 @@ -- | Given an initial bracketing of a root -- (an interval (a,b) for which f(a) f(b) <= 0),--- produce a bracket of arbitrary smallness.+-- produce a bracket (c,d) for which |c-d| < desired accuracy. bracketRoot :: (Ord a, Fractional a) => a -- ^ desired accuracy -> (a -> a) -- ^ function
src/Physics/Learn/RungeKutta.hs view
@@ -1,10 +1,10 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.RungeKutta-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
+ src/Physics/Learn/Schrodinger1D.hs view
@@ -0,0 +1,415 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE Trustworthy #-}++{- | +Module : Physics.Learn.Schrodinger1D+Copyright : (c) Scott N. Walck 2015-2018+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+ , xRange+ , listForm+ )+ where++import Data.Complex+ ( Complex(..)+ , magnitude+ )+import Graphics.Gloss+ ( Picture(..)+ , yellow+ , black+ , Display(..)+ , display+ )+-- import Math.Polynomial.Hermite+-- ( evalPhysHermite+-- )+import Numeric.LinearAlgebra+ ( R+ , C+ , Vector+ , Matrix+ , (|>)+ , (<.>)+ , fromLists+ , toList+ , size+ )+import Physics.Learn.QuantumMat+ ( probVector+ , timeEv+ )++--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+ :: R -- ^ length scale = sqrt(hbar / m omega)+ -> C -- ^ parameter z+ -> R -> C -- ^ wavefunction+coherent l z x+ = ((1/(pi*l**2))**0.25 * exp(-x**2/(2*l**2)) :+ 0)+ * exp(-z**2/2 + (sqrt(2/l**2) * x :+ 0) * z)++gaussian+ :: R -- ^ width parameter+ -> R -- ^ center of wave packet+ -> R -> C -- ^ wavefunction+gaussian a x0 x = exp(-(x-x0)**2/(2*a**2)) / sqrt(a * sqrt pi) :+ 0++movingGaussian+ :: R -- ^ width parameter+ -> R -- ^ center of wave packet+ -> R -- ^ l0 = hbar / p0+ -> R -> C -- ^ wavefunction+movingGaussian a x0 l0 x = exp((0 :+ x/l0) - ((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++-- options for representing wave functions+-- 1. A function R -> C+-- 2. ([R],Vector C), where lengths match+-- 3. [(R,C)]+-- 4. (R,R,Vector C) -- xmin, xmax, state vector (assumes even spacing)++-- 2,4 are best for evolution++listForm :: (R,R,Vector C) -> ([R],Vector C)+listForm (xmin,xmax,v)+ = let dt = (xmax - xmin) / fromIntegral (size v - 1)+ in ([xmin, xmin + dt .. xmax],v)+++{-+-- | 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 --+------------------++-- n is number of points+-- n-1 is number of intervals+xRange :: R -> R -> Int -> [R]+xRange xmin xmax n+ = let dt = (xmax - xmin) / fromIntegral (n - 1)+ in [xmin, xmin + dt .. xmax]+++{-+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? No.+-- 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)
src/Physics/Learn/SimpleVec.hs view
@@ -3,7 +3,7 @@ {- | Module : Physics.Learn.SimpleVec-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental@@ -17,24 +17,9 @@ easier for a person just learning Haskell. -} --- 2011 Apr 10--- Placed the code common to SimpleVec and CarrotVec in CommonVec---- 2011 Mar 19--- Add support for sumV, so that the interface matches CarrotVec.hs---- This uses the same internal data representation as SimpleVector,--- but uses an interface to match Conal Elliott's operators for--- vectors. (A similar interface to CarrotVector and SimpleCarrotVector.)--- The notation--- zeroV, negateV, (^+^), (^-^)--- is borrowed from Data.AdditiveGroup, and--- (*^), (^*), (^/), (<.>), magnitude--- is borrowed from Data.VectorSpace.--- Cross product operator is my own.- module Physics.Learn.SimpleVec ( Vec+ , R , xComp , yComp , zComp@@ -58,6 +43,7 @@ import Physics.Learn.CommonVec ( Vec(..)+ , R , vec , iHat , jHat@@ -95,23 +81,23 @@ -- | Scalar multiplication, where the scalar is on the left -- and the vector is on the right.-(*^) :: Double -> Vec -> Vec+(*^) :: R -> Vec -> Vec c *^ Vec ax ay az = Vec (c*ax) (c*ay) (c*az) -- | Scalar multiplication, where the scalar is on the right -- and the vector is on the left.-(^*) :: Vec -> Double -> Vec+(^*) :: Vec -> R -> Vec Vec ax ay az ^* c = Vec (c*ax) (c*ay) (c*az) -- | Division of a vector by a scalar.-(^/) :: Vec -> Double -> Vec+(^/) :: Vec -> R -> Vec Vec ax ay az ^/ c = Vec (ax/c) (ay/c) (az/c) -- | Dot product of two vectors.-(<.>) :: Vec -> Vec -> Double+(<.>) :: Vec -> Vec -> R Vec ax ay az <.> Vec bx by bz = ax*bx + ay*by + az*bz -- | Magnitude of a vector.-magnitude :: Vec -> Double+magnitude :: Vec -> R magnitude v = sqrt(v <.> v)
src/Physics/Learn/StateSpace.hs view
@@ -1,10 +1,10 @@ {-# OPTIONS_GHC -Wall -fno-warn-orphans #-} {-# LANGUAGE FlexibleContexts, FlexibleInstances, TypeFamilies #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.StateSpace-Copyright : (c) Scott N. Walck 2014+Copyright : (c) Scott N. Walck 2014-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/Surface.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Surface-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental
src/Physics/Learn/Visual/VisTools.hs view
@@ -12,13 +12,13 @@ ) where -import SpatialMath+import Linear ( V3(..)- , Euler(..) ) import Vis ( VisObject(..) , Color+ , Euler(..) ) import Physics.Learn.CarrotVec ( Vec
src/Physics/Learn/Volume.hs view
@@ -1,10 +1,10 @@ {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Safe #-} {- | Module : Physics.Learn.Volume-Copyright : (c) Scott N. Walck 2012-2014+Copyright : (c) Scott N. Walck 2012-2019 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck <walck@lvc.edu> Stability : experimental