Large/Infinite results

Tags:
1. Jun 18, 2015

angelos_physik

Hello to all!
I have written this 1D FDTD code but the results I take are very large and after some timestep they get to infinite/NaN... I really don't know what's wrong. The expressions I have used are the exact given solution for the Maxwell equations.

Code (Fortran):
program yee
implicit none

real :: Delta, tau,omega,n,mu,time,Eprev,Hprev
real, dimension(:), allocatable :: sigma,eps,A,B,C,D,J
real, dimension(:), allocatable :: E,H
real, parameter :: pi=3.1415926
integer :: i,m=10000, Length=5000,ls,k

Allocate(E(Length),H(Length),sigma(Length),eps(Length))
Allocate(A(Length),B(Length),C(Length),D(Length))
Allocate(J(Length))

!Initialization
open(unit=1,file="Elektrfeld.txt")
Delta=0.02d0
tau=0.9d0*Delta
omega=2.0d0*pi
n=1.46d0
mu=1
ls=1000
E(:)=0.0d0
H(:)=0.0d0
A(:)=0.0d0
B(:)=0.0d0
C(:)=0.0d0
D(:)=0.0d0
J(:)=0.0d0
Eprev=0.0d0
Hprev=0.0d0

do i=1,Length
if (i.le.6/Delta) then
sigma(i)=1.0d0
else if (i.gt.6/Delta .and. i.lt.((Length*Delta-6.0d0)/Delta))then
sigma(i)=0.0d0
else if (i.ge.(Length*Delta-6.0d0)/Delta .and. i.le.Length) then
sigma(i)=1.0d0
end if
end do

do i=1,Length
if (i.lt.Length*0.5d0) then
eps(i)=1.
else if (i.ge.1100 .and. i.lt.(1100+2)) then
eps(i)=n**2.
else if (i.ge.(1100+2)) then
eps(i)=1.
end if
end do

do i=1,Length
A(i)=(1.-sigma(i)*tau/(2.*mu))/(1.+sigma(i)*tau/(2.*mu))
B(i)=(tau/mu)/(1.+(sigma(i)*tau)/(2.*mu))
C(i)=(1.-sigma(i)*tau/(2.*eps(i)))/(1.+sigma(i)*tau/(2.*eps(i)))
D(i)=(tau/eps(i))/(1.+sigma(i)*tau/(2.*eps(i)))
end do

write(*,*) "time"

!----------------Main calculation----------------
do k=1,m        !time
E(1)=0.0d0
E(Length)=0.0d0
J(ls)=sin(2.*pi*k*tau)*exp(-((k*tau-30.)/10.)**2.)
do i=2,Length
Eprev=E(i)
E(i)=D(i)*(H(i)-H(i-1))/Delta+C(i)*E(i)-D(i)*J(i)
end do
do i=1,Length
Hprev=H(i)
H(i)=B(i)*(E(i)-E(i-1))/Delta+A(i)*Hprev
end do
if (k.eq.time ) then
write(1,*) "-----",k
do i=800,1200
write(1,*) E(i)
end do
end if
end do

close(1)
end program yee

Last edited: Jun 18, 2015
2. Jun 18, 2015

phinds

we have "code" tags that get your code properly formatted. It would be a good idea if you modified your post to use them.

3. Jun 18, 2015

gsal

Well, at quick glance, the second i-loop goes from 1 to Length and it contains the term E(i-1), this means that for i=1, you are calling on E(0) which is not defined...it is outside your program altogether given the way you allocated the arrays, their indexes start at 1.