- #1
s_hy
- 61
- 0
hi all,
i have wrote codes for 2d fdtd in different permittivity (epsilon). in this code, cell size is 200 x 200, start with eps=1 from center, and different permittivity started at boundary (i1,i2) = (25,75) = (j1,j2), epsilon = 2. the problem is, when the wave propagates and approach (i1,i2) = (25,75) = (j1,j2), i have detected some reflections happened to the wave at timesteps, n after 50, the reflections become obvious...i just wondering, why is that happened. i attached here the animation of the wave and my codes. i really have no idea how to solve this problem.
i have wrote codes for 2d fdtd in different permittivity (epsilon). in this code, cell size is 200 x 200, start with eps=1 from center, and different permittivity started at boundary (i1,i2) = (25,75) = (j1,j2), epsilon = 2. the problem is, when the wave propagates and approach (i1,i2) = (25,75) = (j1,j2), i have detected some reflections happened to the wave at timesteps, n after 50, the reflections become obvious...i just wondering, why is that happened. i attached here the animation of the wave and my codes. i really have no idea how to solve this problem.
Code:
subroutine test7
implicit none
double precision :: f0,A0
double precision :: delta, deltat
double precision :: eps0,epsr,miu0,miur
double precision :: lamda
double precision :: omega
double precision :: S,k,c
integer :: i,j
integer :: ie,je,ic,jc,i1,i2,j1,j2
integer :: n
double precision,dimension(202,202) :: Ez,Hy,Hx
double precision,dimension(202,202) :: Ca,Cb,Da,Db,Dax,Dbx,Day,Dby
!double precision :: Ca,Cb,Da,Db
double precision,dimension(202,202) :: miu,eps,sigma, sigmat
double precision, parameter :: pi = 3.14159265
character(len=20) :: filename
! parameters
f0 = 3.0e7
eps0 = 8.8e-12
miu0 = 4*pi*1e-7
c = 3.e8
lamda = c/f0
omega = 2*pi*f0
print *, 'lamda=',lamda
A0 = 5.0
k = 2*pi / lamda
omega = c*k
!courant stability factor
!S = 1/sqrt(2.0)
! spatial and time grid steps
delta = lamda/10
deltat = delta/(2*c)
print *, 'deltat=',deltat
!model size
ie = 100
je = 100
!source location
ic = ie/2
jc = je/2
!boundary for different permittivity
i1 = ie/4
i2 = ie*3/4
j1 = je/4
j2 = je*3/4
print*, 'see me?'
!initialize Ez, Hx,Hy to zero at t=0
do i = 1,ie
do j = 1,je
Hx(i,j) = 0.0
Hy(i,j) = 0.0
Ez(i,j) = 0.0
end do
end do
!medium's properties
miur = 1.0
epsr = 1.0
do i = 1,ie
do j = 1,je
!initialize permittivity and permeability
eps(i,j) = epsr*eps0
miu(i,j) = miur*miu0
end do
end do
do i = i1,i2
do j = j1,j2
eps(i,j) = 2*eps0
miu(i,j) = 1*miu0
end do
end do
!medium's properties
! initialize electric conductivity & magnetic conductivity
do i = 1,ie
do j = 1,je
sigma(i,j) = 0.0
sigmat(i,j) = 0.0
Ca(i,j) = (2*eps(i,j)-deltat*sigma(i,j))/(2*eps(i,j)+deltat*sigma(i,j))
Cb(i,j) = (2*deltat)/(2*eps(i,j)+deltat*sigma(i,j))
Da(i,j) = (2*miu(i,j)-deltat*sigmat(i,j))/(2*miu(i,j)+deltat*sigmat(i,j))
Db(i,j) = (2*deltat)/(2*miu(i,j)+deltat*sigmat(i,j))
end do
end do
! beginning of timesteps
do n = 1,300
write (filename, "('data',I3.3,'.dat')") n
open (unit=130,file=filename)
!initiate sinusoidal wavepulse at center
Ez(ic,jc) = A0*sin(2*pi*f0*n*deltat)
! calculate ez-field
do i = 2,ie-1
do j = 2,je-1
Ez(i,j) = Ca(i,j)*Ez(i,j) + Cb(i,j)*(Hy(i,j) - Hy(i-1,j) - Hx(i,j) + Hx(i,j-1))
write (130,*) i,j,Ez(i,j)
if (j == 99) write (130,*) ' '
! !print *,'i=',i,'j=',j,'Ez(i,j)=',Ez(i,j)
end do
end do
do i = 1,ie-1
do j = 1,je-1
Hx(i,j) = Da(i,j)*Hx(i,j) + Db(i,j)*(Ez(i,j) - Ez(i,j+1))
!print *,'i=',i,'j=',j,'Hx(i,j)=',Hx(i,j)
end do
end do
do i = 1,ie-1
do j = 1,je-1
Hy(i,j) = Da(i,j)*Hy(i,j) + Db(i,j)*(Ez(i+1,j) - Ez(i,j))
!print *,'i=',i,'j=',j,'Hy(i,j)=',Hy(i,j)
end do
end do
close (unit=130)
end do !n
end SUBROUTINE test7