onyx2.tex

计算光子晶体能带的程序

\documentclass{article}
\usepackage{graphics}
\pagestyle{myheadings}
\markright{Order N program Documentation}
\oddsidemargin 0mm
\textwidth 16cm
\newcommand{\Exp}[1]{ {e^{#1}} }
%\textheight 24cm
\begin{document}

\title{Order N program Documentation}
\author{A. J. Ward \\ Imperial College London}

\begin{titlepage}
\maketitle
\thispagestyle{empty}
\begin{abstract}
The following document is intended 
to accompany the second version of the Order N method (Finite Difference 
Time Domain) electromagnetism program that we have been developing. 
The program is written in Fortran 90.

The code can perform three basic types
of calculation. The first is  a band structure calculation,
 which starts off with an
arbitrary set of initial fields corresponding to some wavevector $\mathbf{k}$,
integrates in time storing the fields at some random sampling points and then
Fourier transforms into the frequency domain. The peaks in the frequency 
spectrum correspond to the frequencies of the allowed modes.
The second type of calculation finds the time dependent Green's function
by setting the initial fields to be a delta function in space for one of the
field components, zero for all the others. The fields are integrated in time
and stored at each time step to give $g(t,\mathbf{r},\mathbf{r'})$. Fourier
transforming gives $G(\omega,\mathbf{r},\mathbf{r'})$.
The third type of calculation finds the transmission and reflection
coefficients of a finite system. This is achieved by projecting the scattered
wavefield onto a complete basis set of plane waves generated using the
transfer matrix on the discrete lattice. The current contribution from each
mode is then summed and divided by the current in the incident wave to give
the transmission and reflection coefficients.

Program code is presented in \texttt{typewriter} font. Fortran 90 keywords
 are in \textbf{bold}.
\end{abstract}
\end{titlepage}

\tableofcontents
\newpage

\section{Background Theory}
\subsection{Deriving the discrete Maxwell's equations}

We wish to derive discrete versions of Maxwell's equations which
can be used as the basis for a finite difference time domain calculation
of a form similar to that of Yee~\cite{Yee}.
Starting from Maxwell's equations,
\[
\nabla\times\mathbf{H}=\varepsilon_{0}\varepsilon_{r}
\frac{\partial\mathbf{E}}{\partial t}
\hspace{1cm}
\nabla\times\mathbf{E}=-\mu_{0}\mu_{r}
\frac{\partial\mathbf{H}}{\partial t}
\]

Next we want to use the results for writing Maxwell's equations in
generalised co-ordinates that were developed in~\cite{AJW2} and
explained in detail in chapter four of my thesis~\cite[pages
66-93]{mythesis}.
The basic
results from that work are very simple and mean that in the generalised
co-ordinates we can write,
\[
\nabla_{q}\times\hat{\mathbf{H}}=\varepsilon_{0}\hat{\varepsilon}(\mathbf{r})
\frac{\partial\hat{\mathbf{E}}}{\partial t}
\hspace{1cm}
\nabla_{q}\times\hat{\mathbf{E}}=-\mu_{0}\hat{\mu}(\mathbf{r})
\frac{\partial\hat{\mathbf{H}}}{\partial t}
\]
The renormalised fields are,
\[
\hat{E}_{i}=Q_{i}E_{i}\hspace{2cm}
\hat{H}_{i}=Q_{i}H_{i}
\]
with,
\[
Q_{i}=\sqrt{
 \left(\frac{\partial x}{\partial q_{i}}\right)^{2}
+\left(\frac{\partial y}{\partial q_{i}}\right)^{2}
+\left(\frac{\partial z}{\partial q_{i}}\right)^{2}}
\]
The effective permittivity $\hat{\varepsilon}$ and permeability $\hat{\mu}$
can be shown to be tensors of the following form,
\[
\hat{\varepsilon}^{ij}(\mathbf{r})=\varepsilon(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}}
\]
\[
\hat{\mu}^{ij}(\mathbf{r})=\mu(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}}
\]
Where the vectors $\mathbf{u}_i$ are unit vectors pointing along the
axes of the generalised co-ordinate system.
The point of this result is that all the details of the co-ordinate
system have been `hidden' inside the effective $\varepsilon$ and $\mu$.
All we need do is ensure that we use $\hat{\varepsilon}^{ij}$ and
 $\hat{\mu}^{ij}$ at each point.

Now we want to write down the discrete versions of these time dependent
equations. To this end we replace the temporal and spatial derivatives
with finite differences. We choose to use forward differences for 
derivatives of the E-field and backwards differences for the H-fields. 
This choice leads to conditionally stable equations.
So\ldots
\[
\frac{\partial\hat{E}}{\partial q_{i}}\mapsto 
\Delta^{+}_{i}\hat{E}=
\hat{E}(\mathbf{r+Q_{i}},t)-\hat{E}(\mathbf{r},t)
\]
\[
\frac{\partial\hat{H}}{\partial q_{i}}\mapsto 
\Delta^{-}_{i}\hat{H}=
\hat{H}(\mathbf{r},t)-\hat{H}(\mathbf{r-Q_{i}},t)
\]
and
\[
\frac{\partial\hat{E}}{\partial t}\mapsto 
\frac{\Delta^{+}_{\tau}\hat{E}}{\delta t}=
\frac{\hat{E}(\mathbf{r},t+\delta t)-\hat{E}(\mathbf{r},t)}{\delta t}
\]
\[
\frac{\partial\hat{H}}{\partial t}\mapsto 
\frac{\Delta^{-}_{\tau}\hat{H}}{\delta t}=
\frac{\hat{H}(\mathbf{r},t)-\hat{H}(\mathbf{r},t-\delta t)}{\delta t}
\]
Where $\mathbf{Q}_i=Q_i\mathbf{u}_i$.
We can use the discrete spatial derivatives to define a discrete version
of the curl,
\[
\left[\nabla_{q}^{+}\times\mathcal{F}\right]_{3}=
\Delta_{1}^{+}\mathcal{F}_{2}-\Delta_{2}^{+}\mathcal{F}_{1}
\]
etc. Substituting into Maxwell's equations gives,
\[
\nabla_{q}^{+}\times\hat{\mathbf{E}}=-\frac{\mu_{0}}{\delta t}
\hat{\mu}(\mathbf{r})\Delta^{+}_{\tau}\hat{\mathbf{H}}
\]
\[
\nabla_{q}^{-}\times\hat{\mathbf{H}}=+\frac{\varepsilon_{0}}{\delta t}
\hat{\varepsilon}(\mathbf{r})\Delta^{-}_{\tau}\hat{\mathbf{E}}
\]
Here we should note that on the discrete mesh the $Q_i$ factors simply
correspond to the spacings between the mesh points in each direction.

This choice of finite differences is equivalent to making the following
approximations to $\mathbf{k}$ and $\omega$. For the electric field
terms,
\[
k_j\mapsto k^{+}_j=\frac{\exp{[ik_j Q_j]}-1}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^{+}=\frac{\exp{[-i\omega\delta t]}-1}{-i\delta t}
\]
And for the magnetic field terms,
\[
k_j\mapsto k^{-}_j=\frac{1-\exp{[-ik_j Q_j]}}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^{-}=\frac{1-\exp{[i\omega\delta t]}}{-i\delta t}
\]
For the free-space dispersion relation, $\omega/k=c_0$ these 
approximations give,
\[
\frac{4}{\delta t^2}\sin^2\left(\frac{\omega \delta t}{2}\right)=
4c_0^2\left\{
 \frac{1}{Q_1^2}\sin^2\left(\frac{Q_1\ k_x}{2}\right)
+\frac{1}{Q_2^2}\sin^2\left(\frac{Q_2\ k_y}{2}\right)
+\frac{1}{Q_3^2}\sin^2\left(\frac{Q_3\ k_z}{2}\right)\right\}
\]

For reasons of numerical stability it is desirable to make the electric
and magnetic fields of roughly equal magnitude. This can be done by 
introducing a rescaled magnetic field.
\[
\hat{\mathbf{H}}'=\frac{\delta t}{\varepsilon_{0}Q_{0}}\hat{\mathbf{H}}
\]
Here we introduce the constant $Q_0$ to have the dimensions of length
and be roughly equal in length to the $Q_i$'s.
If we also introduce a rescaling to $\hat{\varepsilon}$ and $\hat{\mu}$
then Maxwell's equations become,
\[
\Delta^{+}_{\tau}\hat{\mathbf{E}}(\mathbf{r},t)=
\left[\hat{\varepsilon}(\mathbf{r})\right]^{-1}
\nabla^{-}_{q}\times \hat{\mathbf{H}}'(\mathbf{r},t)
\]
\[
\Delta^{-}_{\tau}\hat{\mathbf{H}}'(\mathbf{r},t)=
-\left\{\frac{\delta t\ c_{0}}{Q_{0}}\right\}^{2}
\left[\hat{\mu}(\mathbf{r})\right]^{-1}
\nabla^{+}_{q}\times \hat{\mathbf{E}}(\mathbf{r},t)
\]
where $\hat{\varepsilon}$ and $\hat{\mu}$ are now,
\[
\hat{\varepsilon}^{ij}(\mathbf{r})=\varepsilon(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}Q_{0}}
\]
\[
\hat{\mu}^{ij}(\mathbf{r})=\mu(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}Q_{0}}
\]
So if we know the electric and magnetic fields at time $t$ we can
calculate the fields at $t+\delta t$
\begin{eqnarray*}
\hat{E}_{1}(\mathbf{r},t+\delta t)=\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)=\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)=\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)&=&\hat{H}'_{1}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\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\}\right.\\
&+&\left.
\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\}\right.\\
&+&\left.
\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)&=&\hat{H}'_{2}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\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\}\right.\\
&+&\left.
\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\}\right.\\
&+&\left.
\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)&=&\hat{H}'_{3}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\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\}\right.\\
&+&\left.
\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\}\right.\\
&+&\left.
\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*}
where we have written $\mathbf{a}\equiv\mathbf{Q}_1$, 
$\mathbf{b}\equiv\mathbf{Q}_2$ and $\mathbf{c}\equiv\mathbf{Q}_3$.

\subsection{Introducing the current terms}

All that remains now is to include the current terms.
This is done by means of a finite electrical conductivity. For historical
reasons a ficticious magnetic
current is also introduced by including a magnetic conductivity.