kuma4578 の回答履歴

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平面についても同様に書きました。

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平面についても同様に書きました。