- #1

angelos_physik

- 1

- 0

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.

Thank you for your help.

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.

Thank you for your help.

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"
read(*,*) 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: