c---------------------------------------------------------------------|
c FDTD3dm1.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 Conducting Wall Bounary Conditionを用いる |
c プログラム作成者： 北條仁士， 筑波大学 2004年 11月 |
c---------------------------------------------------------------------|
implicit real*8(a-h,o-z)
dimension BX(501,201,201), BY(501,201,201), BZ(501,201,201)
dimension EX(501,201,201), EY(501,201,201), EZ(501,201,201)
dimension PX(501,201,201), PY(501,201,201), PZ(501,201,201)
dimension EX1(501,201,201), EY1(501,201,201), EZ1(501,201,201)
dimension DE(501)
common /BLK1/w0,B0,OMEGA
pi = 3.141592654d0
c---------------------------------------------------------------------|
c 数値計算パラメータの設定 |
c w0 = 3.0d11 のとき，長さの規格化 c/w0 = 1 mm , dx=dy=dz= 0.1 mm |
c w0 = 3.0d10 のとき，長さの規格化 c/w0 = 1 cm , dx=dy=dz= 1.0 mm |
c---------------------------------------------------------------------|
w0 = 2.5d10
B0 = 1.0d0
dt = 0.05d0
dx = 0.10d0
dy = 0.10d0
dz = 0.10d0
wx = dt / dx
wy = dt / dy
wz = dt / dz
c---------------------------------------------------------------------|
c システムパラメータの設定 |
c 入力周波数：OMEGA は GHz で表示 |
c システム長：Xmin, Zmin, Xmax, Zmax はミリ単位で表示 |
c---------------------------------------------------------------------|
itmax = 5000
OMEGA = 2.45d0
FREQ = 2.0d9*pi*OMEGA / w0
Xmin = 0.0d0
Ymin = 0.0d0
Zmin = 0.0d0
Xmax = 20.0d0
Ymax = 10.0d0
Zmax = 10.0d0
IMX = idint( (Xmax - Xmin) / dx + 0.5d0 )
JMX = idint( (Ymax - Ymin) / dy + 0.5d0 )
KMX = idint( (Zmax - Zmin) / dz + 0.5d0 )
IMDX = IMX / 100
JMDY = JMX / 100
KMDZ = KMX / 100
c---------------------------------------------------------------------|
c 電磁場の初期値の設定 |
c---------------------------------------------------------------------|
do k=1,KMX+1
do j=1,JMX+1
do 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
EX1(i,j,k) = 0.0d0
EY1(i,j,k) = 0.0d0
EZ1(i,j,k) = 0.0d0
end do
end do
end do
c---------------------------------------------------------------------|
c 計算のスタート:時間 ： time = ( it ) * dt |
c---------------------------------------------------------------------|
it = 0
20 continue
it = it + 1
c---------------------------------------------------------------------|
c z=Zminで励起する入射波の条件設定 |
c---------------------------------------------------------------------|
time = dt * dfloat(it)
do k=1,KMX+1
do j=1,JMX+1
z = dfloat(k-1)*dz + Zmin
y = dfloat(j-1)*dy + Ymin
p0 = 2.d0*pi*OMEGA / 30.d0
C1 = - p0*pi/10.d0 / ( (pi/10.d0)**2 + (pi/10.d0)**2 )
C2 = p0*pi/10.d0 / ( (pi/10.d0)**2 + (pi/10.d0)**2 )
c BX(1,j,k)= dcos(pi*y/10.d0)*dcos(pi*z/10.d0)*dsin(-FREQ*time)
c EY(1,j,k)= C1*dcos(pi*y/10.d0)*dsin(pi*z/10.d0)*dcos(-FREQ*time)
c EZ(1,j,k)= C2*dsin(pi*y/10.d0)*dcos(pi*z/10.d0)*dcos(-FREQ*time)
p0 = 2.d0*pi*OMEGA / 30.d0
C2 = p0*pi/10.d0 / ( (pi/10.d0)**2 + (pi/10.d0)**2 )
EX(1,j,k)= dsin(pi*y/10.d0)*dsin(pi*z/10.d0)*dsin(-FREQ*time)
c BZ(1,j,k)= C2*dcos(pi*y/10.d0)*dsin(pi*z/10.d0)*dcos(-FREQ*time)
end do
end do
c---------------------------------------------------------------------|
c 電磁場の磁場成分 BX, BY, BZ の計算 |
c---------------------------------------------------------------------|
do i=1,IMX+1
do j=1,JMX
do k=1,KMX
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
end do
end do
end do
do i=1,IMX
do j=1,JMX+1
do k=1,KMX
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 )
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
end do
end do
end do
do i=1,IMX
do j=1,JMX
do k=1,KMX+1
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 )
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
end do
end do
end do
c---------------------------------------------------------------------|
c 誘起電流成分 PX, PY, PZ の計算 |
c---------------------------------------------------------------------|
do k=2,KMX
do j=2,JMX
do 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)
PXT = C1*PX(i,j,k) + D1*dt*F*EX(i,j,k)
PXT = PXT - D1*dt*A*( PY(i,j,k)*GZ - PZ(i,j,k)*GY )
PYT = C1*PY(i,j,k) + D1*dt*F*EY(i,j,k)
PYT = PYT - D1*dt*A*( PZ(i,j,k)*GX - PX(i,j,k)*GZ )
PZT = C1*PZ(i,j,k) + D1*dt*F*EZ(i,j,k)
PZT = PZT - D1*dt*A*( PX(i,j,k)*GY - PY(i,j,k)*GX )
PX(i,j,k) = PXT
PY(i,j,k) = PYT
PZ(i,j,k) = PZT
end do
end do
end do
c---------------------------------------------------------------------|
c 電磁場の電場成分 EX, EY, EZ の計算 |
c---------------------------------------------------------------------|
do i=1,IMX
do j=2,JMX
do k=2,KMX
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) )
end do
end do
end do
do i=2,IMX
do j=1,JMX
do k=2,KMX
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)
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) )
end do
end do
end do
do i=2,IMX
do j=2,JMX
do k=1,KMX
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)
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) )
end do
end do
end do
c---------------------------------------------------------------------|
c YZ面（x=Xmax）でのMurの1次の吸収境界条件 |
c---------------------------------------------------------------------|
CT = ( wx - 1.d0 ) / ( wx + 1.d0 )
do j=1,JMX
do k=1,KMX+1
EY(IMX+1,j,k) = EY1(IMX,j,k) + CT*( EY(IMX,j,k)-EY1(IMX+1,j,k) )
end do
end do
do j=1,JMX+1
do k=1,KMX
EZ(IMX+1,j,k) = EZ1(IMX,j,k) + CT*( EZ(IMX,j,k)-EZ1(IMX+1,j,k) )
end do
end do
c---------------------------------------------------------------------|
c YZ面（x=Xmax）での EY=EZ=0 |
c---------------------------------------------------------------------|
c do j=1,JMX
c do k=1,KMX+1
c EY(IMX+1,j,k) = 0.d0
c end do
c end do
c do j=1,JMX+1
c do k=1,KMX
c EZ(IMX+1,j,k) = 0.d0
c end do
c end do
c---------------------------------------------------------------------|
c YZ面（x=Xmin）での EY=EZ=0 |
c---------------------------------------------------------------------|
c do j=1,JMX
c do k=1,KMX+1
c EY(1,j,k) = 0.d0
c end do
c end do
c do j=1,JMX+1
c do k=1,KMX
c EZ(1,j,k) = 0.d0
c end do
c end do
c---------------------------------------------------------------------|
c XZ面（y=Ymax）での EX=EZ=0 |
c---------------------------------------------------------------------|
do i=1,IMX+1
do k=1,KMX
EZ(i,JMX+1,k) = 0.d0
end do
end do
do i=1,IMX
do k=1,KMX+1
EX(i,JMX+1,k) = 0.d0
end do
end do
c---------------------------------------------------------------------|
c XZ面（y=Ymin）での EX=EZ=0 |
c---------------------------------------------------------------------|
do i=1,IMX+1
do k=1,KMX
EZ(i,1,k) = 0.d0
end do
end do
do i=1,IMX
do k=1,KMX+1
EX(i,1,k) = 0.d0
end do
end do
c---------------------------------------------------------------------|
c XY面（z=Zmax）での EX=EY=0 |
c---------------------------------------------------------------------|
do i=1,IMX
do j=1,JMX+1
EX(i,j,KMX+1) = 0.d0
end do
end do
do i=1,IMX+1
do j=1,JMX
EY(i,j,KMX+1) = 0.d0
end do
end do
c---------------------------------------------------------------------|
c XY面（z=Zmin）での EX=EY=0 |
c---------------------------------------------------------------------|
do i=1,IMX
do j=1,JMX+1
EX(i,j,1) = 0.d0
end do
end do
do i=1,IMX+1
do j=1,JMX
EY(i,j,1) = 0.d0
end do
end do
c---------------------------------------------------------------------
do k=1,KMX+1
do j=1,JMX+1
do 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)
end do
end do
end do
c---------------------------------------------------------------------|
c 電磁場の十分に小さい値はゼロで置き換える |
c it = itmax になれば，計算を終了する |
c---------------------------------------------------------------------|
do k=1,KMX+1
do j=1,JMX+1
do 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
end do
end do
end do
IF( it.EQ.500 ) go to 501
IF( it.EQ.1000 ) go to 502
IF( it.EQ.1500 ) go to 503
IF( it.EQ.2000 ) go to 504
IF( it.EQ.2500 ) go to 505
IF( it.EQ.3000 ) go to 506
IF( it.EQ.3500 ) go to 507
IF( it.EQ.4000 ) go to 508
IF( it.EQ.4500 ) go to 509
if( it.eq.itmax ) go to 510
go to 20
c---------------------------------------------------------------------|
c 電磁場の値のデータファイルを作る |
c---------------------------------------------------------------------|
501 open( 11, FILE='exp1', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(11,911) x, EX(i,51,51)
911 format(1h , e15.5, 2x, e15.5 )
end do
close(11)
open( 21, FILE='eyp1', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(21,911) x, EY(i,51,51)
end do
close(21)
go to 20
c--------------------------------------------------
502 open( 12, FILE='exp2', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(12,911) x, EX(i,51,51)
end do
close(12)
open( 22, FILE='eyp2', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(22,911) x, EY(i,51,51)
end do
close(22)
go to 20
c--------------------------------------------------
503 open( 13, FILE='exp3', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(13,911) x, EX(i,51,51)
end do
close(13)
open( 23, FILE='eyp3', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(23,911) x, EY(i,51,51)
end do
close(23)
go to 20
c--------------------------------------------------
504 open( 14, FILE='exp4', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(14,911) x, EX(i,51,51)
end do
close(14)
open( 24, FILE='eyp4', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(24,911) x, EY(i,51,51)
end do
close(24)
go to 20
c--------------------------------------------------
505 open( 15, FILE='exp5', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(15,911) x, EX(i,51,51)
end do
close(15)
open( 25, FILE='eyp5', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(25,911) x, EY(i,51,51)
end do
close(25)
go to 20
c--------------------------------------------------
506 open( 16, FILE='exp6', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(16,911) x, EX(i,51,51)
end do
close(16)
open( 26, FILE='eyp6', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(26,911) x, EY(i,51,51)
end do
close(26)
go to 20
c--------------------------------------------------
507 open( 17, FILE='exp7', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(17,911) x, EX(i,51,51)
end do
close(17)
open( 27, FILE='eyp7', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(27,911) x, EY(i,51,51)
end do
close(27)
go to 20
c--------------------------------------------------
508 open( 18, FILE='exp8', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(18,911) x, EX(i,51,51)
end do
close(18)
open( 28, FILE='eyp8', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(28,911) x, EY(i,51,51)
end do
close(28)
go to 20
c--------------------------------------------------
509 open( 19, FILE='exp9', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(19,911) x, EX(i,51,51)
end do
close(19)
open( 29, FILE='eyp9', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(29,911) x, EY(i,51,51)
end do
close(29)
go to 20
c-------------------------------------------------
510 open( 20, FILE='exp10', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(20,911) x, EX(i,51,51)
end do
close(20)
open( 30, FILE='eyp10', STATUS='unknown' )
do i=1,IMX+1
x = dfloat(i-1)*dx + Xmin
write(30,911) x, EY(i,51,51)
end do
close(30)
open( 31, FILE='exdata10', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(31,931) (EX(131,j,k),k=1,KMX+1,KMDZ)
931 format(1h , e15.5, 100(2x, e15.5) )
end do
close(31)
open( 32, FILE='eydata10', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(32,932) (EY(131,j,k),k=1,KMX+1,KMDZ)
932 format(1h , e15.5, 100(2x, e15.5) )
end do
close(32)
open( 33, FILE='ezdata10', STATUS='unknown' )
do j=1,JMX+1,JMDY
write(33,933) (EZ(131,j,k),k=1,KMX+1,KMDZ)
933 format(1h , e15.5, 100(2x, e15.5) )
end do
close(33)
open( 34, FILE='dexp10', STATUS='unknown' )
do i=1,IMX
x = dfloat(i-1)*dx + Xmin
DE(i) = EX(i,51,51) - EX(i+1,51,51)
write(34,934) x, DE(i)
934 format(1h , e15.5, 2x, e15.5 )
end do
close(34)
open( 35, FILE='SLEX10', STATUS='unknown' )
do i=1,IMX+1,2
write(35,935) (EX(i,51,k),k=1,KMX+1,KMDZ)
935 format(1h , e15.5, 100(2x, e15.5) )
end do
close(35)
open( 36, FILE='SLEY10', STATUS='unknown' )
do i=1,IMX+1,2
write(36,936) (EY(i,51,k),k=1,KMX+1,KMDZ)
936 format(1h , e15.5, 100(2x, e15.5) )
end do
close(36)
open( 37, FILE='SLEZ10', STATUS='unknown' )
do i=1,IMX+1,2
write(37,937) (EZ(i,51,k),k=1,KMX+1,KMDZ)
937 format(1h , e15.5, 100(2x, e15.5) )
end do
close(37)
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
pi = 3.141592654d0
DENTY = 3.0d11
cc DENTY = 6.75d11
BEX0 = 0.0d0
dmp = 0.08d0
dmp = 0.04d0
F0 = 3.181d9 * DENTY / ( w0*w0 )
A = 1.76d11 * B0 / w0
cc IF( x .LT. 12.d0 ) F = 0.0d0
cc IF( x .GE. 12.d0 ) F = F0*(x-12.d0)/8.d0
cc IF( x .GE. 12.d0 ) F = F0
IF( x.LT.8.0d0 ) F = 0.0d0
IF( x.LT.8.0d0 ) go to 999
F = F0 * dtanh( (x - 8.0d0)/2.0d0 )
cc F = F0*( 1.d0 + dtanh( (x-12.d0)/3.d0 ) )/2.d0
F = F *( dtanh((z-2.d0)/1.5d0) - dtanh((z-8.d0)/1.5d0) )/2.d0
F = F *( dtanh((y-2.d0)/1.5d0) - dtanh((y-8.d0)/1.5d0) )/2.d0
999 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
pi = 3.141592654d0
Sr = 1.d0
Si = 0.d0
Qr = 1.d0
IF( (x-7.d0)*(x-8.d0) .LT. 0.d0 ) Sr = 4.0d0
IF( x .EQ. 7.d0 ) Sr = 2.5d0
IF( x .EQ. 8.d0 ) Sr = 2.5d0
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