diff --git a/examples/src/BCircularLoop.hs b/examples/src/BCircularLoop.hs new file mode 100644 --- /dev/null +++ b/examples/src/BCircularLoop.hs @@ -0,0 +1,30 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Physics.Learn +import Vis + +loopCurve :: Curve +loopCurve = Curve (\phi -> cyl 1 phi 0) 0 (2*pi) + +loop :: CurrentDistribution +loop = LineCurrent 20 loopCurve + +samplePoints :: [Position] +samplePoints = [cyl s phi z | + s <- [0.25,0.75..1.75] + , phi <- [pi/6,pi/2..2*pi] + , z <- [-1.5,-1..1.5]] + +arrows :: VisObject Double +arrows = displayVectorField blue 5e-5 samplePoints (bField loop) + +drawFun :: VisObject Double +drawFun = VisObjects [curveObject red loopCurve, arrows] + +myOptions :: Options +myOptions = defaultOpts {optWindowName = "Magnetic Field from a Current Loop"} + +main :: IO () +main = display myOptions drawFun diff --git a/examples/src/LorentzForceSimulation.hs b/examples/src/LorentzForceSimulation.hs new file mode 100644 --- /dev/null +++ b/examples/src/LorentzForceSimulation.hs @@ -0,0 +1,64 @@ +module Main where + +import Physics.Learn +import Vis +import SpatialMath + ( Euler(..) + ) + +drawFunction :: SimpleState -> VisObject Double +drawFunction (_t,r,_v) + = RotEulerDeg (Euler 270 0 0) $ RotEulerDeg (Euler 0 180 0) $ + VisObjects [ Axes (0.5, 15) + , Trans (v3FromPos r) (Sphere 0.1 Solid red) + ] + +statePropagationFunction :: Float -> SimpleState -> SimpleState +statePropagationFunction t' (t,r,v) = rungeKutta4 newton2 (realToFrac t' - t) (t,r,v) + +-- Newton's Second Law +newton2 :: SimpleState -> Diff SimpleState +newton2 (t,r,v) = (1,v,force (t,r,v) ^/ m) + +-- Lorentz Force Law +force :: SimpleState -> Vec +force (_t,r,v) = q *^ (electricField r ^+^ v >< magneticField r) + +myOptions :: Options +myOptions = defaultOpts {optWindowName = "Particle Experiencing Electromagnetic Force"} + +main :: IO () +main = simulate + myOptions + 0.01 + (0,initialPosition,initialVelocity) + drawFunction + statePropagationFunction + +-- particle mass +m :: Double +m = 1 + +-- particle charge +q :: Double +q = 1 + +-- Electric Field +electricField :: VectorField +electricField r = vec 0 2 0 + where + (x,y,z) = cartesianCoordinates r + +-- Magnetic Field +magneticField :: VectorField +magneticField r = vec 0 0 4 + where + (x,y,z) = cartesianCoordinates r + +-- Initial displacement +initialPosition :: Position +initialPosition = cart 0 0 0 + +-- Initial velocity +initialVelocity :: Vec +initialVelocity = vec 0 0 0 diff --git a/examples/src/NMR.hs b/examples/src/NMR.hs new file mode 100644 --- /dev/null +++ b/examples/src/NMR.hs @@ -0,0 +1,16 @@ +{-# OPTIONS_GHC -Wall #-} + +-- ^ Nuclear Magnetic Resonance on the Bloch Sphere + +module Main where + +import Physics.Learn.QuantumMat + ( zm + ) +import Physics.Learn.BlochSphere + ( hamRabi + , evolutionBlochSphere + ) + +main :: IO () +main = evolutionBlochSphere zm (hamRabi 10 1 10) diff --git a/examples/src/PlaneWave.hs b/examples/src/PlaneWave.hs new file mode 100644 --- /dev/null +++ b/examples/src/PlaneWave.hs @@ -0,0 +1,48 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Vis + ( animate + , VisObject(..) + , red + , blue + , Options(..) + , defaultOpts + ) +import Physics.Learn.CarrotVec + ( vec + ) +import Physics.Learn + ( Position + , VectorField + , displayVectorField + , cart + , cartesianCoordinates + ) + +samplePoints :: [Position] +samplePoints = [cart x y z | x <- [-2,0,2], y <- [-2,0,2], z <- [-4,-3.6..4]] + +drawFun :: Float -> VisObject Double +drawFun time = VisObjects [displayVectorField blue 1 samplePoints (eField t) + ,displayVectorField red 1 samplePoints (bField t) + ] + where + t = realToFrac time + +eField :: Double -> VectorField +eField t r = vec (cos (z - t)) 0 0 + where + (_,_,z) = cartesianCoordinates r + +bField :: Double -> VectorField +bField t r = vec 0 (cos (z - t)) 0 + where + (_,_,z) = cartesianCoordinates r + +myOptions :: Options +myOptions = defaultOpts {optWindowName = "Plane Wave"} + +main :: IO () +main = animate myOptions drawFun diff --git a/examples/src/eFieldLine3D.hs b/examples/src/eFieldLine3D.hs new file mode 100644 --- /dev/null +++ b/examples/src/eFieldLine3D.hs @@ -0,0 +1,48 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Vis + ( display + , VisObject(..) + , red + , blue + , Options(..) + , defaultOpts + ) +import Physics.Learn.Visual.VisTools + ( curveObject + , displayVectorField + ) +import Physics.Learn.Position + ( Position + , cart + ) +import Physics.Learn.Curve + ( Curve(..) + ) +import Physics.Learn.Charge + ( ChargeDistribution(..) + , eField + ) + +curve1 :: Curve +curve1 = Curve (\t -> cart t 0 0) (-4) 4 + +lineCharge :: ChargeDistribution +lineCharge = LineCharge (const 1e-9) curve1 + +samplePoints :: [Position] +samplePoints = [cart x y z | x <- [-8,-6..8], y <- [-4,-2..4], z <- [-4,-2..4], abs y + abs z > 0.5 || abs x > 4.5] + +arrows :: VisObject Double +arrows = displayVectorField blue 10 samplePoints (eField lineCharge) + +drawFun :: VisObject Double +drawFun = VisObjects [curveObject red curve1, arrows] + +myOptions :: Options +myOptions = defaultOpts {optWindowName = "Electric Field from a Line Charge"} + +main :: IO () +main = display myOptions drawFun diff --git a/learn-physics.cabal b/learn-physics.cabal --- a/learn-physics.cabal +++ b/learn-physics.cabal @@ -1,5 +1,5 @@ Name: learn-physics -Version: 0.6.2 +Version: 0.6.3 Synopsis: Haskell code for learning physics Description: A library of functions for vector calculus, calculation of electric field, electric flux, @@ -13,67 +13,68 @@ Build-type: Simple Cabal-version: >=1.8 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.Visual.PlotTools - Physics.Learn.Visual.GlossTools - Physics.Learn + Physics.Learn.Position Physics.Learn.QuantumMat - Physics.Learn.Ket - Physics.Learn.BeamStack --- Physics.Learn.BlochSphere --- Physics.Learn.Visual.VisTools + 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.Volume Build-depends: base >= 4.7 && < 4.12, vector-space >= 0.8.4, hmatrix >= 0.17, gloss >= 1.8, - gnuplot >= 0.5 && < 0.6 --- not-gloss >= 0.5.0.4, --- spatial-math >= 0.2, + gnuplot >= 0.5 && < 0.6, + not-gloss >= 0.5.0.4, + spatial-math >= 0.2 Hs-source-dirs: src Source-repository head type: git location: https://github.com/walck/learn-physics --- Executable learn-physics-PlaneWave --- Main-is: examples/src/PlaneWave.hs --- Build-depends: not-gloss >= 0.7.4, --- base >= 4.5 && < 4.12, --- learn-physics +Executable learn-physics-PlaneWave + Main-is: examples/src/PlaneWave.hs + Build-depends: not-gloss >= 0.7.4, + base >= 4.5 && < 4.12, + learn-physics --- Executable learn-physics-eFieldLine3D --- Main-is: examples/src/eFieldLine3D.hs --- Build-depends: not-gloss >= 0.7.4, --- base >= 4.5 && < 4.12, --- learn-physics +Executable learn-physics-eFieldLine3D + Main-is: examples/src/eFieldLine3D.hs + Build-depends: not-gloss >= 0.7.4, + base >= 4.5 && < 4.12, + learn-physics --- Executable learn-physics-LorentzForceSimulation --- Main-is: examples/src/LorentzForceSimulation.hs --- Build-depends: not-gloss >= 0.7.4, --- spatial-math >= 0.2, --- base >= 4.5 && < 4.12, --- learn-physics +Executable learn-physics-LorentzForceSimulation + Main-is: examples/src/LorentzForceSimulation.hs + Build-depends: not-gloss >= 0.7.4, + spatial-math >= 0.2, + base >= 4.5 && < 4.12, + learn-physics --- Executable learn-physics-BCircularLoop --- Main-is: examples/src/BCircularLoop.hs --- Build-depends: not-gloss >= 0.7.4, --- base >= 4.5 && < 4.12, --- learn-physics +Executable learn-physics-BCircularLoop + Main-is: examples/src/BCircularLoop.hs + Build-depends: not-gloss >= 0.7.4, + base >= 4.5 && < 4.12, + learn-physics Executable learn-physics-sunEarth Main-is: examples/src/sunEarthRK4.hs @@ -93,7 +94,7 @@ base >= 4.5 && < 4.12, learn-physics --- Executable learn-physics-NMR --- Main-is: examples/src/NMR.hs --- Build-depends: base >= 4.5, --- learn-physics +Executable learn-physics-NMR + Main-is: examples/src/NMR.hs + Build-depends: base >= 4.5, + learn-physics diff --git a/src/Physics/Learn.hs b/src/Physics/Learn.hs --- a/src/Physics/Learn.hs +++ b/src/Physics/Learn.hs @@ -166,12 +166,12 @@ , arrow , thickArrow -- ** Vis library - -- , v3FromVec - -- , v3FromPos - -- , visVec - -- , oneVector - -- , displayVectorField - -- , curveObject + , v3FromVec + , v3FromPos + , visVec + , oneVector + , displayVectorField + , curveObject ) where @@ -275,14 +275,14 @@ , shiftVolume , volumeIntegral ) --- import Physics.Learn.Visual.VisTools --- ( v3FromVec --- , v3FromPos --- , visVec --- , oneVector --- , displayVectorField --- , curveObject --- ) +import Physics.Learn.Visual.VisTools + ( v3FromVec + , v3FromPos + , visVec + , oneVector + , displayVectorField + , curveObject + ) import Physics.Learn.StateSpace ( StateSpace(..) , (.-^) diff --git a/src/Physics/Learn/BlochSphere.hs b/src/Physics/Learn/BlochSphere.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/BlochSphere.hs @@ -0,0 +1,226 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE CPP #-} + +{- | +Module : Physics.Learn.BlochSphere +Copyright : (c) Scott N. Walck 2016 +License : BSD3 (see LICENSE) +Maintainer : Scott N. Walck +Stability : experimental + +This module contains functions for displaying the +state of a spin-1/2 particle or other quantum two-level +system as a point on the Bloch Sphere. +-} + +module Physics.Learn.BlochSphere + ( VisObj + , toPos + , 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 SpatialMath + ( Euler(..) + ) +import Vis + ( VisObject(..) + , Flavour(..) + , Options(..) + , Camera0(..) + , 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 diff --git a/src/Physics/Learn/Schrodinger1D.hs b/src/Physics/Learn/Schrodinger1D.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Schrodinger1D.hs @@ -0,0 +1,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 +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 -- ^ , 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) diff --git a/src/Physics/Learn/SimpleVec.hs b/src/Physics/Learn/SimpleVec.hs --- a/src/Physics/Learn/SimpleVec.hs +++ b/src/Physics/Learn/SimpleVec.hs @@ -3,7 +3,7 @@ {- | Module : Physics.Learn.SimpleVec -Copyright : (c) Scott N. Walck 2012-2014 +Copyright : (c) Scott N. Walck 2012-2018 License : BSD3 (see LICENSE) Maintainer : Scott N. Walck 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 @@ -72,6 +57,8 @@ infixl 7 ^/ infixl 7 <.> +type R = Double + -- | The zero vector. zeroV :: Vec zeroV = vec 0 0 0 @@ -95,23 +82,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) diff --git a/src/Physics/Learn/Visual/VisTools.hs b/src/Physics/Learn/Visual/VisTools.hs new file mode 100644 --- /dev/null +++ b/src/Physics/Learn/Visual/VisTools.hs @@ -0,0 +1,126 @@ +{-# OPTIONS_GHC -Wall #-} + +-- | Some tools related to the not-gloss 3D graphics and animation library. + +module Physics.Learn.Visual.VisTools + ( v3FromVec + , v3FromPos + , visVec + , oneVector + , displayVectorField + , curveObject + ) + where + +import SpatialMath + ( V3(..) + , Euler(..) + ) +import Vis + ( VisObject(..) + , Color + ) +import Physics.Learn.CarrotVec + ( Vec + , xComp + , yComp + , zComp +-- , magnitude + , (^/) + ) +import Physics.Learn.Position + ( Position + , cartesianCoordinates + , VectorField + ) +import Physics.Learn.Curve + ( Curve(..) + ) + +-- | Make a 'V3' object from a 'Vec'. +v3FromVec :: Vec -> V3 Double +v3FromVec v = V3 x y z + where + x = xComp v + y = yComp v + z = zComp v + +-- | Make a 'V3' object from a 'Position'. +v3FromPos :: Position -> V3 Double +v3FromPos r = V3 x y z + where + (x,y,z) = cartesianCoordinates r + +-- | Display a vector field. +displayVectorField :: Color -- ^ color for the vector field + -> Double -- ^ scale factor + -> [Position] -- ^ list of positions to show the field + -> VectorField -- ^ vector field to display + -> VisObject Double -- ^ the displayable object +displayVectorField col unitsPerMeter samplePts field + = VisObjects [Trans (v3FromPos r) $ visVec col (e ^/ unitsPerMeter) | r <- samplePts, let e = field r] + +-- | A displayable VisObject for a curve. +curveObject :: Color -> Curve -> VisObject Double +curveObject color (Curve f a b) + = Line' Nothing [(v3FromPos (f t), color) | t <- [a,a+(b-a)/1000..b]] + +-- | Place a vector at a particular position. +oneVector :: Color -> Position -> Vec -> VisObject Double +oneVector c r v = Trans (v3FromPos r) $ visVec c v + +data Cart = Cart Double Double Double + deriving (Show) + +data Sph = Sph Double Double Double + deriving (Show) + +sphericalCoords :: Cart -> Sph +sphericalCoords (Cart x y z) = Sph r theta phi + where + r = sqrt (x*x + y*y + z*z) + s = sqrt (x*x + y*y) + theta = atan2 s z + phi = atan2 y x + +-- | A VisObject arrow from a vector +visVec :: Color -> Vec -> VisObject Double +visVec c v = rotZ phi $ rotY theta $ Arrow (r,20*r) (V3 0 0 1) c + where + x = xComp v + y = yComp v + z = zComp v + Sph r theta phi = sphericalCoords (Cart x y z) + +{- +rotX :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotX alpha = RotEulerRad (Euler 0 0 alpha) +-} + +rotY :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotY alpha = RotEulerRad (Euler 0 alpha 0) + +rotZ :: Double -- ^ in radians + -> VisObject Double + -> VisObject Double +rotZ alpha = RotEulerRad (Euler alpha 0 0) + + +{- +adjacentDistance :: [Position] -> Double +adjacentDistance [] = 0 +adjacentDistance rs'@(_:rs) = minimum (map magnitude $ zipWith displacement rs' rs) + +visVectorField :: Color -> [Position] -> VectorField -> VisObject Double +visVectorField c rs vf = let prs = [(r,vf r) | r <- rs] + bigV = maximum [magnitude (snd pr) | pr <- prs] + disp = adjacentDistance rs + scaleFactor = disp / bigV + newPrs = [(r, scaleFactor *^ v) | (r,v) <- prs] + vecs = [oneVector c r v' | (r,v') <- newPrs] + in VisObjects vecs +-}