c---------------------------------------------------------------------|
c FDTD3dp.f ( Maxwell Equation Simulator ) |
c |
c FDTD法によるスラブ座標系(x,y,z)でのプラズマ中における |
c マクスウェル方程式の２次元数値計算 |
c サブルーチン：PLASMA |
c 位置(x,y,z)での密度 F ，磁場 GX, GY, GZ，減衰率 dmp を与える |
c サブルーチン：YUDENTAI |
c 位置(x,y,z)を与えて，比誘電率(Sr, Si), 比透磁率 Qr を与える |
c サブルーチン：REISIN |
c x=Xminでの入射波形の時間分布を与える |
c サブルーチン：PROFILE |
c x=Xminでの入射波形のz方向分布を与える |
c x=Xmin, x=Xmax, z=Zmin, z=ZmaxでMurの２次吸収境界条件を用いる |
c プログラム作成者： 北條仁士， 筑波大学 2003年 1月 |
c---------------------------------------------------------------------|
implicit real*8(a-h,o-z)
dimension BX(1001,201,201), BY(1001,201,201), BZ(1001,201,201)
dimension EX(1001,201,201), EY(1001,201,201), EZ(1001,201,201)
dimension PX(1001,201,201), PY(1001,201,201), PZ(1001,201,201)
dimension EX1(1001,201,201),EY1(1001,201,201),EZ1(1001,201,201)
common /BLK1/w0,B0,OMEGA,R0
pi = 3.141592654d0
c---------------------------------------------------------------------|
c 数値計算パラメータの設定 |
c w0 = 3.0d11 のとき，長さの規格化 c/w0 = 1 mm , dx=dz= 0.1 mm |
c---------------------------------------------------------------------|
w0 = 3.0d10
B0 = 1.0d0
R0 = 4.0d0
dt = 0.05d0
dx = 0.1d0
dy = 0.1d0
dz = 0.1d0
wx = dt / dx
wy = dt / dy
wz = dt / dz
c---------------------------------------------------------------------|
c システムパラメータの設定 |
c 入力周波数：OMEGA は GHz で表示 |
c システム長：Xmin, Zmin, Xmax, Zmax はミリ単位で表示 |
c---------------------------------------------------------------------|
itmax = 3200
OMEGA = 2.45d0
Xmin = 0.0d0
Ymin = -5.0d0 + 1.0d-3
Zmin = -5.0d0
Xmax = 20.0d0
Ymax = 5.0d0
Zmax = 5.0d0
IMX = idint( (Xmax - Xmin) / dx + 0.5d0 )
JMX = idint( (Ymax - Ymin) / dy + 0.5d0 )
KMX = idint( (Zmax - Zmin) / dz + 0.5d0 )
IMDX = IMX / 200
JMDY = JMX / 100
KMDZ = KMX / 100
c---------------------------------------------------------------------|
c 電磁場の初期値の設定 |
c---------------------------------------------------------------------|
do 100 k=1,KMX+1
do 100 j=1,JMX+1
do 100 i=1,IMX+1
BX(i,j,k) = 0.0d0
BY(i,j,k) = 0.0d0
BZ(i,j,k) = 0.0d0
EX(i,j,k) = 0.0d0
EY(i,j,k) = 0.0d0
EZ(i,j,k) = 0.0d0
PX(i,j,k) = 0.0d0
PY(i,j,k) = 0.0d0
PZ(i,j,k) = 0.0d0
100 continue
do 110 k=1,KMX+1
do 110 j=1,JMX+1
do 110 i=1,IMX+1
EX1(i,j,k) = 0.0d0
EY1(i,j,k) = 0.0d0
EZ1(i,j,k) = 0.0d0
110 continue
c---------------------------------------------------------------------|
c 計算のスタート:時間 ： time = ( it ) * dt |
c---------------------------------------------------------------------|
it = 0
20 continue
it = it + 1
c---------------------------------------------------------------------|
c z=Zminで励起する入射波の条件設定 |
c---------------------------------------------------------------------|
time = dt * dfloat(it)
do k=2,KMX
do j=2,JMX
z = dfloat(k-1)*dz + Zmin
y = dfloat(j-1)*dy + Ymin
R = dsqrt( z*z + y*y )
IF( R<=R0 ) THEN
EX(1,j,k) = (1.0d0 - R/R0)**2*dsin( 2.0d9*pi*OMEGA/w0 * time )
ELSE
EX(1,j,k) = 0.0d0
END IF
end do
end do
c---------------------------------------------------------------------|
c 電磁場の磁場成分 BX, BY, BZ の計算 |
c---------------------------------------------------------------------|
do 200 k=1,KMX
do 200 j=1,JMX
do 200 i=1,IMX
x = dfloat(i-1)*dx + Xmin + dx/2.d0
y = dfloat(j-1)*dy + Ymin + dy/2.d0
z = dfloat(k-1)*dz + Zmin + dz/2.d0
CALL YUDENTAI( x, y, z, Sr, Si, Qr )
BX(i,j,k) = BX(i,j,k) + wz*( EY(i,j,k+1) - EY(i,j,k) ) / Qr
BX(i,j,k) = BX(i,j,k) - wy*( EZ(i,j+1,k) - EZ(i,j,k) ) / Qr
BY(i,j,k) = BY(i,j,k) + wx*( EZ(i+1,j,k) - EZ(i,j,k) ) / Qr
BY(i,j,k) = BY(i,j,k) - wz*( EX(i,j,k+1) - EX(i,j,k) ) / Qr
BZ(i,j,k) = BZ(i,j,k) + wy*( EX(i,j+1,k) - EX(i,j,k) ) / Qr
BZ(i,j,k) = BZ(i,j,k) - wx*( EY(i+1,j,k) - EY(i,j,k) ) / Qr
200 continue
c---------------------------------------------------------------------|
c 誘起電流成分 PX, PY, PZ の計算 |
c---------------------------------------------------------------------|
do 250 k=2,KMX
do 250 j=2,JMX
do 250 i=2,IMX
x = dfloat(i-1)*dx + Xmin
y = dfloat(j-1)*dy + Ymin
z = dfloat(k-1)*dz + Zmin
CALL PLASMA( x, y, z, dmp, F, A, GX, GY, GZ )
C1 = ( 1.d0 - dmp*dt/2.d0 ) / ( 1.d0 + dmp*dt/2.d0)
D1 = 1.d0 / ( 1.d0 + dmp*dt/2.d0)
PX(i,j,k) = C1*PX(i,j,k) + D1*dt*F*EX(i,j,k)
PX(i,j,k) = PX(i,j,k) - D1*dt*A*( PY(i,j,k)*GZ - PZ(i,j,k)*GY )
PY(i,j,k) = C1*PY(i,j,k) + D1*dt*F*EY(i,j,k)
PY(i,j,k) = PY(i,j,k) - D1*dt*A*( PZ(i,j,k)*GX - PX(i,j,k)*GZ )
PZ(i,j,k) = C1*PZ(i,j,k) + D1*dt*F*EZ(i,j,k)
PZ(i,j,k) = PZ(i,j,k) - D1*dt*A*( PX(i,j,k)*GY - PY(i,j,k)*GX )
250 continue
c---------------------------------------------------------------------|
c 電磁場の電場成分 EX, EY, EZ の計算 |
c---------------------------------------------------------------------|
do 300 k=2,KMX
do 300 j=2,JMX
do 300 i=2,IMX
x = dfloat(i-1)*dx + Xmin
y = dfloat(j-1)*dy + Ymin
z = dfloat(k-1)*dz + Zmin
CALL YUDENTAI( x, y, z, Sr, Si, Qr )
V = ( OMEGA*2.0d9*pi/w0 ) * Si / Sr
C = ( 1.d0 - V*dt/2.d0 ) / ( 1.d0 + V*dt/2.d0)
D = 1.d0 / Sr / ( 1.d0 + V*dt/2.d0)
EX(i,j,k) = C*EX(i,j,k) - D*dt*PX(i,j,k)
EX(i,j,k) = EX(i,j,k) + D*wy*( BZ(i,j,k) - BZ(i,j-1,k) )
EX(i,j,k) = EX(i,j,k) - D*wz*( BY(i,j,k) - BY(i,j,k-1) )
EY(i,j,k) = C*EY(i,j,k) - D*dt*PY(i,j,k)
EY(i,j,k) = EY(i,j,k) + D*wz*( BX(i,j,k) - BX(i,j,k-1) )
EY(i,j,k) = EY(i,j,k) - D*wx*( BZ(i,j,k) - BZ(i-1,j,k) )
EZ(i,j,k) = C*EZ(i,j,k) - D*dt*PZ(i,j,k)
EZ(i,j,k) = EZ(i,j,k) + D*wx*( BY(i,j,k) - BY(i-1,j,k) )
EZ(i,j,k) = EZ(i,j,k) - D*wy*( BX(i,j,k) - BX(i,j-1,k) )
300 continue
c---------------------------------------------------------------------|
c YZ面（x=Xmin）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
ux = ( wx - 1.d0 )/( wx + 1.d0 )
uy = ( wy - 1.d0 )/( wy + 1.d0 )
uz = ( wz - 1.d0 )/( wz + 1.d0 )
d1 = 0.005d0
ds = 1.d0 - d1
do 400 k=1,KMX+1
do 400 j=1,JMX+1
EY(1,j,k)= ds*EY1(2,j,k) + ux*(EY(2,j,k)-EY1(1,j,k))
EZ(1,j,k)= ds*EZ1(2,j,k) + ux*(EZ(2,j,k)-EZ1(1,j,k))
400 continue
c---------------------------------------------------------------------|
c YZ面（x=Xmax）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
do 410 k=1,KMX+1
do 410 j=1,JMX+1
EY(IMX+1,j,k)= ds*EY1(IMX,j,k) + ux*(EY(IMX,j,k)-EY1(IMX+1,j,k))
EZ(IMX+1,j,k)= ds*EZ1(IMX,j,k) + ux*(EZ(IMX,j,k)-EZ1(IMX+1,j,k))
410 continue
c---------------------------------------------------------------------|
c XZ面（y=Ymax）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
do 420 k=1,KMX+1
do 420 i=1,IMX+1
EX(i,JMX+1,k)= ds*EX1(i,JMX,k) + uy*(EX(i,JMX,k)-EX1(i,JMX+1,k))
EZ(i,JMX+1,k)= ds*EZ1(i,JMX,k) + uy*(EZ(i,JMX,k)-EZ1(i,JMX+1,k))
420 continue
c---------------------------------------------------------------------|
c XZ面（y=Ymin）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
do 430 k=1,KMX+1
do 430 i=1,IMX+1
EX(i,1,k)= ds*EX1(i,2,k) + uy*(EX(i,2,k)-EX1(i,1,k))
EZ(i,1,k)= ds*EZ1(i,2,k) + uy*(EZ(i,2,k)-EZ1(i,1,k))
430 continue
c---------------------------------------------------------------------|
c XY面（z=Zmax）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
do 440 j=1,JMX+1
do 440 i=1,IMX+1
EX(i,j,KMX+1)= ds*EX1(i,j,KMX) + uz*(EX(i,j,KMX)-EX1(i,j,KMX+1))
EY(i,j,KMX+1)= ds*EY1(i,j,KMX) + uz*(EY(i,j,KMX)-EY1(i,j,KMX+1))
440 continue
c---------------------------------------------------------------------|
c XY面（z=Zmin）でのHigdonの1次の吸収境界条件 |
c---------------------------------------------------------------------|
do 450 j=1,JMX+1
do 450 i=1,IMX+1
EX(i,j,1)= ds*EX1(i,j,2) + uz*(EX(i,j,2)-EX1(i,j,1))
EY(i,j,1)= ds*EY1(i,j,2) + uz*(EY(i,j,2)-EY1(i,j,1))
450 continue
c----------------------------------------------------------------------
do 490 k=1,KMX+1
do 490 j=1,JMX+1
do 490 i=1,IMX+1
EX1(i,j,k) = EX(i,j,k)
EY1(i,j,k) = EY(i,j,k)
EZ1(i,j,k) = EZ(i,j,k)
490 continue
c---------------------------------------------------------------------|
c 電磁場の十分に小さい値はゼロで置き換える |
c it = itmax になれば，計算を終了する |
c---------------------------------------------------------------------|
do 500 k=1,KMX+1
do 500 j=1,JMX+1
do 500 i=1,IMX+1
if( dabs(EX(i,j,k)) .lt. 1.0d-80 ) EX(i,j,k) = 0.d0
if( dabs(EY(i,j,k)) .lt. 1.0d-80 ) EY(i,j,k) = 0.d0
if( dabs(EZ(i,j,k)) .lt. 1.0d-80 ) EZ(i,j,k) = 0.d0
if( dabs(BX(i,j,k)) .lt. 1.0d-80 ) BX(i,j,k) = 0.d0
if( dabs(BY(i,j,k)) .lt. 1.0d-80 ) BY(i,j,k) = 0.d0
if( dabs(BZ(i,j,k)) .lt. 1.0d-80 ) BZ(i,j,k) = 0.d0
500 continue
if( it.eq.itmax ) go to 600
go to 20
c---------------------------------------------------------------------|
c 電磁場の値のデータファイルを作る |
c---------------------------------------------------------------------|
600 continue
open( 11, FILE='exdata1', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(11,911) (EX(61,j,k),k=1,KMX+1,KMDZ)
911 format(1h , e15.5, 100(2x, e15.5) )
end do
close(11)
open( 12, FILE='exdata2', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(12,912) (EX(81,j,k),k=1,KMX+1,KMDZ)
912 format(1h , e15.5, 100(2x, e15.5) )
end do
close(12)
open( 13, FILE='exdata3', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(13,913) (EX(101,j,k),k=1,KMX+1,KMDZ)
913 format(1h , e15.5, 100(2x, e15.5) )
end do
close(13)
open( 14, FILE='exdata4', STATUS='unknown' )
do i=1,IMX+1,IMDX
write(14,914) (EX(i,51,k),k=1,KMX+1,KMDZ)
914 format(1h , e15.5, 100(2x, e15.5) )
end do
close(14)
c open( 15, FILE='xbunpu', STATUS='unknown' )
c do i=1,IMX+1
c x = dfloat(i-1)*dx + Xmin
c write(15,915) x, EY(i,51,51)
c 915 format(1h , e15.5, 2x, e15.5 )
c end do
c close(15)
c---------------------------------------------------------------------|
c メインプログラムの終了 |
c---------------------------------------------------------------------|
9999 stop
end
c---------------------------------------------------------------------|
c サブルーチン：PLASMA |
c 位置(x,z)での密度 F ，磁場 GX, GY, GZ，減衰率 dmp を与える |
c 最大密度 ：DENTY は cm-3 で表示 |
c 磁場 ：B0, BEX0 は Tesla で表示 |
c---------------------------------------------------------------------|
SUBROUTINE PLASMA( x, y, z, dmp, F, A, GX, GY, GZ )
implicit real*8(a-h, o-z)
common /BLK1/w0,B0,OMEGA
pai = 3.141592654d0
DENTY = 2.0d11
BEX0 = 0.0d0
dmp = 0.0d0
F0 = 3.181d9 * DENTY / ( w0*w0 )
A = 1.76d11 * B0 / w0
FNC = ( 1.d0 + dtanh((x-10.d0)/3.d0) )/2.d0
F = F0*FNC
c IF( x .GE. 10.d0 ) F = F0*( x - 10.d0 )/10.d0
c IF( x .LT. 10.d0 ) F = 0.d0
GX = 0.0d0
GZ = 0.0d0
GY = BEX0 / B0
return
end
c---------------------------------------------------------------------|
c サブルーチン：YUDENTAI |
c 位置(x,z)での比誘電率(Sr, Si), 比透磁率 Qrを与える |
c---------------------------------------------------------------------|
SUBROUTINE YUDENTAI( x, y, z, Sr, Si, Qr )
implicit real*8(a-h,o-z)
common /BLK1/w0,B0,OMEGA,R0
pi = 3.141592654d0
Sr = 1.d0
Si = 0.d0
Qr = 1.d0
R = dsqrt( z*z + y*y )
IF( R>=R0 ) THEN
Sr = 1.d0
Si = 0.d0
END IF
999 return
end
c---------------------------------------------------------------------|
c サブルーチン：REISIN |
c x=Xminでの入射波の波形の時間分布を与える |
c---------------------------------------------------------------------|
c SUBROUTINE REISIN( time, Bunpt )
c implicit real*8(a-h,o-z)
c common /BLK1/w0,B0,OMEGA
c pai = 3.141592654d0
c Freq = 2.0d9*pai*OMEGA / w0
c Bunpt = dsin( Freq * time )
c return
c end
c---------------------------------------------------------------------|
c サブルーチン：PROFILE |
c x=Xexcでの入射波の波形のz方向分布を与える |
c---------------------------------------------------------------------|
c SUBROUTINE PROFILE( y, z, Bunp )
c implicit real*8(a-h,o-z)
c pai = 3.141592654d0
c yy = y - 5.d0
c zz = z - 5.d0
c R = dsqrt( yy**2 + zz**2 )
c Bunp = dexp( - ( R / 2.0d0 )**2 )
c IF( R.LE.2.d0 ) Bunp = ( 1.d0 - (R/2.d0)**2 )**2
c IF( R.GT.2.d0 ) Bunp = 0.d0
c return
c end