!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
c = 1
E0 = 1
deltat = (delta/c)/sqrt(2.0)
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
 Ez(1,p-1) = 0
 Ez(m-1,p-1) = 0

 do m = 1,2*ilast
 Hx(m+1,p) = 0
 Hy(2*m,p) = 0
 Ez(m+1,p-1) = 0
 Ez(m+1,1) = 0
 Ez(m+1,plast+2) = 0
 end do
end do

!medium's properties

  Ca(m,p) = 1.0
  Cb(m,p) = deltat/(1-0.5*deltax)
  !Print *,'Ca, Cb=',Ca, ' ',Cb

  Da(m,p) = 1.0
  Db(m,p) = deltat/deltax
  !Print *, 'Da, Db=',Db, ' ',Db
  !print*, 'see me?1' 

!initiate sinusoidal wavepulse at center
do n = 1,100
 do m = 1,2*ilast
  do p = 1,2*jlast

   Ez(101,101) = sin((n+1)*omega*deltat)

  end do
 end do
end do

do n = 1,200
   
  write (filename, "('ez',I3.3,'.dat')") n
  open (unit=130,file=filename)
 
  do m = 1,2*ilast-1
    do p = 1,2*jlast-1
    
      if (Mod(m,2)/=0) then
	if (Mod(p,2)/=0) then
	  Ez(m,p) = Ca*Ez(m,p-1)+Cb*(Hy(m+1,p-1)-Hy(i+1,p-1)+Hx(m+1,p+1)-Hx(m,p-1))
        end if
      end if
    end do
  end do

  do m = 1,2*ilast-1
    do p = 1,2*jlast-1

      if (Mod(m,2) == 0) then
	if (Mod(p,2) == 0) then
	  Hx(m-1,p+2) = Da*Hx(m-1,p+2)+Db*(Ez(m-1,p+2)-Ez(m-1,p+2))
          Hy(m,p+1) = Da*Hy(m,p+1)+Db*(Ez(m+1,p+1)-Ez(m-1,p+1))
        end if
      end if
     end do
  end do

 ! set boundary condition
 do m = 1,2*ilast-1
  do p = 1, 2*jlast-1
  if (p == 0) then
    Ez(m,p) = 0
   if (p == plast) then 
     Ez(1,p-1) = 0
     Ez(mlast+1,p+1) = 0
   end if
  end if
  end do
 end do

 close (unit=130)

end do !n

end SUBROUTINE fd2d

subroutines.f90:120.15:

  Ca(m,p) = 1.0
               1
Error: Unexpected STATEMENT FUNCTION statement at (1)
subroutines.f90:121.33:

  Cb(m,p) = deltat/(1-0.5*deltax)
                                 1
Error: Unexpected STATEMENT FUNCTION statement at (1)
subroutines.f90:124.15:

  Da(m,p) = 1.0
               1
Error: Unexpected STATEMENT FUNCTION statement at (1)
subroutines.f90:125.25:

  Db(m,p) = deltat/deltax
                         1
Error: Unexpected STATEMENT FUNCTION statement at (1)