- #1

s_hy

- 61

- 0

attached here is my code for 2d fdtd in polar coordinates, from 'numerical electromagnetic: the fdtd method (umran s inan, pg 94-96) written in fortran90. I have try a few approach I could think about to troubleshoot this code but the output is still infinity. Anybody here can give me advice. Thank you

Code:

```
!2d polar coordinates
!reference: Numerical electromagnetic: fdtd method (Umran S Inan, R Marshall)
subroutine test3
implicit none
double precision :: f0,miu0,eps0
double precision :: delta, deltat,delr,delphi
double precision :: c0
integer :: i,j
integer :: n,M
double precision,dimension(600,600):: Ez,Hp,Hr !Ez,Hy,Hx
double precision, parameter :: pi = 3.14159265
character(len=20) :: filename
f0 = 1.e6
c0 = 3.e8
miu0 = 4.0*pi*1.0e-7
eps0 = 8.854e-12
! deltat = 1/sqrt(2.0)*(delta/c0)
delphi = 0.0314 !radian
delr = 15
!delt = 1/sqrt(2.0)*delr*delphi/c0
deltat = delr*delphi/(c0*2)
!print*,'delt=',delt
M = 200
!print*, 'see me?'
!initialize Ez, Hx,Hy to zero at t=0
do i = 1,M
do j = 1,M
Hp(i,j) = 0
Hr(i,j) = 0
Ez(i,j) = 0
end do
end do
do n = 1,500
write (filename, "('data',I3.3,'.dat')") n
open (unit=130,file=filename)
!initiate sinusoidal wavepulse at center
!Ez(1,1) = sin((n+1)*omega*deltat)
Ez(100,100) = Ez(100,100) + sin(2*pi*f0*n*deltat)
!print *, 'Ez(100,100)=',Ez(100,100)
!update electric field equation
do i = 2,M
do j = 2,M
Ez(i,j) = Ez(i,j)+(deltat/eps0)*((delr*(i)*Hp(i,j)-delr*(i-1)*Hp(i-1,j))/(delr*(i)*delr))-(Hr(i,j)-Hr(i,j-1)/(delr*(i)*delphi))
write (130,*) i,j,Ez(i,j)
if (j == 199) write (130,*) ' '
!print *, 'i=',i,'j=',j,'Ez(i,j)=',Ez(i,j)
end do
end do
!update magnetic field equation
do i = 1,M
do j = 1,M
Hr(i,j) = Hr(i,j) - (deltat/miu0)*((Ez(i,j+1) - Ez(i,j))/delr*(i)*delphi)
print *,'i=',i,'j=',j,'Hr(i,j)=',Hr(i,j)
end do
end do
do i = 1,M
do j = 1,M
Hp(i,j) = Hp(i,j) + (deltat/miu0)*((-Ez(i+1,j)+Ez(i,j))/delr)
!print *,'i=',i,'j=',j,'Hp(i,j)=',Hp(i,j)
end do
end do
close (unit=130)
end do !n
```