FDTD 遠方界のプログラム
FDTD法で遠方界を計算するプログラムを作ったのですがコンパイルはうまくいっても実行結果がちゃんと表示されません。プログラムに問題があるのか、大きな値を与えるとbus errorになってしまいます。どなたか遠方界のプリグラムに詳しい方がいらっしゃったら教えて下さい。よろしくお願いします。
プログラム内容
ubroutine far_field
implicit none
call far_field_yz_plane
call far_field_zx_plane
call far_field_xy_plane
return
end subroutine
!===================================================================================
! y-z 平面に対して
!===================================================================================
subroutine far_field_yz_plane
use consts
use fdtd
implicit none
integer :: i, j, k,ie
real(8) :: plane_yz, fac_yz
plane_yz = dy*dz
fac_yz = plane_yz/(4.0d0*pi*c*dt)
ie = nx-6
i = ie+1
do j = 5, ny-5
do k = 5, nz-5
x = (i-0.5d0-ic)*dx
y = (j-jc)*dy
z = (k-kc)*dz
eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1))
ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k))
hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k))
hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1))
do l = 1, nang
time_e = time-dt
time_h = time-dt/2.0d0
timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c
tc_e = time_e+timesh
tc_h = time_h+timesh
m_e = int(tc_e/dt+0.5d0)
m_h = int(tc_h/dt+0.5d0)
a_e = (0.5d0+tc_e/dt-m_e)*fac_yz
b_e = (0.5d0-tc_e/dt+m_e)*fac_yz
ab_e = a_e-b_e
a_h = (0.5d0+tc_h/dt-m_h)*fac_yz
b_h = (0.5d0-tc_h/dt+m_h)*fac_yz
ab_h = a_h-b_h
wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h
wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h
uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e
uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e
wy(l,m_h) = wy(l,m_h)-hzs*ab_h
wz(l,m_h) = wz(l,m_h)+hys*ab_h
uy(l,m_e) = uy(l,m_e)+ezs*ab_e
uz(l,m_e) = uz(l,m_e)-eys*ab_e
wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h
wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h
uy(l,m_e+1) = uy(l,m_e+1)-ezs*a_e
uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e
enddo
enddo
enddo
ie = 5
i = ie
do j = 5, ny-5
do k = 5, nz-5
x = (i-0.5d0-ic)*dx
y = (j-jc)*dy
z = (k-kc)*dz
eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1))
ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k))
hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k))
hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1))
do l = 1, nang
time_e = time-dt
time_h = time-dt/2.0d0
timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c
tc_e = time_e+timesh
tc_h = time_h+timesh
m_e = int(tc_e/dt+0.5d0)
m_h = int(tc_h/dt+0.5d0)
a_e = (0.5d0+tc_e/dt-m_e)*fac_yz
b_e = (0.5d0-tc_e/dt+m_e)*fac_yz
ab_e = a_e-b_e
a_h = (0.5d0+tc_h/dt-m_h)*fac_yz
b_h = (0.5d0-tc_h/dt+m_h)*fac_yz
ab_h =a_h-b_h
wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h
wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h
uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e
uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e
wy(l,m_h) = wy(l,m_h)-hzs*ab_h
wz(l,m_h) = wz(l,m_h)+hys*ab_h
uy(l,m_e) = uy(l,m_e)+ezs*ab_e
uz(l,m_e) = uz(l,m_e)-eys*ab_e
wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h
wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h
uy(l,m_h+1) = uy(l,m_e+1)-ezs*a_e
uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e
enddo
enddo
enddo
return
end subroutine
zx平面、xy平面についても同様に書きました。
