Numerical solution for 1d fdtd maxwell equation using yee algorithm

[s_hy]

Homework Statement

to compute 1d fdtd maxwell equation using yee algorithm with fortran 90

Homework Equations

1D discretization for maxwell equation (TEM mode) :

electric field vector:
Ez(i-1/2,n+1/2) = Ca*(Ez(i-1/2,n-1/2) + Cb(Hy(i,n)-Hy(i-1,n)

magnetic field
Hy(i,n+1) = Da*(Hy(i,n) + Db(Ez(i+1/2,n+1/2)-Ez(i-1/2,n+1/2)

with source (gaussian pulse)
Ez(1/2+(ilast-iinit)/2,n) = E0*sin (2*pi*f0*n*tdelta)

The Attempt at a Solution

The fortran 90 code:

!1d fdtd Simulation in free space 
subroutine fd1d01(f0,miu,delta,S,E0)
implicit none

real        :: f0 
real        :: miu
real        :: delta 
real        :: S 
real        :: E0
integer     :: iinit 
integer     :: ilast 
real        :: Ca
real        :: Da
integer     :: i
integer     :: n
real        :: tdelta 
real        :: c
real,dimension(:,:),allocatable   :: Ez
real,dimension(:,:),allocatable   :: Hy
real, parameter :: pi = 3.14159265
real        :: Cb
real        :: Db
real        :: lambda
real        :: alpha
character(len=20)  ::filename

allocate (Ez(-1:103,-1:503))
allocate (Hy(-1:103,-1:503))

f0 = 1.0
miu = 1.0
delta = 1.0 
S = 1.0
E0 = 1.0
iinit = 0 
ilast = 100
 
c = 3.e8
lambda = c/f0

alpha = 0.04*lambda

tdelta = 1.0*alpha/(S*c)!initialization
do i = iinit,ilast
  do n  = 1,500
     Ez(iinit+1/2,n) = 0
     Hy(iinit+1,n) = 0
  end do
end doCa = 1.0
Cb = tdelta/(delta*alpha)

Da = 1.0
Db = tdelta/(miu*alpha)do n = 1,500
 write (filename, "('ez',I3.3,'.dat')") n
 open (unit=130,file=filename)

  do i = iinit+1,ilast
  Ez(i-1/2,n+1/2) = Ca*(Ez(i-1/2,n-1/2)) + Cb*(Hy(i,n)-Hy(i-1,n)) 
  Hy(i,n+1) = Da*(Hy(i,n)) + Db*(Ez(i+1/2,n+1/2)-Ez(i-1/2,n+1/2)) 
!  Print*, 'Ez(i-1/2,n+1/2)=',Ez(i-1/2,n+1/2)
  Write (130,*) i-1./2, Ez(i-1/2,n+1/2)
  end do !i

  !plane wave source

  Ez(1/2+(ilast-iinit)/2,n) = E0*sin (2*pi*f0*n*tdelta)

 close (unit=130)

end do !n

end SUBROUTINE fd1d01

the problem i am facing right now is the output .dat file for Ez are all 0's. I had try few ways to find the problem source but unable to find out what is wrong with my code. Can anyone here that expert electromagnetic as well as fortran help me.

thank you very much

 

[LawrenceC]

While I am not an expert in electromagnetics, I have spent a great deal of time over the years writing code in the finite element arena. In your case, the first thing I would do is add a write statement on some preliminary calculation to assist ferreting out why you are winding up with zeros for results. Move the write statement along towards the end of the program, recompile, and execute again. In that way you can quickly determine where something may not be defined.

The worst error I ever had was one where a large program of mine would run with the extra write statement included but would get wrong answers when the statement was not present. It turned out I was exceeding an array specification. The addition of the write statements changed how storage was allocated so something that was getting overwritten had already been used in the case where the write statements were included.