- #1

s_hy

- 61

- 0

Currently I am doing 2d fdtd (TMz mode) modeling. I expect to have output for Ez, Hx and Hy in one file, which is mean that the output will be in 3 column. But I don't have any idea how format it. Anyone can help me? Below is the code

Code:

```
!problem 3.9, Taflove chapter 3
!to model 2d TMz mode fdtd Simulation using Yee grid.
!time step deltat = (1/sqrt(2)*delta/c) is used
subroutine fd2d
implicit none
double precision :: f0 !frequency
double precision :: miu
double precision :: miur !permittivity
double precision :: miu0
double precision :: E0
double precision :: delta
double precision :: deltat
double precision :: deltax
double precision :: deltay
double precision :: sigma
double precision :: sigmastar
double precision :: Eps
double precision :: Eps0
double precision :: Epsr
double precision :: lambda
double precision :: lambdax
double precision :: lambday
double precision :: nlambdax
double precision :: nlambday
double precision :: omega
double precision :: S
double precision :: k
double precision :: c
integer :: i,j
integer :: m,p
integer :: mlast,plast
integer :: n
integer :: iinit
integer :: ilast
integer :: jinit
integer :: jlast
double precision,dimension(202,202):: Ez
double precision,dimension(202,202):: Hy
double precision,dimension(202,202):: Hx
double precision :: EzAmp
double precision :: HyAmp
double precision :: HxAmp
double precision :: Ca
double precision :: Cb
double precision :: Da
double precision :: Db
double precision, parameter :: pi = 3.14159265
character(len=20) :: filename
f0 = 100.0
miu0 = 1.0
miur = 1.0
Eps0 = 1.0
Epsr = 1.0
!delta = 1.e-3
delta = 0.1
!lambda = c/f0
lambda = 1
c = 1.0
E0 = 1.0
deltat = 1/sqrt(2.0)*(delta/c)
k = 2*pi / lambda
omega = c*k
nlambdax = 5
nlambday = 5
lambdax = lambda
lambday = lambda
deltax = delta
deltay = delta
S = 1.0001
iinit = 1
ilast = 100
jinit = 1
jlast = 100
mlast = 2*ilast-1
plast = 2*jlast-1
print*, 'see me?'
!initialize Ez, Hx,Hy to zero at t=0
do p = 1,2*jlast
do m = 1,2*ilast
if (Mod(m,2) == 0) then
if (Mod(p,2) == 0) then
Hx(m-1,p+2) = 0
Hy(m,p+1) = 0
else
Ez(m-1,p+1) = 0
end if
end if
end do
end do
print*, 'see me?2'
!medium's properties
Ca = 1.0
Cb = deltat/(1-0.5*deltax)
Print *,'Ca=',Ca,'Cb=',Cb
Da = 1.0
Db = deltat/deltax
Print *, 'Da=',Da, 'Db=',Db
!initiate sinusoidal wavepulse at center
do n = 1,200
write (filename, "('ez',I3.3,'.dat')") n
open (unit=130,file=filename)
!initiate sinusoidal wavepulse at center
Ez(101,101) = sin((n+1)*omega*deltat)
do m = 2,2*ilast-1
do p = 2,2*jlast-1
if (Mod(m,2)/=0) then
if (Mod(p,2)/=0) then
Ez(m,p+2) = Ca*Ez(m,p+2)+Cb*(Hy(m+1,p+2)-Hy(m-1,p+2)+Hx(m,p+1)-Hx(m,p+3))
write (130,1) Ez(m,p+2)
!print *, Ez(m,p+2)
end if
end if
end do
end do
do m = 2,2*ilast-1
do p = 2,2*jlast-1
if (Mod(m,2) == 0) then
if (Mod(p,2) == 0) then
Hy(m,p+1) = Da*Hy(m,p+1)+Db*(Ez(m+1,p+1)-Ez(m-1,p+1))
Hx(m-1,p+2) = Da*Hx(m-1,p+2)+Db*(Ez(m-1,p+1)-Ez(m-1,p+3))
write (130,1) Hx(m-1,p+2)
write (130,1) Hy(m,p+1)
end if
end if
end do
end do
1 format (130(3E15.6))
close (unit=130)
end do !n
end SUBROUTINE fd2d
```

Anyway, format I used to write file is not working.

Thank you