|
|
２次元-電子動力学法の雛形プログラム(005wavePacketQED2D.bas)を公開します。
ポテンシャル中の電子の波動関数psi(x,y,t)とポテンシャルV(x,y)からdt時間後のpsi(x,y,t+dt)を求めます。
これを繰り返して時間発展させていきます。この波動関数から、電子の存在確率やエネルギーなどが解ります。
２次元になると、複素関数の表現方法が難しく、ここでは存在確率<psi|psi>しか表現できていません。
何か良い方法があれば、教えてください。主プログラム中のCALL drawWave(2)をCALL drawWave(4)に変更すると、
位相が表現されていますが、成功しているとは思えません。
この雛形プログラムでは、放物型ポテンシャル中の電子を追跡します。古典的な放物型ポテンシャル中のボールの動きと
類似しているように見えます。全エネルギーは一定になるはずですが、わずかに変動が見られます。
１次元と比較すると、およそ100倍ほど計算量が多くなるため、ほとんど動いて見えません。なんとか速くしたいものです。
表示の説明：
緑色は放物型ポテンシャルを、オリーブ色は電子の存在確率を表します。
試験環境：
本プログラムは十進BASIC 6.6.0/macOS 10.7.5と 十進BASIC Ver 7.7.8/windows 10でテストしました。
!
! ========= Quantum Electron Dynamics 2D ==========
!
! 005wavePacketQED2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.02.19 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB qed2d.setInitialCondition, qed2d.evolveWave, qed2d.drawWave
CALL setInitialCondition
FOR it=1 TO 1000
CALL evolveWave(0) !0:no loss, 1: loss energy (steepest-descent)
CALL drawWave(2) !0:psi^2, 1:cloud, 2:grid3D, 3:prob current, 5:phase
NEXT it
END
! ---------- QED2D module (QED: Quantum Electron Dynamics) ----------
!
! - time dependent Schrodinger equation: i(d/dt)psi(r,t) = H psi(r,t)
! - time evolution
! psi(r,t+dt) = exp(-i dt H) psi(r,t), (H:Hamiltonian of the system)
! H = -delta/2 + V(r), delta = d^2/dx^2 + d^2/dy^2 + d^2/dz^2
! psi(r,t+dt) = exp(-i dt H) psi(r,t) nearly=
! {exp(-i(dt/2)V} {exp(i dt(delta/2)} {exp(-i(dt/2)V} psi(r,t)
! - algorism: {exp(i dt(delta/2)}
! QED: Watanabe's algorism (semi-implicit method)
! Naoki Watanabe, Masaru Tsukada; arXiv:physics/0011068v1
! (Published from Physical Review E. 62, 2914, (2000).)
!
! Cayley's form : exp(i dt delta/2) nearly= (1 + i dt delta/4)/(1 - i dt delta/4)
! psi(r,t+dt) = exp(i dt delta/2) psi(r,dt)
! (1 - i dt delta/4) psi(r,t+dt) = (1 + i dt delta/4) psi(r,t)
!
! difference form psi(r,t) --> psi(j,n)
! psi(j,n+1) - i (dt/dx^2)/4 {psi(j-1,n+1))-2psi(j,n+1)+psi(j+1,n+1)}
! = psi(j,n) + i (dt/dx^2)/4 {psi(j-1,n))-2psi(j,n)+psi(j+1,n)}
! x i(4dx^2/dt) by each term
! psi(j-1,n+1) + A Psi(j,n+1) + psi(j+1,n+1) = -psi(j-1,n) + B Psi(j,n) -psi(j+1,n)
! where A=(i4dx^2/dt)-2, B=(i4dx^2/dt)+2
! bnj = -psi(j-1,n) + B Psi(j,n) -psi(j+1,n) is calculated using known psi(j,n)
! psi(j-1,n+1) + A Psi(j,n+1) + psi(j+1,n+1) = bnj
!
! solve tri-diagonal equation A X = B
! | a1 1 0 0 | | x1 | | b1 |
! | 1 a2 1 0 | | x2 | = | b2 |
! | 0 1 a3 1 | | x3 | | b3 |
! | 0 0 1 a4 | | x4 | | b4 |
!
! u(1) = 1.0/a(1) ! u() : work vector
! x(1) = b(1)*u(1)
!
! FOR i=2 TO N-2 ! forward elimination
! u(i) = 1/(a(i)-u(i-1))
! x(i) = (b(i)-x(i-1))*u(i)
! NEXT i
!
! FOR i=N-3 TO 1 STEP -1 ! backward substitution
! x(i) -= x(i+1)*u(i)
! NEXT i
!
MODULE qed2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, evolveWave, drawWave
PUBLIC FUNCTION systemTime, norm
SHARE NUMERIC sysTime, NNx, NNy, dx, dy, dt, auLength, auTime, auEnergy
SHARE NUMERIC cosAx, sinAx, cosAy, sinAy, cx0, cy0, cz0 !--- 3D graphics
SHARE NUMERIC psi(320,320,0 TO 1) ! wave function (,,0):real part (,,1):imaginary part
SHARE NUMERIC vv(320,320) ! external potential
SHARE NUMERIC wrk(320,320,0 TO 1) ! work space in lossEnergy
SHARE NUMERIC bRe(400),bIm(400) ! b vector in kxStep,kyStep
SHARE NUMERIC uRe(400),uIm(400) ! u vector in kxStep,kyStep
SHARE NUMERIC psicol(160,160) ! MAT PLOT CELL matrix for psi(x,y)
SHARE NUMERIC vcol(160,160) ! MAT PLOT CELL matrix for potential V(x,y)
SHARE NUMERIC srnd(1000) ! 1000 RND orderd series 0 to 1,use drawCloud
LET sysTime = 0.0 ! (au) au : atomic unit (hbar=1, me=1, e=1)
LET NNx = 160 ! max number of psi(x,,)
LET NNy = 160 ! max number of psi(,y,)
LET dx = 0.5 ! (au) x division
LET dy = 0.5 ! (au) y division
LET dt = 1.0*dx*dx ! (au) time division dt/(dx*dx)>3 ~ unstable?
LET auLength = 5.29177e-11 ! (m) 1(au) = auLength (m)
LET auTime = 2.41888e-17 ! (s) 1(au) = auTime (s)
LET auEnergy = 4.38975e-18 ! (J) 1(au) = auEnergy (J) (= 27.2114 eV)
! ---------- set initial condition
EXTERNAL SUB setInitialCondition !public
DECLARE EXTERNAL SUB setGaussianWave,setHarmonicPotential,initDraw
RANDOMIZE
LET sysTime = 0.0
CALL setGaussianWave(NNx*dx/2,NNy*dy/4,3,1,0) !(xPos,yPos,waveWidth,kx,ky)
CALL setHarmonicPotential(2.0) ! v(x)=k*(x-x0)^2
CALL initDraw !set color pallet and vcol(,)
! set window
LET xMargin = 60
LET yMargin = 120
SET WINDOW -xMargin,500-xMargin,-yMargin,500-yMargin
END SUB
EXTERNAL SUB setGaussianWave(xPos,yPos,waveWidth,kx,ky)
DECLARE EXTERNAL SUB normalize
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET x = i*dx
LET y = j*dy
LET phAb = EXP(-((x-xPos)*(x-xPos)+(y-yPos)*(y-yPos))/(4*waveWidth*waveWidth) )
LET phPh = kx*x+ky*y
LET psi(i,j,0) = phAb*COS(phPh)
LET psi(i,j,1) = phAb*SIN(phPh)
NEXT j
NEXT i
FOR i=1 TO NNx-1
LET psi(i,0,0) = 0
LET psi(i,0,1) = 0
LET psi(i,NNy-1,0) = 0
LET psi(i,NNy-1,1) = 0
NEXT i
FOR j=0 TO NNy-1
LET psi(0,j,0) = 0
LET psi(0,j,1) = 0
LET psi(NNx-1,j,0) = 0
LET psi(NNx-1,j,1) = 0
NEXT j
CALL normalize
END SUB
EXTERNAL SUB setHarmonicPotential(k0) !--- V(r)= k0*r^2
LET aa = k0/(NNx*dx*NNx*dx/4)
LET x0 = NNx*dx/2.0
LET y0 = NNy*dy/2.0
FOR i=0 TO NNx-1
FOR j=0 TO NNy-1
LET x = i*dx
LET y = j*dy
LET vv(i,j) = aa*((x-x0)*(x-x0)+(y-y0)*(y-y0))
NEXT j
NEXT i
END SUB
EXTERNAL SUB initDraw !--- set set color pallet and MAT PLOT matrix vcol(,)
FOR i = 0 TO 100 !--- set color pallet
SET COLOR MIX(40+i) 0.01*i,0.01*i,0 ! yellow for <psi|psi>
SET COLOR MIX(150+i) 0,0.01*i,0 ! green for V(x)
NEXT i
FOR i=0 TO NNx-1 !--- set vcol(,)
FOR j=0 TO NNy-1
LET col = 0.25*vv(i,j)
IF col>1 THEN LET col = 1
IF col<0 THEN LET col = 0
LET vcol(i,j) = 150+INT(col*100)
NEXT j
NEXT i
END SUB
! ---------- evolve Wave
EXTERNAL SUB evolveWave(lossSW) !public
DECLARE EXTERNAL SUB phaseStep,kxStep,kyStep,lossEnergy
CALL phaseStep
CALL kxStep
CALL kyStep
CALL phaseStep
LET sysTime = sysTime + dt
IF (lossSW=1) THEN CALL lossEnergy(0.01)
END SUB
EXTERNAL SUB phaseStep !--- {exp(-i(dt/2)V} : evolve 0.5dt
FOR i=1 TO NNx-2
FOR j=1 TO NNx-2
LET th = 0.5*dt*vv(i,j)
LET cs = COS(th)
LET sn = SIN(th)
LET phr = psi(i,j,0)
LET phi = psi(i,j,1)
LET psi(i,j,0) = cs*phr+sn*phi
LET psi(i,j,1) = -sn*phr+cs*phi
NEXT j
NEXT i
END SUB
EXTERNAL SUB kxStep !--- {exp(i dt(delta/2)} semi-implicit method : evolve dt
LET a = 4.0*dy*dy/dt
LET aaAb = 4.0+a*a
FOR j=1 TO NNy-2
FOR i=1 TO NNx-2
LET bRe(i) = 2*psi(i,j,0)-a*psi(i,j,1) - psi(i+1,j,0) - psi(i-1,j,0)
LET bIm(i) = 2*psi(i,j,1)+a*psi(i,j,0) - psi(i+1,j,1) - psi(i-1,j,1)
NEXT i
LET uRe(1) = -2/aaAb
LET uIm(1) = -a/aaAb
LET psi(1,j,0) = bRe(1)*uRe(1) - bIm(1)*uIm(1)
LET psi(1,j,1) = bIm(1)*uRe(1) + bRe(1)*uIm(1)
FOR i=2 TO NNx-2 !--- forward elimination
LET auAb = (-2-uRe(i-1))*(-2-uRe(i-1))+(a-uIm(i-1))*(a-uIm(i-1))
LET uRe(i) = (-2-uRe(i-1))/auAb
LET uIm(i) = -(a-uIm(i-1))/auAb
LET psi(i,j,0) = (bRe(i)-psi(i-1,j,0))*uRe(i) - (bIm(i)-psi(i-1,j,1))*uIm(i)
LET psi(i,j,1) = (bRe(i)-psi(i-1,j,0))*uIm(i) + (bIm(i)-psi(i-1,j,1))*uRe(i)
NEXT i
FOR i=NNx-3 TO 1 STEP -1 !--- backward substitution
LET psi(i,j,0) = psi(i,j,0) - (psi(i+1,j,0)*uRe(i) - psi(i+1,j,1)*uIm(i))
LET psi(i,j,1) = psi(i,j,1) - (psi(i+1,j,0)*uIm(i) + psi(i+1,j,1)*uRe(i))
NEXT i
NEXT j
END SUB
EXTERNAL SUB kyStep !--- {exp(i dt(delta/2)} semi-implicit method : evolve dt
LET a = 4*dy*dy/dt
LET aaAb = 4+a*a
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET bRe(j) = 2*psi(i,j,0)-a*psi(i,j,1) - psi(i,j+1,0) - psi(i,j-1,0)
LET bIm(j) = 2*psi(i,j,1)+a*psi(i,j,0) - psi(i,j+1,1) - psi(i,j-1,1)
NEXT j
LET uRe(1) = -2/aaAb
LET uIm(1) = -a/aaAb
LET psi(i,1,0) = bRe(1)*uRe(1) - bIm(1)*uIm(1)
LET psi(i,1,1) = bIm(1)*uRe(1) + bRe(1)*uIm(1)
FOR j=2 TO NNy-2 !--- forward elimination
LET auAb = (-2-uRe(j-1))*(-2-uRe(j-1))+(a-uIm(j-1))*(a-uIm(j-1))
LET uRe(j) = (-2-uRe(j-1))/auAb
LET uIm(j) = -(a-uIm(j-1))/auAb
LET psi(i,j,0) = (bRe(j)-psi(i,j-1,0))*uRe(j) - (bIm(j)-psi(i,j-1,1))*uIm(j)
LET psi(i,j,1) = (bRe(j)-psi(i,j-1,0))*uIm(j) + (bIm(j)-psi(i,j-1,1))*uRe(j)
NEXT j
FOR j=NNy-3 TO 1 STEP -1 !--- backward substitution
LET psi(i,j,0) = psi(i,j,0) - (psi(i,j+1,0)*uRe(j) - psi(i,j+1,1)*uIm(j))
LET psi(i,j,1) = psi(i,j,1) - (psi(i,j+1,0)*uIm(j) + psi(i,j+1,1)*uRe(j))
NEXT j
NEXT i
END SUB
EXTERNAL SUB lossEnergy(dumpingFactor) ! steepest-descent
DECLARE EXTERNAL SUB normalize
LET h2 = 2.0*dx*dx
LET ee = kineticEnergy + potentialEnergy
FOR i=1 TO NNx-2
FOR j=1 TO NNx-2
LET wrk(i,j,0) = -(psi(i+1,j,0)+psi(i-1,j,0)+psi(i,j+1,0)+psi(i,j-1,0)-4*psi(i,j,0))/h2+(vv(i,j)-ee)*psi(i,j,0)
LET wrk(i,j,1) = -(psi(i+1,j,1)+psi(i-1,j,1)+psi(i,j+1,1)+psi(i,j-1,1)-4*psi(i,j,1))/h2+(vv(i,j)-ee)*psi(i,j,1)
NEXT j
NEXT i
FOR i=1 TO NNx-2
FOR j=1 TO NNx-2
LET psi(i,j,0) = psi(i,j,0) - dumpingFactor*wrk(i,j,0)
LET psi(i,j,1) = psi(i,j,1) - dumpingFactor*wrk(i,j,1)
NEXT j
NEXT i
CALL normalize
END SUB
! ---------- utility
EXTERNAL FUNCTION systemTime !public
LET systemTime = sysTime
END FUNCTION
EXTERNAL FUNCTION norm !public <psi|psi>
LET p = 0
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET p = p + (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))
NEXT j
NEXT i
LET norm = p*dx*dy
END FUNCTION
EXTERNAL SUB normalize
DECLARE EXTERNAL FUNCTION norm
LET a = 1/SQR(norm)
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET psi(i,j,0) = psi(i,j,0)*a
LET psi(i,j,1) = psi(i,j,1)*a
NEXT j
NEXT i
END SUB
EXTERNAL FUNCTION potentialEnergy !--- <psi| V(x) |psi>
LET p = 0
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET p = p + (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*vv(i,j)
NEXT j
NEXT i
LET potentialEnergy = p*dx*dy
END FUNCTION
EXTERNAL FUNCTION kineticEnergy !--- <psi| -0.5(d^2/dx^2+d^2/dy^2) |psi>
LET h2 = dx*dx
LET p = 0
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET d2phRe = (psi(i+1,j,0)+psi(i-1,j,0)+psi(i,j+1,0)+psi(i,j-1,0)-4.0*psi(i,j,0))/h2
LET d2phIm = (psi(i+1,j,1)+psi(i-1,j,1)+psi(i,j+1,1)+psi(i,j-1,1)-4.0*psi(i,j,1))/h2
LET p = p + (psi(i,j,0)*d2phRe+psi(i,j,1)*d2phIm)
NEXT j
NEXT i
LET kineticEnergy = -0.5*p*dx*dy
END FUNCTION
EXTERNAL FUNCTION xProbabilityCurrentDensity(i,j) !--- Re(<psi| -i(d/dx) |psi>)
LET pRe = (psi(i+1,j,1)-psi(i-1,j,1))/(2*dx)
LET pIm = (-psi(i+1,j,0)+psi(i-1,j,0))/(2*dx)
LET xProbabilityCurrentDensity = psi(i,j,0)*pRe + psi(i,j,1)*pIm
END FUNCTION
EXTERNAL FUNCTION yProbabilityCurrentDensity(i,j) !--- Re(<psi| -i(d/dy) |psi>)
LET pRe = (psi(i,j+1,1)-psi(i,j-1,1))/(2*dy)
LET pIm = (-psi(i,j+1,0)+psi(i,j-1,0))/(2*dy)
LET yProbabilityCurrentDensity = psi(i,j,0)*pRe + psi(i,j,1)*pIm
END FUNCTION
! ---------- 3D graphics
EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
LET cosAx = COS(angleX)
LET sinAx = SIN(angleX)
LET cosAy = COS(angleY)
LET sinAy = SIN(angleY)
LET cx0 = xCenter
LET cy0 = yCenter
LET cz0 = zCenter
END SUB
EXTERNAL SUB plotLines3D(x,y,z) !shift*xRotateAx*yRotateAy*(shift^-1)
LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
PLOT LINES: x1+cx0, y1+cy0; !z=z1+cz0
END SUB
! ---------- draw wave
EXTERNAL SUB drawWave(drawMode)
DECLARE EXTERNAL FUNCTION kineticEnergy,potentialEnergy
DECLARE EXTERNAL SUB drawPsi2,drawCloud,draw3DPsi2,drawProbabilityCurrent,drawPsiPhase
SET DRAW MODE HIDDEN
CLEAR
IF drawMode=0 THEN CALL drawPsi2
IF drawMode=1 THEN CALL drawCloud
IF drawMode=2 THEN CALL draw3DPsi2(-PI/4,-PI/12) !(xRotateAngle,yRotateAngle)
IF drawMode=3 THEN CALL drawProbabilityCurrent
IF drawMode=4 THEN CALL drawPsiPhase
!--- caption
LET ke = kineticEnergy
LET pe = potentialEnergy
SET TEXT HEIGHT 10
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,-20 ,USING "time =#####.## (au) =####.#### (femto s)":sysTime,sysTime*auTime*1e15
PLOT TEXT, AT 0,-35 ,USING "norm of wave function =#.###############":norm
PLOT TEXT, AT 0,-50 ,USING "kinetic energy =#.###### (au) =##.### (eV)":ke,ke*27.2114
PLOT TEXT, AT 0,-65 ,USING "potential energy =#.###### (au) =##.### (eV)":pe,pe*27.2114
PLOT TEXT, AT 0,-80 ,USING "total energy =#.###### (au) =##.### (eV)":ke+pe,(ke+pe)*27.2114
PLOT TEXT, AT 0,-95 ,USING "Box =###.# x ###.# (au)":NNx*dx,NNx*dx
PLOT TEXT, AT 0,-110 :"wave packet - quantum electron dynamics 2D"
SET DRAW MODE EXPLICIT
END SUB
EXTERNAL SUB drawPsi2 !drawMode=0
MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
FOR j=0 TO NNy-1
LET col = (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*100
IF col>0.2 THEN
IF col>1 THEN LET col = 1
IF col<0 THEN LET col = 0
LET psicol(i,j) = 40+INT(col*100)
END IF
NEXT j
NEXT i
MAT PLOT CELLS,IN 0,0; 320,320 :psicol
END SUB
EXTERNAL SUB drawCloud !drawMode=1
MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
LET srnd(0) = RND
FOR i = 1 TO 1000
LET srnd(i) = srnd(i-1) + RND
NEXT i
FOR i = 0 TO 1000
LET srnd(i) = srnd(i)/srnd(1000)
NEXT i
LET m = 0
LET s = 0
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET cij = 0
LET s = s + (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*dx*dy
DO WHILE (s>srnd(m) AND m<1000)
LET m = m + 1
LET cij = cij + 1
LOOP
IF cij>0 THEN LET psicol(i,j) = 6 ! yellow
IF cij>3 THEN LET psicol(i,j) = 4 ! red
NEXT j
NEXT i
MAT PLOT CELLS,IN 0,0; 320,320 :psicol
END SUB
EXTERNAL SUB draw3DPsi2(xRotateAngle,yRotateAngle) !drawMode=2
DECLARE EXTERNAL SUB setRotateXYParameters,plotLines3D
LET sc = 320/NNx
LET zMag = 20
CALL setRotateXYParameters(xRotateAngle,yRotateAngle,NNx/2,NNy/2,0)
FOR j=0 TO NNy-1 STEP 2
FOR i=0 TO NNx-1
LET psi2 = (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*300
IF psi2>0.2 THEN
SET LINE COLOR 13 !psi^2:dark yellow
ELSEIF vv(i,j)>2 THEN
SET LINE COLOR 10 !potential:dark green
ELSE
SET LINE COLOR 3 !potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi2 + vv(i,j))) !(x,y,z)
NEXT i
PLOT LINES
NEXT j
FOR i=0 TO NNx-1 STEP 2
FOR j=0 TO NNy-1
LET psi2 = (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*300
IF psi2>0.2 THEN
SET LINE COLOR 13 !psi^2:dark yellow
ELSEIF vv(i,j)>2 THEN
SET LINE COLOR 10 !potential:dark green
ELSE
SET LINE COLOR 3 !potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi2 + vv(i,j))) !(x,y,z)
NEXT j
PLOT LINES
NEXT i
END SUB
EXTERNAL SUB drawProbabilityCurrent !drawMode=3
DECLARE EXTERNAL FUNCTION xProbabilityCurrentDensity,yProbabilityCurrentDensity
DECLARE EXTERNAL SUB drawPsi2
LET sc = 320/NNx
LET mag = 20000
CALL drawPsi2
SET LINE COLOR 4 !red
FOR i=1 TO NNx-2
FOR j=1 TO NNy-2
LET psi2 = (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*100
IF psi2>0.4 AND MOD(i,4)=0 AND MOD(j ,4)=0 THEN
LET xp = xProbabilityCurrentDensity(i,j)*dt*mag
LET yp = yProbabilityCurrentDensity(i,j)*dt*mag
PLOT LINES: i*sc, j*sc; i*sc+xp,j*sc+yp
END IF
NEXT j
NEXT i
PLOT LINES
END SUB
EXTERNAL SUB drawPsiPhase !drawMode=4
MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
FOR j=0 TO NNy-1
LET col = (psi(i,j,0)*psi(i,j,0)+psi(i,j,1)*psi(i,j,1))*100
IF col>0.2 AND col>RND+0.2 THEN
IF abs(psi(i,j,0))>abs(psi(i,j,1)) THEN
IF psi(i,j,0)>=0 THEN LET psicol(i,j) = 6 ELSE LET psicol(i,j) = 5
ELSE
IF psi(i,j,1)>=0 THEN LET psicol(i,j) = 13 ELSE LET psicol(i,j) = 8
END IF
END IF
NEXT j
NEXT i
MAT PLOT CELLS,IN 0,0; 320,320 :psicol
END SUB
END MODULE
|
|