onyx.tex

order发计算光子晶体的能带结构

$k_{z}$. Remember that \texttt{akx} \emph{etc.} are defined such
\texttt{akx=(Q(1)*ixmax)*kx}. If one is calculating
a Green's function there is the possibility to loop over \texttt{ipol}, the
variable which determines which field component contains the initial delta
function.
Depending on whether one is performing a band structure calculation or
a density of states calculation separate versions of some subroutines are
provided. Whilst this does result in a certain amount of code duplication
it does mean that switching from one type can be done entirely by changing
subroutine calls in the main program and removes the error prone process
of searching through the whole program trying to find all the places where
changes must be made.

\subsection{\texttt{subroutine driver}}

\begin{center}

\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{driver}}\\
\hline
input & complex,pointer & e\_cur(:,:,:,:) 
& Electric Field at the current time\\
input & complex,pointer & h\_cur(:,:,:,:) 
& Magnetic Field at the current time\\
input & complex,pointer & eps\_inv(:,:,:,:,:) 
& Inverse of renormalised permittivity\\
input & complex,pointer & mu\_inv(:,:,:,:,:) 
& Inverse of renormalised permeability\\
input & complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
input & complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
input & real,pointer & sigma(:,:,:) & Renormalised electric conductivity\\
input & real,pointer & sigma\_m(:,:,:) & Renormalised magnetic conductivity\\
output & complex,pointer & fft\_data(:,:) & Array for storing the fields\\
input & integer,pointer & store\_pts(:,:) 
& Positions at which we store the fields\\
\hline
internal & complex,pointer & e\_prev(:,:,:,:) 
& Electric Field at the previous time\\
internal & complex,pointer & h\_prev(:,:,:,:) 
& Electric Field at the previous time\\
internal & complex,pointer & div(:,:,:) & Divergence of the field\\
internal & real,pointer & rho(:,:,:) & Energy density\\
internal & integer & i & Loop counter\\
internal & integer & iblock & Loop counter\\
\hline
\end{tabular}
\end{center}
This subroutine handles the main part of the calculation which
steps the electromagnetic fields forward in time. The total
number of time steps \texttt{itmax} is divided up into \texttt{n\_blocks}
blocks each of \texttt{block\_size} time steps each (\texttt{n\_blocks}
and \texttt{block\_size} are parameters set from the input file). The
subroutine \texttt{int} performs the actual time integration for the
\texttt{block\_size} time steps in each block and at the end of each
block there is the option to perform several tests on the fields,
charge conservation, energy conservation in order to make sure 
nothing is going wrong.

The routine also does a few other odds and ends; it \textbf{allocates}
storage to the arrays to store the fields at the previous time-step,
and various arrays used  to store the divergence of the fields and the
energy density
when testing the calculation. These arrays are defined here because they
are only required during the time integration loop and are not needed 
by the rest of the program.
 
\subsection{\texttt{subroutine int}}

\begin{center}

\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{int}}\\
\hline
in/out & complex,pointer & e\_cur(:,:,:,:) 
& Electric Field at the current time\\
in/out & complex,pointer & e\_prev(:,:,:,:) 
& Electric Field at the previous time\\
in/out & complex,pointer & h\_cur(:,:,:,:) 
& Magnetic Field at the current time\\
in/out & complex,pointer & h\_prev(:,:,:,:) 
& Magnetic Field at the previous time\\
input & complex,pointer & eps\_inv(:,:,:,:,:) 
& 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\\
\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}.

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