onyx2.tex

计算光子晶体能带的程序

& Inverse of renormalised permittivity\\
input & complex,pointer & mu\_inv(:,:,:,:,:) 
& Inverse of renormalised permeability\\
input & real,pointer & sigma(:,:,:) & Renormalised electric conductivity\\
input & real,pointer & sigma\_m(:,:,:) & Renormalised magnetic conductivity\\
input & integer & nt & Number of time steps to integrate for\\
input & integer & iblock & Current time step block\\
output & complex,pointer & fft\_data(:,:) & Array for storing the fields\\
input & integer & store\_pts(:,:) & Positions at which we store the fields\\
in/out & complex,pointer & intcurle\_1(:,:,:) & 
Integral of $(\nabla\times\mathbf{E})_{z}$ on left side of system\\
in/out & complex,pointer & intcurlh\_1(:,:,:) & 
Integral of $(\nabla\times\mathbf{H})_{z}$ on left side of system\\
in/out & complex,pointer & intcurle\_2(:,:,:) & 
Integral of $(\nabla\times\mathbf{E})_{z}$ on right side of system\\
in/out & complex,pointer & intcurlh\_2(:,:,:) & 
Integral of $(\nabla\times\mathbf{H})_{z}$ on right side of system\\
input & real,pointer & wz\_1(:) & 
Strength of absorber on left side of system\\
input & real,pointer & wz\_2(:) & 
Strength of absorber on right side of system\\
\hline
internal & integer & ix,iy,iz,it,i\_pts & Loop counters \\
internal & complex & curl(3) & Work space for the curl\\
internal & real & dtcq2 & $(\delta t\ c_0/Q_0)^2$\\
internal & complex,pointer & temp(:,:,:,:) & Temporary pointer\\
\hline
\end{tabular}
\end{center}
The subroutine forms the heart of the finite difference time domain
calculation. It advances the electric and magnetic fields forward in
time by \texttt{nt} time-steps. The equations used to advance the fields 
are the finite difference time dependent Maxwell equations that we derived
 earlier.

\begin{eqnarray*}
\hat{E}_{1}(\mathbf{r},t+\delta t)=\left[1-\hat{\sigma}\right]
\hat{E}_{1}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{11}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{12}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{13}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{E}_{2}(\mathbf{r},t+\delta t)=\left[1-\hat{\sigma}\right]
\hat{E}_{2}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{21}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{22}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{23}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{E}_{3}(\mathbf{r},t+\delta t)=\left[1-\hat{\sigma}\right]
\hat{E}_{3}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{31}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{32}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{33}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{1}(\mathbf{r},t+\delta t)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{1}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{11}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{12}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{13}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{2}(\mathbf{r},t+\delta t)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{2}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{21}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{22}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{23}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{3}(\mathbf{r},t+\delta t)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{3}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{31}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{32}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{33}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}

There are a few other points that should be mentioned about the \texttt{int}
subroutine. Two sets of fields are stored \texttt{e\_cur} and \texttt{h\_cur},
which refer to the fields at the current time step, and 
\texttt{e\_prev} and \texttt{h\_prev} which refer to the previous time step. 
The first thing the routine does at each time step is to switch round the
current and previous fields so that the current fields become the previous 
ones and vice versa. 
This can be done very efficiently because the
arrays in which the fields are stored are accessed as Fortran 90 
\textbf{pointers} so all you have to do is switch over which pointer points 
to which array.
\begin{verbatim}
temp=>e_cur
e_cur=>e_prev
e_prev=>temp

temp=>h_cur
h_cur=>h_prev
h_prev=>temp
\end{verbatim}
This means the arrays \texttt{e\_cur} and 
\texttt{h\_cur} can now be used to hold the updated fields at the next
time step as they are calculated using the equations given earlier.
By doing this we are able to hold on to the fields at two successive
time steps which as we shall see later we need when we want to
calculate the energy density.

Normal Bloch boundary conditions or perfect metal boundary conditions
are provided by the set subroutines \texttt{bc\_xmax\_bloch} 
or \texttt{bc\_xmax\_metal} \emph{etc}.
Notice that the equations given above are only used to update the fields in
the region of the system which is not part of any absorbing layer. Calls are
made to the subroutines \texttt{PML\_e} and \texttt{PML\_h} to update the 
fields in the absorbing layers.

At the end of each time step some components of the fields are stored in
the array \texttt{fft\_data}. The array \texttt{store\_pts} determines which 
components are to be sampled and the parameter \texttt{n\_pts\_store} fixes
the number of fields components to be stored. The array \texttt{store\_pts}
works as follows; for the $i$th point to be stored \texttt{store\_pts(i,2)},
\texttt{store\_pts(i,3)} and \texttt{store\_pts(i,4)} contain the x, y and z
co-ordinates of the mesh point at which the fields are to be sampled.
\texttt{store\_pts(i,1)} determines which of the six fields components is
to be stored; 1 corresponding to $E_x$, 2 is $E_y$, \emph{etc} up to
6 being $H_z$.
The stored field can also be damped by an exponential term corresponding
to adding a small imaginary part to the energy. This comes in handy when
calculating a Green's function as it ensures the sum in the Fourier 
transform converges to give us the causal Green's function.
The size of this damping is set by the run-time parameter \texttt{damping}.
Another technical detail - it turns out that in order to correctly 
calculate the Green's function we must store the current value for any
electric field, but the \emph{previous} value for the magnetic field.
This arises from the fact that we took the forward time difference for the
electric field and the backwards one for the magnetic field.

\section{Subroutines to initialise the calculation}

\subsection{\texttt{subroutine setparam}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{setparam}}\\
\hline
internal & integer & nbits & Bit size of an integer\\
internal & integer & set\_bits & Counter for number of set bits\\
internal & integer & i & Loop counter\\
\hline
\end{tabular}
\end{center}
This subroutine reads in the run-time parameters from an input
file called \texttt{infile.dat} which it assumes has already 
been opened and connected to
unit number \texttt{infile}. The parameters are then written out
to make up the header blocks for the output files, 
\texttt{log.dat}, \texttt{out.dat} and \texttt{fields.dat}. 
If the subroutine can't make
sense of the input file it calls the subroutine \texttt{err}.
The parameters read in from the input file are read into the variables 
in the module \texttt{parameters} and are described in
section~\ref{sec:input}.

\subsection{\texttt{subroutine initfields}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutines \texttt{initfields} and
\texttt{initfields\_dos}}\\
\hline
input & real & g1(3),g2(3),g3(3) & Reciprocal lattice vectors\\
output & complex,pointer & e(:,:,:,:) & Electric field\\
output & complex,pointer & h(:,:,:,:) & Magnetic field\\
input & complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
input & complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
input & integer & ix\_cur,iy\_cur,iz\_cur & x, y, z co-ordinates\\
input & integer & i\_pol & Field component label\\
\hline
internal & integer & ix,iy,iz,jx,jy,jz,i & Loop counters \\
internal & real & G(3) & General reciprocal lattice vector\\
internal & real & k(3) & General wavevector\\
internal & real & v(3) & Arbitrary vector with components $\approx 1$\\
internal & real & k\_plus\_G(3) & $\mathbf{k+G}$\\
internal & complex & vec(3) & Temporary vector\\
internal & complex & h0(3) & Temporary magnetic field\\
internal & complex & expo & $\exp{[i\mathbf{(k+G)\cdot r}]}$\\
internal & complex,pointer & div(:,:,:) & Array to store divergence\\
internal & real & dtcq2 & $(\delta t\ c_0/Q_0)^2$\\
\hline
\end{tabular}
\end{center}
The \texttt{initfields} subroutine defines the initial values of the
electric and magnetic fields at the first time step and returns them in
the arrays \texttt{e} and \texttt{h}. Code is supplied for several possible
initial fields; a Gaussian pulse, a delta function pulse and a sum of
plane waves similar to the starting conditions used by Chan, Yu and
Ho~\cite{ordern}. The delta function is centred on the point 
(\texttt{ix\_cur,iy\_cur,iz\_cur}) and the parameter \texttt{i\_pol} determines
which component of the fields contains the delta function. 
The subroutine also takes as arguments \texttt{eps\_hat} and \texttt{mu\_hat}
in order that tests on the divergence of the initial fields can be made
if required.

A version of this subroutine \texttt{initfields\_dos} is
provided, and should be used when calculating the density of states, since it
sets the starting fields to be the required delta function. In addition,
another version, \texttt{initfields\_trans} is provided for the 
transmission/reflection calculation and sets the initial fields to be some
Gaussian pulse. Notice that it is important for a successful 
transmission/reflection calculation that the initial fields contain only
one mode at each frequency. If the initial fields were to contain more than
one mode then it would become impossible to work out which part of the
scattered field corresponded to which incident mode. It is only by virtue
of the system being linear that we can calculate transmission and reflection
for all frequencies at once. If the system were non-linear then frequency
mixing could occur and it would be impossible to determine which
scattered waves corresponded to which incident frequency.

Whilst it is hoped that the subroutines given here will serve as some sort
of guide, it is clear that, in general, the user will have to write his own
subroutine to set the initial conditions appropriately for the calculation
he has in mind.

\subsection{\texttt{subroutine defcell}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{defcell}}\\
\hline
output & complex,pointer & eps(:,:,:) & Physical permittivity\\
output & complex,pointer & mu(:,:,:) & Physical permeability\\
output & real,pointer & sigma(:,:,:) & Electric conductivity\\
output & real,pointer & sigma\_m(:,:,:) & Magnetic conductivity\\
input & real & u1(3),u2(3),u3(3) & Unit lattice vectors\\
\hline
\end{tabular}
\end{center}
This subroutine defines the physical properties of the unit cell by
setting the values of $\varepsilon$ $\mu$ and the conductivity
$\sigma$ at each point on the discrete mesh. For historical reasons a
magnetic conductivity $\sigma_{m}$ is also defined. When defining a
new unit cell this is the subroutine which must be replaced.

\subsection{\texttt{subroutine init\_store\_pts\_band},
 \texttt{init\_store\_pts\_dos},
\texttt{init\_store\_pts\_trans}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicol