|
|
> No.4251[元記事へ]
mikeさんへのお返事です。
２進モードで、今度は見やすいです。
走り書きに近く、どこか壊れているかも知れませんので、点検して下さい。
とりあえず、波動関数の、らせん形が、見えます。
! ========= Quantum Electron Dynamics 1D ==========
!
! wavePacketQED1D.bas
!
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.01.31 created
!
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DECLARE FUNCTION systemTime,norm
DECLARE FUNCTION kineticEnergy,potentialEnergy
!
! ---------- QED1D 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).)
!
!--------------
DIM pV(1 TO 4),Axys(1 TO 4,1 TO 4),P3D(1 TO 4,1 TO 4)
DIM rotx(1 TO 4,1 TO 4),shxyz(1 TO 4,1 TO 4),shres(1 TO 4,1 TO 4)
MAT Axys=IDN
MAT rotx=IDN
MAT shxyz=IDN
MAT shres=IDN
LET pV(4)=1
!
LET Sx=1
LET Sy=1
LET Sz=1
LET Ax=0 !開始のｚ軸方向( 画面垂直０度からｘ軸回転成分)
LET Ay=0 ! 〃 〃 ( 〃 〃 ｙ軸回転成分)
LET Az=0 !ｚ軸回転成分
!
DIM phRe(500),phIm(500) ! wave function
DIM vv(500) ! external potential
DIM bRe(500),bIm(500) ! b vector in kxStep
DIM uRe(500),uIm(500) ! u vector in kxStep
LET NNx = 400 ! max number of phRe(),phIm(),..
LET sysTime = 0.0 ! (au) au : atomic unit (hbar=1, me=1, e=1)
LET dx = 0.5 ! (au) x 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)
!
CALL setInitialCondition(1) !(menu)
!
FOR it=1 TO 1000
FOR i9=1 TO 10
IF mlb=0 THEN CALL evolveWave
NEXT i9
CALL drawWave
!IF (MOD(it,10)=0) THEN
! PRINT USING "time =#####.## (au) norm =#.###############":systemTime,norm
!END IF
mouse poll mx,my,mlb,mrb
IF mrb=1 THEN EXIT FOR
! LET Ay=Ay+1e-2 !trim. ｙ軸で、自動で回す場合
NEXT it
! ----- set initial condition
!
SUB setInitialCondition(menu) !public
LET sysTime = 0.0
IF menu=1 THEN
CALL setGaussianWave(50.0,10.0,0.0) !(wavePos,waveWidth,momentum)
CALL setparabolicPotential(0.5*(NNx-1)*dx,0.0001) ! v(x)=k*(x-x0)^2
ELSEIF menu=2 THEN
CALL setGaussianWave(50.0,10.0,1.0) !(wavePos,waveWidth,momentum)
CALL setparabolicPotential(0.5*(NNx-1)*dx,0.0) ! v(x)=k*(x-x0)^2
END IF
! set window
LET w = NNx*dx
LET b = 20
SET WINDOW -b,w+b,-b,w+b
!--
CALL mat_shxyz(shxyz,-w/2,-w/2, 0) !画面中心を原点へ移動する行列 shxyz 作成
CALL mat_shxyz(shres, w/2, w/2, 0) !画面中心を原点から元へ戻す行列 shres 作成
END SUB
SUB setGaussianWave(wavePos,waveWidth,momentum)
FOR i=1 TO NNx-2
LET x = i*dx
LET phAb = EXP(-((x-wavePos)/(2.0*waveWidth))^2)
LET phPh = momentum*x
LET phRe(i) = phAb*cos(phPh)
LET phIm(i) = phAb*sin(phPh)
NEXT i
LET phRe(0) = 0.0
LET phIm(0) = 0.0
LET phRe(NNx-1) = 0.0
LET phIm(NNx-1) = 0.0
CALL normalize
END SUB
SUB setparabolicPotential(x0,k) ! v(x)=k*(x-x0)^2
FOR i=0 TO NNx-1
LET x = i*dx
LET vv(i) = k*(x-x0)^2
NEXT i
END SUB
! ----- evolve Wave
!
SUB evolveWave !public
DECLARE EXTERNAL SUB phaseStep,kxStep
CALL phaseStep
CALL kxStep
CALL phaseStep
LET sysTime = sysTime + dt
END SUB
SUB kxStep !{exp(i dt(delta/2)} semi-implicit method
LET ai = 4*dx*dx/dt
LET aaAb = 4+ai*ai
FOR i=1 TO NNx-2 !set b-vector
LET bRe(i) = 2*phRe(i)-ai*phIm(i) - phRe(i+1) - phRe(i-1)
LET bIm(i) = 2*phIm(i)+ai*phRe(i) - phIm(i+1) - phIm(i-1)
NEXT i
! forward elimination
LET uRe(1) = -2/aaAb
LET uIm(1) = -ai/aaAb
LET phRe(1) = bRe(1)*uRe(1) - bIm(1)*uIm(1)
LET phIm(1) = bIm(1)*uRe(1) + bRe(1)*uIm(1)
FOR i=2 TO NNx-2
LET auAb = (-2-uRe(i-1))*(-2-uRe(i-1))+(ai-uIm(i-1))*(ai-uIm(i-1))
LET uRe(i) = (-2-uRe(i-1))/auAb
LET uIm(i) = -(ai-uIm(i-1))/auAb
LET phRe(i) = (bRe(i)-phRe(i-1))*uRe(i) - (bIm(i)-phIm(i-1))*uIm(i)
LET phIm(i) = (bRe(i)-phRe(i-1))*uIm(i) + (bIm(i)-phIm(i-1))*uRe(i)
NEXT i
! backward substitution
FOR i=NNx-3 TO 1 STEP -1
LET phRe(i) = phRe(i) - (phRe(i+1)*uRe(i) - phIm(i+1)*uIm(i))
LET phIm(i) = phIm(i) - (phRe(i+1)*uIm(i) + phIm(i+1)*uRe(i))
NEXT i
END SUB
SUB phaseStep !{exp(-i(dt/2)V} evolve 0.5dt
FOR i=1 TO NNx-2
LET th = 0.5*dt*vv(i)
LET cs = COS(th)
LET sn = SIN(th)
LET phr = phRe(i)
LET phi = phIm(i)
LET phRe(i) = cs*phr+sn*phi
LET phIm(i) = -sn*phr+cs*phi
NEXT i
END SUB
! ----- utility
!
FUNCTION systemTime
LET systemTime = sysTime
END FUNCTION
SUB normalize
LET a = norm
FOR i=1 TO NNx-2
LET phRe(i) = phRe(i)/a
LET phIm(i) = phIm(i)/a
NEXT i
END SUB
FUNCTION norm !public SQR(<psi|psi>)
LOCAL i,a2
LET a2 = 0.0
FOR i=1 TO NNx-2
LET a2 = a2 + (phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
NEXT i
LET norm = SQR(a2)
END FUNCTION
FUNCTION kineticEnergy !<psi| d^2/dx^2 |psi>
LET s = 0.0
FOR i=1 TO NNx-2
LET hphRe = (2.0*phRe(i)-phRe(i+1)-phRe(i-1))/(2.0*dx*dx)
LET hphIm = (2.0*phIm(i)-phIm(i+1)-phIm(i-1))/(2.0*dx*dx)
LET s = s + (phRe(i)*hphRe + phIm(i)*hphIm)*dx
NEXT i
LET kineticEnergy = s
END FUNCTION
FUNCTION potentialEnergy !<psi| V(x) |psi>
LET s = 0.0
FOR i=1 TO NNx-2
LET s = s + vv(i)*(phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
NEXT i
LET potentialEnergy = s
END FUNCTION
! ----- draw wave
!
SUB drawWave
LET yh = (NNx*dx+40)*0.4
LET sth = SIN(PI/12)
!
SET DRAW MODE HIDDEN
CLEAR
CALL control_
!---plot V(x),<psi|psi>,|psi>
SET LINE COLOR 1 ! black : PLOT potential V(x);
FOR i=0 TO NNx-1
CALL line3D(dx*i, vv(i)*100+yh, 0) !PLOT LINES: dx*i, vv(i)*100+yh;
NEXT i
PLOT LINES
!--
SET LINE COLOR 4 ! red : plot probability <psi|psi>
FOR i=0 TO NNx-1
LET ph2 = (phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
CALL line3D( dx*i, ph2*1000+yh, 0) !PLOT LINES: dx*i, ph2*1000+yh;
NEXT i
PLOT LINES
!--
FOR i=0 TO NNx-1 !plot wave function |psi(x,t)>
LET x = i*dx
LET y = phRe(i)*50
LET z = phIm(i)*50
IF phIm(i)>=0 THEN SET LINE COLOR 2 ELSE SET LINE COLOR 10
CALL line3D( x, y+yh, z) !PLOT LINES: x-sth*z, y+yh;
!PLOT LINES: dx*i-sth*phIm(i)*50, phRe(i)*50+yh;
NEXT i
PLOT LINES
!--
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0, yh+22 :"potential V(x)"
SET TEXT COLOR 4 ! red
PLOT TEXT, AT 0, yh+6 :"<psi|psi>"
SET TEXT COLOR 2 ! blue
PLOT TEXT, AT 0, yh-10 :"|psi(x,t)>"
!---caption
LET ke = kineticEnergy
LET pe = potentialEnergy
SET TEXT HEIGHT 5
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,48 ,USING "time =#####.## (au) =####.#### (femto s)":sysTime,sysTime*auTime*1e15
PLOT TEXT, AT 0,40 ,USING "norm of wave function =#.###############":norm
PLOT TEXT, AT 0,32 ,USING "kinetic energy =#.###### (au) =##.### (eV)":ke,ke*27.2114
PLOT TEXT, AT 0,24 ,USING "potential energy =#.###### (au) =##.### (eV)":pe,pe*27.2114
PLOT TEXT, AT 0,16 ,USING "total energy =#.###### (au) =##.### (eV)":ke+pe,(ke+pe)*27.2114
PLOT TEXT, AT 0, 8 ,USING "xLength =####.# (au) =###.### (nm)":NNx*dx,NNx*dx*auLength*1e9
PLOT TEXT, AT 0, 0 :"wave packet - quantum electron dynamics 1D"
SET DRAW MODE EXPLICIT
END SUB
SUB control_
PLOT label,AT 20, 210:"左 click 一時停止、drag 手動回転。右 click 終了。"
!-----click_drag-----
IF mlb=1 THEN
!!LET Ax= -(my-mybak)*PI/100 !※Ｙ軸回転のみにする。 !ドラッグ方向から、軸方向と回転量
LET Ay= +(mx-mxbak)*PI/100
END IF
LET mxbak=mx
LET mybak=my
!-----
LET ar0=SQR(Ax^2+Ay^2) !回転の角度(∝マウス・ドラッグの長さ)
IF ar0<>0 THEN
LET DIRar0=ANGLE(Ax,Ay) !軸の角度
CALL mat_rotx(rotx, ar0)
MAT Axys=Axys*ROTATE(-DIRar0)*rotx*ROTATE(DIRar0) !ドラッグ累積 (方向,回転)
LET Ax=0
LET Ay=0
END IF
MAT P3D=shxyz*ROTATE(Az)*Axys*shres
END SUB
SUB axes3D(x1,y1,z1, x2,y2,z2, a$ )
CALL line3D(x1,y1,z1)
CALL line3D(x2,y2,z2)
PLOT TEXT,AT pV(1),pV(2) :a$ ! PEN-off
END SUB
SUB line3D(x,y,z)
LET pV(1)=x*Sx !描画目盛は、全方向等しくないと、回転で、形が保てない。
LET pV(2)=y*Sy !スケール Sx,Sy,Sz の違いは、入力の倍率として、行なう。
LET pV(3)=z*Sz !入力 z 座標は 出力 x,y に反映。出力ｚは 描画不可。
LET pV(4)=1 ! shxyzM …shxyzP で必要。
MAT pV=pV*P3D
PLOT LINES: pV(1),pV(2); ! PEN-on
END SUB
!---------------------------------
! ｘ軸で 回転する行列 → 配列引数
!(x,y,z,1)| 1, 0, 0, 0 |
! | 0, cos(a), sin(a), 0 |
! | 0,-sin(a), cos(a), 0 |
! | 0, 0, 0, 1 |
!---------------------------------
SUB mat_rotx(m(,), a)
LET m(2,2)=COS(a)
LET m(3,2)=-SIN(a)
LET m(2,3)=SIN(a)
LET m(3,3)=COS(a) !他の要素は、呼出し側で管理
END SUB
!-----------------------------
! 平行移動。(sx,sy,sz)
!(x,y,z,1)| 1, 0, 0, 0 |
! | 0, 1, 0, 0 |
! | 0, 0, 1, 0 |
! | sx, sy, sz, 1 |
!-----------------------------
SUB mat_shxyz(m(,), sx,sy,sz)
LET m(4,1)=sx
LET m(4,2)=sy
LET m(4,3)=sz !他の要素は、呼出し側で管理
END SUB
END
|
|