|
|
電磁波（光）の伝搬（FDTD法）の雛形プログラムを公開します。
雛形なので、なるべくFDTD法の本質的な部分だけをわかりやすく書いたつもりです。
（高速化や汎用性をねらうと、mikeの力量ではプログラムが読み難くなりますので）
光の伝搬はMaxwellの電磁場の方程式を用いて現在の電場Ezと磁場Hx,Hyからdt時間後の電場Ezと磁場Hx,Hyを求めます。
これを繰り返して時間発展させていきます。
３次元（Ex,Ey,Ez,Hx,Hy,Hz)だと計算時間が長く(Ez(200,200,200) etc.)なり、
波が動いてのを見ることが難しいので、2次元にします(Ez(200,200) etc.)。
Ezと相互に結びつきの強いHx,Hyの３つのみで、矛盾なく２次元にすることができます。
この雛形プログラムでは、パルス状の平面波を誘電体(誘電率εr=4 or 屈折率n=2)にあて、そのふるまいを観察します。
誘電体での波の反射，透過、散乱が見られます。透過のとき、波の速度が1/2のなることがわかります。
表示の説明：
青のグリッドは電場Ezを鳥瞰したものです。青色は真空中（あるいは空気中）、赤い部分はガラスのような屈折率2の誘電体を表します。
長さや時間の間には、時間刻みdt、空間刻みdxの間にdx=c*dt（c=1）の関係があります。
（たとえばdt=1psとするとc=3e8 m/sとしてdx=0.3mmになります）
試験環境：
本プログラムは十進BASIC 0.6.6.0 / macOS 10.7.5と
十進BASIC Ver 7.7.8 / windows 10でテストしました。
!
! ========= FDTD 2D ==========
!
! waveFDTD2D.bas
! Mitsuru Ikeuchi (c) copyleft
!
! ver 0.0.1 2017.02.06
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB fdtd2d.setInitialCondition, fdtd2d.evolveField, fdtd2d.drawField
DECLARE EXTERNAL FUNCTION fdtd2d.systemTime
!
!CALL setInitialCondition(epsilon,mue,sigma)
CALL setInitialCondition(4.0,1.0,0.0) !dielectric material
!CALL setInitialCondition(1.0,4.0,0.0) !Magnetic material
!CALL setInitialCondition(1.0,1.0,0.01) !resistor
!CALL setInitialCondition(1.0,1.0,0.0) !vacuum or air
!CALL setInitialCondition(1000.0,1.0,0.1) ! metal
!
FOR it=1 TO 1000
CALL evolveField
CALL evolveField
CALL drawField
NEXT it
!
END
! ---------- FDTD2D module ----------
!
! time evolution : from Maxwell equation,
! d(mu H)/dt = - rot E (B = mu H)
! d(ep E)/dt = -j + rot H (D=ep E, j=sg E)
! FDTD:Finite-difference time-domain method
! 2D subset: Ez,Hx,Hy - finite difference approximation
!
MODULE fdtd2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, evolveField, drawField
PUBLIC FUNCTION systemTime
SHARE NUMERIC sysTime, dt,epsilon,reg,sgm,Jz,omega,theta,NNx,NNy
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC Ez(320,320) ! electric field
SHARE NUMERIC Hx(320,320),Hy(320,320) ! magnetic field
SHARE NUMERIC ep(320,320) ! dielectric constant
SHARE NUMERIC mu(320,320) ! permeability
SHARE NUMERIC sg(320,320) ! conductivity
LET sysTime = 0 ! system time
LET dt = 1 ! time step
LET dx = 1 ! x-division dx=dy=dz=1 no use
LET omega = PI/32 ! generator wave phase angle change
LET theta = 0 ! generator phase angle
LET NNx = 241 ! max number of Ez(i,) etc.
LET NNy = 241 ! max number of Ez(,j) etc.
!
! set initial condition
!
EXTERNAL SUB setInitialCondition(epsilon,mue,sigma)
LET theta = 0.0
!clear field
FOR i=0 TO NNx
FOR j=0 TO NNy
LET Ez(i,j) = 0
LET Hx(i,j) = 0
LET Hy(i,j) = 0
LET ep(i,j) = 1
LET mu(i,j) = 1
LET sg(i,j) = 0
NEXT j
NEXT i
!set material
FOR i=0 TO NNx
FOR j=0 TO NNy
IF (i>=60 AND i<140 AND j>=50 AND j<150) THEN
LET ep(i,j) = epsilon
LET mu(i,j) = mue
LET sg(i,j) = sigma
END IF
NEXT j
NEXT i
! set window
LET xmargin = 60
LET ymargin = 100
SET WINDOW -xmargin,400-xmargin,-ymargin,400-ymargin
END SUB
!
! evolve Field Ez,Hx,Hy
!
EXTERNAL SUB evolveField
DECLARE EXTERNAL SUB generateEz,evolveEz,evolveHxHy
LET sysTime = sysTime + dt
CALL evolveEz
CALL generateEz(2)
CALL evolveHxHy
END SUB
!
EXTERNAL SUB generateEz(nwave)
LET theta = theta + omega*dt
LET Ezt = SIN(theta)
LET a = 0
IF theta<2*PI*nwave THEN
LET a=1
ELSEIF theta<2*PI*nwave+0.13*PI THEN
LET a=COS(theta)
END IF
IF theta<2*PI*nwave+0.5*PI THEN
FOR j=0 TO NNy
LET Ez(0,j) = a*a*Ezt
NEXT j
END IF
END SUB
!
EXTERNAL SUB evolveEz ! dD/dt=rotH + J ,D=epsi*E
LET a = 1/(2*SQR(2))
FOR i=0 TO NNx-1
FOR j=0 TO NNy-1
LET Ez(i,j)=Ez(i,j)-sg(i,j)*0.5*dt*Ez(i,j)+(a/ep(i,j))*((Hy(i+1,j)-Hy(i,j))-(Hx(i,j+1)-Hx(i,j)))
NEXT j
NEXT i
END SUB
!
EXTERNAL SUB evolveHxHy ! dB/dt=-rotE ,B=mue*H
!no reflection boundary
LET a = 1/(2*SQR(2))
LET c = (1-a)/(1+a)
FOR i=0 TO NNx-1
LET Hx(i,0) = c*Hx(i,0) + Hx(i,1)
LET Hx(i,NNy) = c*Hx(i,NNy) + Hx(i,NNy-1)
NEXT i
FOR j=0 TO NNy-1
LET Hy(0,j) = c*Hy(0,j) + Hy(1,j)
LET Hy(NNx,j) = c*Hy(NNx,j) + Hy(NNx-1,j)
NEXT j
!
FOR i=0 TO NNx-1
FOR j=1 TO NNy-1
LET Hx(i,j) = Hx(i,j) - (a/mu(i,j))*(Ez(i,j)-Ez(i,j-1))
NEXT j
NEXT i
FOR i=1 TO NNx-1
FOR j=0 TO NNy-1
LET Hy(i,j) = Hy(i,j) + (a/mu(i,j))*(Ez(i,j)-Ez(i-1,j))
NEXT j
NEXT i
!no reflection boundary
FOR i=0 TO NNx-1
LET Hx(i,0) = Hx(i,0) - c*Hx(i,1)
LET Hx(i,NNy) = Hx(i,NNy) - c*Hx(i,NNy-1)
NEXT i
FOR j=0 TO NNy-1
LET Hy(0,j) = Hy(0,j) - c*Hy(1,j)
LET Hy(NNx,j) = Hy(NNx,j) - c*Hy(NNx-1,j)
NEXT j
END SUB
!
! utility
!
EXTERNAL FUNCTION systemTime
LET systemTime = sysTime
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 field Ez,Hx,Hy
!
EXTERNAL SUB drawField
DECLARE EXTERNAL SUB drawEz
!
SET DRAW MODE HIDDEN
CLEAR
!CALL setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
CALL drawEz(-PI/6,-PI/12) !(xRotateAngle,yRotateAngle)
!draw caption
SET TEXT HEIGHT 8
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,-20 ,USING "time =#####.# ":sysTime
PLOT TEXT, AT 0,-35 ,USING "field size = #### x ####":NNx,NNy
PLOT TEXT, AT 0,-50 :"2-dimensional FDTD"
!
SET DRAW MODE EXPLICIT
END SUB
!
EXTERNAL SUB drawEz(xRotateAngle,yRotateAngle)
DECLARE EXTERNAL SUB setRotate,plotLines3D
LET sc = 1.0
LET zMag = 20
!CALL setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
CALL setRotateXYParameters(xRotateAngle,yRotateAngle,NNx/2,NNy/2,0)
SET COLOR 2 ! blue : Ez
FOR j=0 TO NNy-1 STEP 5
FOR i=0 TO NNx-1
IF Ez(i,j)>=0 THEN SET COLOR 2 ELSE SET COLOR 9 !>=0:bule <0:dark blue
IF ep(i,j)>1 OR mu(i,j)>1 OR sg(i,j)>0 THEN SET COLOR 4 !something material
CALL plotLines3D(i*sc,j*sc,zMag*Ez(i,j)) !(x,y,z)
NEXT i
PLOT LINES
NEXT j
FOR i=0 TO NNx-1 STEP 5
FOR j=0 TO NNy-1
IF Ez(i,j)>=0 THEN SET COLOR 2 ELSE SET COLOR 9
IF ep(i,j)>1 OR mu(i,j)>1 OR sg(i,j)>0 THEN SET COLOR 4
CALL plotLines3D(i*sc,j*sc,zMag*Ez(i,j))
NEXT j
PLOT LINES
NEXT i
END SUB
!
END MODULE
|
|