1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Numerical solution for 1d fdtd maxwell equation using yee algorithm

  1. Apr 24, 2012 #1
    1. The problem statement, all variables and given/known data
    to compute 1d fdtd maxwell equation using yee algorithm with fortran 90

    2. Relevant 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)

    3. The attempt at a solution

    The fortran 90 code:
    Code (Text):

    !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)

    do i = iinit,ilast
      do n  = 1,500
         Ez(iinit+1/2,n) = 0
         Hy(iinit+1,n) = 0
      end do
    end do

    Ca = 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
  2. jcsd
  3. Apr 27, 2012 #2
    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.
    Last edited: Apr 27, 2012
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook