Transcript FDTD and the Plasmon Cluster

Finite Difference Time Domain
(FDTD) Method
February 14, 2014
Yee Cell (3-D)
• Kane Yee (1966) suggested the following cell
for evaluating discretized Maxwell Equations:
Yee Cell (2-D)
Finite Difference Equations
• Maxwell’s Equations are continuous.
• For all but the simplest geometry, closed form
solutions of Maxwell’s Equations are not able
to be obtained.
• We must convert the continuous equations
(with infinitesimal differences) to discrete
equations with finite differences
Discretize and Conquer
((Ey(i+1, j)-Ey(i, j)) – (Ex(i, j+1) – Ex(i, j)))*(1/dx) =
-(mu_0)*(Hz(i, j)now - Hz(i, j)before)/dt
Hz(i,j) = Hz(i,j) – (dt/mu_0)*(Ey(i,j)-Ey(i+1,j) + Ex(i,j+1) - Ex(i,j))*(1/dx)
(Hz(i, j)-Hz(i-1, j))*(1/dx) = (eps_0)*(Ey(i,j)now - Ey(i,j)before) / dt
Ey(i,j) = Ey(i,j) - (dt/eps_0)*(Hz(i-1,j)-Hz(i,j))*(1/dx)
Convergence
• When we go from continuous to discrete, does the
solution become “incorrect”?
• If you keep making the spatial step smaller and smaller,
eventually the results converge to the correct result
– Rule of thumb: make Dx equal to 1/100 of the wavelength
of the incident light
• The Courant Condition specifies a value for Dt:
Programming
• I wrote examples of 1-D and 2-D FDTD in
Fortran 90, C and Python.
• Maxim uses Fortran 90 for everything and
therefore so should you.
– Ask him for the Absoft Fortran 90 editor. If you try
to compile and run with this editor it will give you
trouble. I suggest using the cluster.
Fortran 90
• Variables must be declared (but not initialized) at the very
beginning
• Constants must be declared and initialized at the beginning
• After variables are declared and constants are declared and
initialized, variables must be set
• Then you can start writing your code.
• End your code with the word “end”.
Code Example
implicit none
!
!~~~ Fundamental Constants ~~~!
!
double precision, parameter :: pi=3.141592653589793D0,c=299792458.0D0
double precision, parameter :: mu0=4.0D-7*pi,eps0=1.0D0/(c*c*mu0)
!
!~~~ Grid, Time Steps, Unit Cells ~~~!
!
double precision, parameter :: dx=1.0d-9,dy=1.0d-9
integer, parameter :: Nt=5000
integer, parameter :: Nx=601
!
!~~~ Source ~~~!
!
integer, parameter :: ms=10
!m position of the incident source
double precision, parameter :: E0=1.0D0 !incident pulse peak amplitude
double precision pulse(Nt)
double precision, parameter :: omega=2.0*pi*2.4180D14*5.00 !laser frequency in
eV
!
!~~~ Time Step ~~~!
!
double precision, parameter :: dt=dx/(2.0D0*c)
!
!~~~ Dielectric ~~~!
!
double precision eps_r(Nx)
double precision, parameter :: eps_rel_dielectric=1.0!4.0
double precision, parameter :: eps_rel_vacuum=1.0
integer, parameter :: xnd = 300 !Grid point where dielectric starts
!
!~~~ EM field components ~~~!
!
double precision Ex(Nx),Hz(Nx-1)
!
!~~~ Boundary condition ~~~!
!
double precision ex_low_m2
double precision ex_low_m1
double precision ex_high_m2
double precision ex_high_m1
!
!~~~ Other parameters and variables ~~~!
!
integer n,i
double precision t
double precision, parameter :: dt_eps0=dt/eps0,dt_mu0=dt/mu0
!
!~~~ Fire off the source ~~~!
!
!Sinewave
do n=1,Nt
t=dt*n
pulse(n)=E0*SIN(omega*t)
enddo
!
!~~~ Structure ~~~!
!
do i=1,Nx
if((i<xnd)) then
eps_r(i)=eps_rel_vacuum
else
eps_r(i)=eps_rel_dielectric
endif
enddo
!
!~~~ Initialize remaining variables ~~~!
!
Ex=0.0
Hz=0.0
ex_low_m2=0.0
ex_low_m1=0.0
ex_high_m2=0.0
ex_high_m1=0.0
do n=1,Nt
!Hz Update
do i=1,Nx-1
Hz(i)=Hz(i)-dt_mu0*(Ex(i+1)-Ex(i))/dx
enddo
!Ex Update
do i=2,Nx-1
if(i==ms) then !Excitation
Ex(i)=Ex(i)-(1/eps_r(i))*dt_eps0*(Hz(i)-Hz(i-1))/dx+ &
pulse(n)
else
Ex(i)=Ex(i)-(1/eps_r(i))*dt_eps0*(Hz(i)-Hz(i-1))/dx
endif
enddo
! Inc Boundary Condition
Ex(1) = ex_low_m2
ex_low_m2 = ex_low_m1
ex_low_m1 = Ex(2)
Ex(Nx) = ex_high_m2
ex_high_m2 = ex_high_m1
ex_high_m1 = Ex(Nx-1)
enddo !Nt
Import into Matlab and type Plot(Ex):
!Record Ex
OPEN(file='Ex.dat',unit=32)
do i=1,Nx
write(32,*) Ex(i)
enddo
CLOSE(unit=32)
end
Results
• As one is learning FDTD,
it is helpful to plot
electric field over the
entire grid
Results
• Graphs of the T and R
coefficients (obtained
from the Poynting
vector) are measured in
experiments
Next Steps
• Look at a single-processor 2-D code with a
dielectric cylinder (ask me or Maxim for it)
• The absorbing boundary conditions (ABC’s)
are implemented via convolutional perfectly
matched layers (CPML)
– It’s there and it works. In general, there’s no need
to mess with it.
Later On
• The best metals for SPP’s are silver and gold
– Metal is represented by the Drude Model
• The interaction of the two-level emitters with the classical
EM field is handled by the Liouville-von Neumann equation
• Longer simulations benefit from message passing interface
(MPI); the simulation is split up between several processors
– Ask me or Maxim for the MPI slits code
• Rather than using point sources, we can define a Total Field
/ Scattered Field boundary to introduce a plane wave into
the system