[005]2次元-電子動力学法プログラムの公開

 投稿者:mike  投稿日:2017年 2月19日(日)07時08分2秒
   2次元-電子動力学法の雛形プログラム(005wavePacketQED2D.bas)を公開します。
 ポテンシャル中の電子の波動関数psi(x,y,t)とポテンシャルV(x,y)からdt時間後のpsi(x,y,t+dt)を求めます。
これを繰り返して時間発展させていきます。この波動関数から、電子の存在確率やエネルギーなどが解ります。
 2次元になると、複素関数の表現方法が難しく、ここでは存在確率<psi|psi>しか表現できていません。
何か良い方法があれば、教えてください。主プログラム中のCALL drawWave(2)をCALL drawWave(4)に変更すると、
位相が表現されていますが、成功しているとは思えません。
 この雛形プログラムでは、放物型ポテンシャル中の電子を追跡します。古典的な放物型ポテンシャル中のボールの動きと
類似しているように見えます。全エネルギーは一定になるはずですが、わずかに変動が見られます。
1次元と比較すると、およそ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



 

Re: [005]2次元-電子動力学法プログラムの公開

 投稿者:mike  投稿日:2017年 2月20日(月)10時45分24秒
  > No.4274[元記事へ]

2次元-電子動力学法の雛形プログラム(005wavePacketQED2D.bas)の動きが遅いので、BASIC Accで加速しようと
コンパイルしようとしたところ、コンパイルエラーを出してしまいました。改訂版を再掲載します。
全計算にかかった時間を計測すると、下記のようになりました。

450.21s  BASIC 6.6.0 / MAC OS 10.7.5 / mac mini Intel Core-i7 (2.7GHz)  4GB
332.56s  BASIC 7.7.8 / windows 10 64bit / FMV-A561D CPU intel core-i5 2400M (2.5GHz) 4GB
123.53s  BASIC Acc 0.9.8.1 / windows 10 64bit / FMV-A561D CPU intel core-i5 2400M (2.5GHz) 4GB
BASIC Accで2~3倍加速されるようです。これで動いて見えます。

----------------

!
! ========= Quantum Electron Dynamics 2D ==========
!
! 005wavePacketQED2D.bas
!    Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1   2017.02.19  created
! ver 0.0.2   2017.02.20  bug fixed (BASIC Acc 0.9.8.1 compiling error fixed)
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB qed2d.setInitialCondition, qed2d.evolveWave, qed2d.drawWave

LET t0 = TIME
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
!PRINT TIME - t0
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 NNy-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*dx*dx/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 FUNCTION kineticEnergy,potentialEnergy
   DECLARE EXTERNAL SUB normalize
   LET h2 = 2.0*dx*dx
   LET ee = kineticEnergy + potentialEnergy
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-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 NNy-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
 

戻る