onyx2.tex

计算光子晶体能带的程序

\right\}
\]
\[
\hat{E_{y}}(\mathbf{r},t+\delta t)=\frac{1}{1+\omega_{z} \delta t}
\left\{
\hat{E_{y}}(\mathbf{r},t)+\frac{1}{\tilde{\varepsilon}^{22}}
\left[\nabla_{q}\times\hat{\mathbf{H}}'(\mathbf{r},t)
\right]_{y}
\right\}
\]
\[
\hat{E_{z}}(\mathbf{r},t+\delta t)=\hat{E_{z}}(\mathbf{r},t)
+\frac{1}{\tilde{\varepsilon}^{33}}
\left[\nabla_{q}\times\hat{\mathbf{H}}'(\mathbf{r},t)\right]_{z}
+\frac{\omega_{z} \delta t}{\tilde{\varepsilon}^{33}} \sum_{n=0}^{\infty}
\left[\nabla_{q}\times\hat{\mathbf{H}}'(\mathbf{r},t-n\delta t)\right]_{z}
\]
\[
\hat{H}'_{x}(\mathbf{r},t)=\frac{1}{1+\omega_{z} \delta t}
\left\{\hat{H}'_{x}(\mathbf{r},t-\delta t)
-\frac{1}{\tilde{\mu}^{11}}\left(\frac{c_{0}\delta t}{Q_{0}}\right)^{2}
\left[\nabla_{q}\times\hat{\mathbf{E}}(\mathbf{r},t)\right]_{x}
\right\}
\]
\[
\hat{H}'_{y}(\mathbf{r},t)=\frac{1}{1+\omega_{z} \delta t}
\left\{\hat{H}'_{y}(\mathbf{r},t-\delta t)
-\frac{1}{\tilde{\mu}^{22}}\left(\frac{c_{0}\delta t}{Q_{0}}\right)^{2}
\left[\nabla_{q}\times\hat{\mathbf{E}}(\mathbf{r},t)\right]_{y}
\right\}
\]
\[
\hat{H}'_{z}(\mathbf{r},t)=\hat{H}'_{z}(\mathbf{r},t-\delta t)
-\frac{1}{\tilde{\mu}^{33}}\left(\frac{c_{0}\delta t}{Q_{0}}\right)^{2}
\left\{\left[\nabla_{q}\times\hat{\mathbf{E}}(\mathbf{r},t)\right]_{z}
+\omega_{z} \delta t \sum_{n=0}^{\infty}
\left[\nabla_{q}\times\hat{\mathbf{E}}(\mathbf{r},t-n\delta t)\right]_{z}
\right\}
\]
where we have made use of the fact that $1/(1-x)=\sum_{n=0}^{\infty}x^{n}$
to transform $1/ \omega^{-}$ into a sum over previous fields. These six
equations can now be used to update the electric and magnetic fields in the
absorbing region and will cause any waves within to be attenuated. The
strength of the attenuation is determined by the parameter $\omega_{z}$.

Note that while this analysis
has been for the purpose of setting up a PML we could use something similar
to treat an absorbing dielectric medium. In that case we would need 
to do something like,
\[
\varepsilon \mapsto \varepsilon_{r} \left( 1+\frac{i\gamma}{\omega^{-}}\right)
\]
and leave the $Q$'s just as the mesh spacings,
\[
Q_{1}=a \hspace{5mm};\hspace{5mm}
Q_{2}=b \hspace{5mm};\hspace{5mm}
Q_{3}=c \]
and then follow through the analysis as before.

\newpage
\section{The ONYX Program}

The ONYX program has been written and structured along a \emph{kit of parts}
philosophy. That is to say that since it is almost impossible to write a
single program which could perform every calculation which one might want
to do, the program has been constructed out of a number of largely
self-contained subroutines which can be bolted together in different ways
to perform different tasks. The subroutines loosely divide up into four 
groups: 
initialisation routines for setting the problem up, time integration routines
which make up the core of the calculation, post-processing routines for making
sense of the results and service routines for providing primitive matrix
operations, error trapping \textit{etc}. In the following sections we will
consider each subroutine in turn in the order in which they appear in the
program. For each routine we will give a table of all the variables used by
that routine, identifying which are used for the input / output of the routine
and which are purely internal variables. We will also try to give a description
of the routine; what it does and how it works.
In order to see how each subroutine fits into the overall program scheme, 
figure~\ref{fig:flow} gives an overview of the organisation of the whole
program.

\begin{figure}[hb]
\hrule\vspace{5mm}
\begin{center}
\include{flow_dia}
\end{center}
\hrule
\caption{Schematic to show the general structure of the ONYX program}
\label{fig:flow}
\end{figure}

\subsection{Preamble}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}
{Variables in modules \texttt{physconsts}, \texttt{files} 
and \texttt{parameters}}\\
\hline
\texttt{physconsts} & real & pi & Value for $\pi$ \\
 & real & c0 & Speed of light in free space (Atomic units)\\
 & real & eps0 & Free space permittivity (Atomic units)\\
 & real & mu0 & Free space permeability (Atomic units)\\
 & real & abohr & Bohr radius (metres)\\
 & real & x & Dummy parameter\\
 & real & emach & Machine accuracy\\
 & complex & ci & $\sqrt{-1}$\\
\hline
\texttt{files} & integer & logfile & Unit number for log.dat\\
 & integer & infile & Unit number for infile.dat\\
 & integer & outfile & Unit number for out.dat\\
 & integer & fields & Unit number for fields.dat\\
\hline
\texttt{parameters} 
 & integer & ixmax & Number of mesh points in the x-direction\\
 & integer & iymax & Number of mesh points in the y-direction\\
 & integer & izmax & Number of mesh points in the z-direction\\
 & integer & itmax & Total number of time steps=n\_block.block\_size\\
 & integer & ikmax & Number of k points\\
 & integer & n\_block & Number of blocks for time integration\\
 & integer & block\_size & Size of each block\\
 & integer & n\_pts\_store & Number of points where we store the fields\\
 & real & Q0 & Scale factor equal to typical mesh spacing\\
 & real & Q(3) & Array holding mesh spacing in each direction\\
 & real & dt & Size of the time step\\
 & real & akx & Rescaled x-component of $\mathbf{k}$. 
\texttt{akx=(Q(1)*ixmax)*kx}\\
 & real & aky & Rescaled y-component of $\mathbf{k}$. 
\texttt{aky=(Q(2)*iymax)*ky}\\
 & real & akz & Rescaled z-component of $\mathbf{k}$. 
\texttt{akz=(Q(3)*izmax)*kz}\\
 & integer & n\_segment & Number of segments for the FFT (leave as 1)\\
 & logical & overlap & True if segments are to overlap (Leave as false)\\
 & integer & fft\_size & If overlap=true then 
   \texttt{itmax=fft\_size*(n\_segment+1)}\\
 &  & & If overlap=false then \texttt{itmax=fft\_size*(2*n\_segment)}\\
 & real & damping & Controls amount fields are damped when stored\\
 & integer & bclayer1 & Thickness of absorbing layer on left-hand side\\
 & integer & bclayer2 & Thickness of absorbing layer on right-hand side\\
 & complex & epsref & Dielectric constant of host material\\
\hline
\end{tabular}
\end{center}
The Order N program begins with a block of comments. Most important are
the lines which identify the version number and the date of that version.
Following the comments are a series of \textbf{modules}. Most of these contain
\textbf{interface} blocks which serve two purposes. Firstly, they allow the
Fortran~90 complier to perform proper type-checking to help prevent errors
caused by passing the wrong type or number of arguments to a subroutine. 
Secondly, this explicit interface syntax is required by Fortran~90 if any
of the parameters passed to a subroutine are the new style Fortran~90
 \textbf{pointers}.
Apart from the interfaces there are three other important modules;
\texttt{files} which puts the unit numbers for the various input and output
files used by the program into sensibly named variables,
\texttt{parameters} which contains variables for all the constants read in
by the program at execution time and \texttt{physconsts} which defines
variables for the physical constants needed by the program.
Notice that the parameters \texttt{bclayer1}, \texttt{bclayer2} and
\texttt{epsref} are not defined in the input file but are set in the
subroutine \texttt{defcell}.

\subsection{\texttt{program onyx}}

\begin{center}

\begin{tabular}{||l|l|l||}
\hline
\multicolumn{3}{||c||}{Variables for program \texttt{onyx}}\\
\hline
integer & ios & Input / output status\\
integer & ik,ik2,ik3 & Loop counters for $\mathbf{k}$\\
integer & ix\_cur,iy\_cur,iz\_cur & x, y, z co-ordinates\\
integer & i\_pol & Field component label\\
integer,pointer & store\_pts(:,:) & Positions at which we store the fields\\
complex,pointer & e(:,:,:,:) & Electric Field at the current time\\
complex,pointer & h(:,:,:,:) & Magnetic Field at the current time\\
complex,pointer & eps(:,:,:) & Physical permittivity\\
complex,pointer & mu(:,:,:) & Physical permeability\\
complex,pointer & eps\_inv(:,:,:,:,:) & Inverse of renormalised permittivity\\
complex,pointer & mu\_inv(:,:,:,:,:) & Inverse of renormalised permeability\\
complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
complex,pointer & fft\_data(:,:) & Array for storing the fields\\
real,pointer & sigma(:,:,:) & Electric conductivity\\
real,pointer & sigma\_m(:,:,:) & Magnetic conductivity\\
real,pointer & spectrum(:) & Array to store the power spectrum or the DOS\\
real & u1(3),u2(3),u3(3) & Unit lattice vectors\\
real & g1(3),g2(3),g3(3) & Reciprocal lattice vectors\\
real & g(3,3) & Metric tensor\\
real & omega & $|\mathbf{u_1\cdot u_2\times u_3}|$\\
\hline
\end{tabular}
\end{center}
The main program itself does very little, other than setting the
calculation up and passing control to other subroutines which do the
hard work. 
The main program does contain some important loops however. There is
the possibility to loop over \texttt{akx}, \texttt{aky} and \texttt{akz}
which are the variables that control the values of $k_{x}$, $k_{y}$ and
$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,
a density of states calculation or a transmission 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\\
internal & complex,pointer & intcurle\_1(:,:,:) & 
Integral of $(\nabla\times\mathbf{E})_{z}$ on left side of system\\
internal & complex,pointer & intcurlh\_1(:,:,:) & 
Integral of $(\nabla\times\mathbf{H})_{z}$ on left side of system\\
internal & complex,pointer & intcurle\_2(:,:,:) & 
Integral of $(\nabla\times\mathbf{E})_{z}$ on right side of system\\
internal & complex,pointer & intcurlh\_2(:,:,:) & 
Integral of $(\nabla\times\mathbf{H})_{z}$ on right side of system\\
internal  & real,pointer & wz\_1(:) & 
Strength of absorber on left side of system\\
internal  & real,pointer & wz\_2(:) & 
Strength of absorber on right side of system\\
\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. 
It also sets up various arrays required by the absorbing boundary layers
including making a call to the subroutine \texttt{set\_wz} which defines the
strengths of the absorbing layers.
All 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(:,:,:,:,:)