- #1

s_hy

- 61

- 0

I wrote programming code for fdtd 2d tfsf in fortran 90 as below:

Code:

```
!simulation for tmz mode incident plane in 2d
subroutine tfsf2d
implicit none
double precision :: f0,miu,eps,delta,S,c,lambda,E0,N0
double precision :: deltat
integer :: minit,pinit
integer :: iinit,ilast
integer :: jinit,jlast
double precision,parameter :: pi = 3.14159265
double precision :: phi,cosphi,sinphi
integer :: h,m
double precision,dimension(300) :: Hm,Em
integer :: y
double precision :: f,G,A,B
double precision,dimension(100) :: k
double precision,dimension(100) :: Vp
double precision :: Ca,Da,Cb,Db
integer :: i,j
integer :: n,p
double precision :: d,dp,dpp
double precision :: Einc,Hinc
double precision,dimension(300,300):: Ezinc,Hxinc,Hyinc
character(len=20) :: filename
double precision :: miu0,eps0
double precision :: deltax, deltay
double precision :: sigma,sigmax,sigmay,sigmat,sigmatx,sigmaty
integer :: iinitp,ilastp
integer :: jinitp,jlastp
double precision,dimension(300,300) :: Cax,Cay,Dax,Day,Cbx,Cby,Dbx,Dby
double precision,dimension(300,300) :: Ez
double precision,dimension(300,300) :: Hy,Hx
!## parameters ##
f0 = 100.0
miu = 1.0
eps = 1.0
delta = 1.e-3
S = 1.0001
c = 1.0
lambda = c/f0
E0 = 1.0D0
N0 = lambda/delta
deltat = (delta/c)/sqrt(2.0)
minit = 1
pinit = 1
phi = pi/4.0
cosphi = cos(phi)
sinphi = sin(phi)
!print *, 'cosphi=',cosphi,'sinphi=',sinphi
miu0 = 1.0
eps0 = 1.0
sigma = 1.0
!print *, 'lambda=',lambda,'deltat=',deltat
deltax = delta
deltay = delta
iinit = 10
ilast = 100
jinit = 10
jlast = 100
iinitp = 1
ilastp = 110
jinitp = 1
jlastp = 110
sigmax = sigma
sigmay = sigma
sigmat = sigma
sigmatx = sigma
sigmaty = sigma
!initialization for 1d grid
do h = 1, 300
Hm(h) = 0
Em(h) = 0
end do
!initialization for incident wave
do p = 1,300
do m = 1,300
Hxinc(m-1,p+2) = 0
Hyinc(m,p+1) = 0
Ezinc(m-1,p+1) = 0
end do
end do
!initialization for field evolution
do p = 1,300
do m = 1,300
Ez(m-1,p+1)=0
Hx(m-1,p+2) = 0
Hy(m,p+1) = 0
end do
end do
!properties coefficient for 2d wave
do y = 1,45
f = 0.5
k(1) = 2*pi
G = (1/f**2)*sin(f*pi/N0)**2
A = (delta*cos(y*pi/180))/2
B = (delta*sin(y*pi/180))/2
k(y+1) = k(y)-(sin(A*k(y))**2+sin(B*k(y))**2)-G/(A*sin(2*A*k(y))+B*sin(2*B*k(y)))
Vp(y) = 2*pi/k(y+1)
Ca = 1.0
Cb = deltat/((Vp(1)/Vp(y))*eps*delta)
Da = 1.0
Db = deltat/((Vp(1)/Vp(y))*miu*delta)
!Print *,'Ca=',Ca,'Cb=',Cb,'Da=',Da,'Db=',Db
end do
!property coefficient
do m = 1,300
do p = 1,300
Cax(m,p) = exp(-sigma*deltat/eps)
Cbx(m,p) = (1-Cax(m,p))/(sigma*deltax)
Cay(m,p) = exp(-sigma*deltat/eps)
Cby(m,p) = (1-Cay(m,p))/(sigma*deltay)
Dax(m,p) = exp(-sigmatx*deltat/miu)
Dbx(m,p) = (1-Dax(m,p))/(sigmatx*deltax)
Day(m,p) = exp(-sigmaty*deltat/miu)
Dby(m,p) = (1-Day(m,p))/(sigmaty*deltay)
end do
end do
!time loop
do n = 1,200
write (filename, "('data',I3.3,'.dat')") n
open (unit=31,file=filename)
!1d background wave
do h = 1, 300
if (Mod(h,2) /= 0) then !j is odd;
Em(h) = Ca*Em(h) + Cb*(Hm(h+1)-Hm(h-1))
!print *,'h=',h,'Em(h)=',Em(h)
end if
end do !h
do h = 1, 300
if (Mod(h,2) == 0) then !j is even;
Hm(h) = Da*Hm(h) + Db*(Em(h+1)-Em(h-1))
!print *,'h=',h,'Hm(h)=',Hm(h)
end if
end do !h
!plane wave source
Em(1) = E0*sin (2*pi*f0*n*deltat)
!incident field
do m = 2,2*ilast+2 !odd
do p = 2,2*jlast+2 !odd
if (Mod(m,2)/=0) then
if (Mod(p,2)/=0) then
d = cosphi*(m - minit) + sinphi*(p - pinit)
dp = d - int(d)
Einc = (1-dp)*Em(int(d))+dp*Em(int(d+1))
Ezinc(m,p+2) = Einc
end if
end if
end do
end do
do m = 2,2*ilast+2
do p = 2,2*jlast+2
if (Mod(m,2) == 0) then !even
if (Mod(p,2) == 0) then !even
d = cosphi*(m - minit) + sinphi*(p - pinit)
dpp = int(d)+0.5
dp = dpp-int(dpp)
Hinc = (1-dp)*Hm(int(dpp-0.5))+dp*Hm(int(dpp+0.5))
Hxinc(m,p+1) = sinphi*Hinc
Hyinc(m-1,p+2) = -cosphi*Hinc
!print *,'m=',m,'p=',p,'Hyinc(m-1,p+2)=',Hyinc(m-1,p+2)
!end if
end if
end if
end do
end do
!##TFSF boundary##
do p = jinit-1,jlast-1
if (Mod(p,2) == 0) then
Hy(10,p)=Hy(10,p)+Dby(10,p)*Ezinc(11,p+1)
!print *,'p=',p,'Hy(10,p)=',Hy(10,p)
end if
end do
do p = jinit-1,jlast-1
if (Mod(p,2) == 0) then
Hy(200,p)=Hy(200,p)+Dby(200,p)*Ezinc(199,p+1)
!print *,'p=',p,'Hy(200),=',Hy(200,p)
end if
end do
do m = iinit-1,ilast-1 !iinit+1,2*ilast-1
if (Mod(m,2) == 0) then
Hx(m,10)=Hx(m,10)+Dbx(m,10)*Ezinc(m+1,11)
!print *,'Hx(m,8)=',Hx(m,8)
end if
end do
do m = iinit-1,ilast-1
if (Mod(m,2) == 0) then
Hx(m,200)=Hx(m,200)-Dbx(m,200)*Ezinc(m+1,199)
!print *,'Hx(m,200)=',Hx(m,200)
end if
end do
do p = jinit-1,jlast-1
if (Mod(p,2)/=0) then
Ez(11,p+1) = Ez(11,p+1)-Cby(11,p)*Hyinc(10,p+1)
!print *,'p=',p,'Ez(11,p+1)=',Ez(11,p+1)
end if
end do
do p = jinit-1,jlast-1
if (Mod(p,2)/=0) then
Ez(199,p+1) = Ez(199,p+1)+Cby(199,p)*Hyinc(201,p+1)
!print *,'p=',p,'Ez(199,p+1)=',Ez(199,p+1)
end if
end do
do m = iinit-1,ilast-1
if (Mod(m,2)/=0) then
Ez(m+1,11) = Ez(m+1,11)+Cbx(m,11)*Hxinc(m+1,9)
!print *,'m=',m,'Ez(m+1,11)=',Ez(m+1,11)
end if
end do
do m = iinit-1,ilast-1
if (Mod(m,2)/=0) then
Ez(m+1,199) = Ez(m+1,199)-Cby(m,199)*Hxinc(m+1,199)
!print *,'m=',m,'Ez(m+1,199)=',Ez(m+1,199)
end if
end do
!## end tfsf boundary##
!## field evolution ##
do m = iinitp,ilastp
do p = jinitp,jlastp
if (Mod(m,2)/=0) then
if (Mod(p,2)/=0) then
Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy(m 1,p+2)+Hx(m,p+1)-Hx(m,p+3))
!print *,'n=',n,'m=',m,'p=',p,'Ez(m,p+2)=',Ez(m,p+2)
!write (31,*) m,p,Ez(m,p+2)
!if (p == 219) write (31,*) ' '
end if
end if
end do
end do
do m = iinitp,ilastp
do p = iinitp,jlastp
if (Mod(m,2) == 0) then
if (Mod(p,2) == 0) then
Hx(m-1,p+2) = Dax(m,p)*Hx(m-1,p+2)+Dbx(m,p)*(Ez(m-2,p+1)-Ez(m-1,p+3))
Hy(m,p+1) = Day(m,p)*Hy(m,p+1)+Dby(m,p)*(Ez(m+1,p+1)-Ez(m-1,p+1))
!print *, 'n=',n,'m=',m,'p=',p,'Hy(m,p+1)=',Hy(m,p+1)
end if
end if
end do
end do
close (unit=31)
end do !n
return
end subroutine
```

i have a problem to call the tfsf boundary condition in field evolution. it's seems that the field evolution didn't react with the boundary (the incident wave = 0)...anyone have any idea which part did i make mistake or other idea how to approach this problem.

thank you