packages feed

learn-physics 0.6.2 → 0.6.3

raw patch · 11 files changed

+1044/−83 lines, 11 filesdep +not-glossdep +spatial-mathdep ~basenew-component:exe:learn-physics-BCircularLoopnew-component:exe:learn-physics-LorentzForceSimulationnew-component:exe:learn-physics-NMRnew-component:exe:learn-physics-PlaneWavenew-component:exe:learn-physics-eFieldLine3D

Dependencies added: not-gloss, spatial-math

Dependency ranges changed: base

Files

+ examples/src/BCircularLoop.hs view
@@ -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
+ examples/src/LorentzForceSimulation.hs view
@@ -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
+ 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)
+ examples/src/PlaneWave.hs view
@@ -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
+ examples/src/eFieldLine3D.hs view
@@ -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
learn-physics.cabal view
@@ -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
src/Physics/Learn.hs view
@@ -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(..)     , (.-^)
+ src/Physics/Learn/BlochSphere.hs view
@@ -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 <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 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
+ 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-2018 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@@ -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)
+ src/Physics/Learn/Visual/VisTools.hs view
@@ -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+-}