!**************************************************************************** ! 初期化 !**************************************************************************** subroutine init use consts use fdtd implicit none integer :: id print *,'**** init ****' write(1,*) '**** init ****' ! ******** 波源情報 ******** freq=4.0d9 ! 周波数 [Hz] ! ******** 格子数 ******** nx=20 ny=100 nz=10 ! ******** タイムステップ、格子間隔 ******** ! Courant の安定条件 ! v dt < 1/sqrt((1/dx)**2+(1/dy)**2+(1/dz)**2) ! v dt < dx/sqrt(3) when dx=dy=dz dx=29.1d-3/10.0d0 dy=dx dz=dx dt=1.0d0/freq/100.0d0 if(dt.gt.1.0d0/sqrt((1.0d0/dx)**2+(1.0d0/dy)**2+(1.0d0/dz)**2)/c) then print *,'Sub(init): Courant stability condition must be satisfied' print *,'dt = ',dt print *,'dt < ',1.0d0/sqrt((1.0d0/dx)**2+(1.0d0/dy)**2+(1.0d0/dz)**2)/c stop end if ntime=500 ! 計算する時間ステップ数 write(1,*) 'dt=',dt write(1,*) 'ntime=',ntime ! ******** 媒質定数設定 ******** nmedia=2 ! 1: 真空 eps(1)=epsilon0 mu(1)=mu0 sig(1)=0.0d0 ! 2: PEC,PMC (これは別のルーチンで処理する) ! ******** フィールド更新係数 ******** ! 1: 真空 cex0=1.0d0 cey0=1.0d0 cez0=1.0d0 cexry0=(dt/epsilon0)/dy cexrz0=(dt/epsilon0)/dz ceyrz0=(dt/epsilon0)/dz ceyrx0=(dt/epsilon0)/dx cezrx0=(dt/epsilon0)/dx cezry0=(dt/epsilon0)/dy chxry0=(dt/mu0)/dy chxrz0=(dt/mu0)/dz chyrz0=(dt/mu0)/dz chyrx0=(dt/mu0)/dx chzrx0=(dt/mu0)/dx chzry0=(dt/mu0)/dy ! 3以上: 一般の媒質 do id=3,nmedia cex(id)=(1.0d0-((sig(id)*dt)/(2.0d0*eps(id)))) & /(1.0d0+((sig(id)*dt)/(2.0d0*eps(id)))) cey(id)=(1.0d0-((sig(id)*dt)/(2.0d0*eps(id)))) & /(1.0d0+((sig(id)*dt)/(2.0d0*eps(id)))) cez(id)=(1.0d0-((sig(id)*dt)/(2.0d0*eps(id)))) & /(1.0d0+((sig(id)*dt)/(2.0d0*eps(id)))) cexry(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dy cexrz(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dz ceyrz(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dz ceyrx(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dx cezrx(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dx cezry(id)=(dt/eps(id))/(1.0d0+((sig(id)*dt)/(2.0d0*eps(id))))/dy chxry(id)=(dt/mu(id))/dy chxrz(id)=(dt/mu(id))/dz chyrz(id)=(dt/mu(id))/dz chyrx(id)=(dt/mu(id))/dx chzrx(id)=(dt/mu(id))/dx chzry(id)=(dt/mu(id))/dy end do ! ******** 物体のモデリング ******** call modeling return end subroutine !---------------------------------------------------------------------------- ! 物体のモデリング !---------------------------------------------------------------------------- subroutine modeling use fdtd implicit none integer :: i,j,k ! すべて真空 (=1) に初期化する do i=1,nx do j=1,ny do k=1,nz media_id(i,j,k)=1 end do end do end do ! 金属壁 PEC (=2) の設定 ! 上下の壁 do i=1,nx do j=1,ny media_id(i,j,1)=2 media_id(i,j,nz)=2 end do end do ! 左右の壁 do j=1,ny do k=1,nz media_id(1,j,k)=2 media_id(nx,j,k)=2 end do end do ! 手前の壁 ! do i=1,nx ! do k=1,nz ! media_id(i,1,k)=2 ! end do ! end do return end subroutine ! ! End of file !