onyx.tex

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

if required.

A special version of this subroutine \texttt{initfields\_dos} is also
provided, for use when calculating the density of states, which only has code
to set the starting fields to be a delta function.

\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} 
and \texttt{subroutine init\_store\_pts\_dos}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutines 
\texttt{init\_store\_pts\_band} and \texttt{init\_store\_pts\_dos}}\\
\hline
output & integer,pointer & store\_pts(:,:) 
& Positions at which we store the fields\\
input & integer & ix\_cur,iy\_cur,iz\_cur & x, y, z co-ordinates\\
input & integer & i\_pol & Field component label\\
\hline
internal & integer & i & Loop counter\\
internal & real & x & Random number\\
\hline
\end{tabular}
\end{center}
These subroutines  initialises the array 
\texttt{store\_pts} to hold information about which field components
are to be stored during the time integration into the array 
\texttt{fft\_data} for later post-processing. 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 parameter \texttt{n\_pts\_store} controls how many
field components are to be stored.
For a band structure calculation \texttt{subroutine init\_store\_pts\_band}
fills the array \texttt{store\_pts} randomly
while for a Green's function calculation 
\texttt{subroutine init\_store\_pts\_dos} sets up \texttt{store\_pts} such that
only the \texttt{i\_pol}'th component
of the fields at the position \texttt{ix\_cur}, \texttt{iy\_cur}, 
\texttt{iz\_cur} is stored.

\section{Subroutines for  post-processing}

\subsection{\texttt{subroutine postproc\_band} and
 \texttt{subroutine postproc\_dos}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutines \texttt{postproc\_band}
and \texttt{postproc\_dos}}\\
\hline
input & complex,pointer & fft\_data(:,:) & Array for storing the fields\\
output & real,pointer & spectrum(:) & Array for the power spectrum\\
& & &                                 or the density of states\\
\hline
internal & integer & i & Loop counter\\
\hline
\end{tabular}
\end{center}
The \texttt{postproc\_band} subroutine provides any necessary post-processing
on the time-dependent fields that have been stored in the array
\texttt{fft\_data}. For the band structure calculation this involves a
Fourier transform by means of the subroutine \texttt{power\_spec} after which
the band structure is written out analysing the derivative of the spectrum
and identifying the peaks. If \texttt{ikmax} is set to one the routine
does not write out the band structure but instead writes the power spectrum
as returned by \texttt{power\_spec}. This was originally done to serve as a
test to make sure the spectrum makes physical sense.

For a calculation of the density of
states the subroutine \texttt{postproc\_dos} provides post-processing by
calling the subroutine \texttt{FFT} to perform a simple Fourier transform
and then adding minus
the imaginary part of the result to the cumulative density
of states stored in the array \texttt{spectrum}.

\subsection{\texttt{subroutine finalproc\_dos}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{finalproc\_dos}}\\
\hline
input & real,pointer & spectrum(:) & Array for the density of states\\
\hline
internal & integer & i & Loop counter\\
\hline
\end{tabular}
\end{center}
The subroutine \texttt{finalproc\_dos} performs the final post-process step
to write out the complete density of states to the output file.

\section{Subroutines to handle generalised co-ordinates}

\subsection{\texttt{subroutine defmetric}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{defmetric}}\\
\hline
input & real & u1(3),u2(3),u3(3) & Unit lattice vectors\\
output & real & g1(3),g2(3),g3(3) & Reciprocal lattice vectors\\
output & real & g(3,3) & Metric tensor\\
output & real & omega & $|\mathbf{u_1\cdot u_2\times u_3}|$\\
\hline
internal & real & tmp(3) & Temporary vector\\
\hline
\end{tabular}
\end{center}
The subroutine \texttt{defmetric} takes the three unit vectors \texttt{u1},
\texttt{u2} and \texttt{u3} and defines some of the quantities required when
setting up the generalised co-ordinate system. The three vectors \texttt{g1},
\texttt{g2} and \texttt{g3} are defined to be the set of vectors reciprocal to
the \texttt{u}'s so that $\mathbf{g_{1}}=(\mathbf{u_{2}\times u_{3}})/
(|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|)$. Therefore 
$\mathbf{u_{\alpha}\cdot g_{\beta}=\delta_{\alpha\beta}}$. The routine
also defines the metric tensor for the co-ordinate system,
$G_{\alpha\beta}=\mathbf{g_{\alpha}\cdot g_{\beta}}$. Finally the
subroutine defines an extra constant \texttt{Q0} which has the dimensions
of a length and a magnitude roughly the same as the spacing between mesh 
points \texttt{Q(1)},\texttt{Q(2)},\texttt{Q(3)}. 
In fact the routine simply sets
\texttt{Q0=Q(1)} 

\subsection{\texttt{subroutine renorm}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{renorm}}\\
\hline
input & complex,pointer & eps(:,:,:) & Physical permittivity\\
input & complex,pointer & mu(:,:,:) & Physical permeability\\
in/out & real,pointer & sigma(:,:,:) & Electric conductivity\\
in/out & real,pointer & sigma\_m(:,:,:) & Magnetic conductivity\\
output & complex,pointer & eps\_inv(:,:,:,:,:) 
& Inverse of renormalised permittivity\\
output & complex,pointer & mu\_inv(:,:,:,:,:) 
& Inverse of renormalised permeability\\
output & complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
output & complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
input & real & g(3,3) & Metric tensor\\
input & real & omega & $|\mathbf{u_1\cdot u_2\times u_3}|$\\
\hline
internal & complex, dimension(3,3) & eps\_tmp,mu\_tmp & Temporary arrays\\
internal & integer & i,j,ix,iy,iz & Loop counters \\
internal & logical & fail & Flag for matrix inversion failure\\
\hline
\end{tabular}
\end{center}
The subroutine \texttt{renorm} takes the real, physical quantities 
$\varepsilon$,
 $\mu$, $\sigma$ and $\sigma_m$ and returns the renormalised 
 $\hat{\varepsilon}$, $\hat{\mu}$, $\hat{\sigma}$ and $\hat{\sigma}_{m}$
according to the formulae we gave previously,
\[
\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}}
\hspace{1cm}
\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}}
\]
 \[
\hat{\sigma}=\frac{\delta t\, \sigma}{\varepsilon_{0}\varepsilon(\mathbf{r})}
\hspace{5cm}
\hat{\sigma}_{m}=\frac{\delta t\, \sigma_{m}}{\mu_{0}\mu(\mathbf{r})}
\]
The
routine then calculates $\hat{\varepsilon}^{-1}$ and $\hat{\mu}^{-1}$ as
it is these quantities which are needed in the equations which update the
fields.

\section{Boundary conditions}

\subsection{Simple boundary conditions}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutines \texttt{bc\_*\_bloch} and
\texttt{bc\_*\_metal}}\\
\hline
in/out & complex,pointer & e(:,:,:,:) & Electric Field at the current time\\
in/out & complex,pointer & h(:,:,:,:) & Magnetic Field at the current time\\
\hline
internal & integer & ix,iy,iz,i & Loop counters \\
internal & complex & fxm1 & Bloch phase - $\exp{[-i\mathbf{k\cdot a}]}$\\
internal & complex & fxp1 & Bloch phase - $\exp{[i\mathbf{k\cdot a}]}$\\
internal & complex & fym1 & Bloch phase - $\exp{[-i\mathbf{k\cdot b}]}$\\
internal & complex & fyp1 & Bloch phase - $\exp{[i\mathbf{k\cdot b}]}$\\
internal & complex & fzm1 & Bloch phase - $\exp{[-i\mathbf{k\cdot c}]}$\\
internal & complex & fzp1 & Bloch phase - $\exp{[i\mathbf{k\cdot c}]}$\\
\hline
\end{tabular}
\end{center}
Two sets of simple boundary condition routines a provided. The first set,
\texttt{bc\_*\_bloch}, use Bloch's law, $F(r+R)=\exp{(ik\cdot R)}F(r)$ to
correctly set the fields along the planes 
\texttt{ix=0}, \texttt{iy=0}, \texttt{iz=0}, \texttt{ix=ixmax+1}, 
\texttt{iy=iymax+1} and \texttt{iz=izmax+1} respectively.
The values for the normalised $k$-vector components \texttt{akx}
,\texttt{aky},\texttt{akz}
are made available by using the module \texttt{parameters}.
The second set of subroutines, \texttt{bc\_*\_metal}, set the fields to
be zero along those planes.

\section{Subroutines to calculate physical quantities}

\subsection{\texttt{subroutine calc\_div\_D} and 
\texttt{subroutine calc\_div\_B}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutines \texttt{calc\_div\_D} and
\texttt{calc\_div\_B}}\\
\hline
input & complex,pointer & e(:,:,:,:) & Electric Field at the current time\\
input & complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
output & complex,pointer & div\_D(:,:,:) 
& Array to store $\nabla\cdot\mathbf{D}$\\
\hline
input & complex,pointer & h(:,:,:,:) & Magnetic Field at the current time\\
input & complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
output & complex,pointer & div\_B(:,:,:) 
& Array to store $\nabla\cdot\mathbf{B}$\\
\hline
internal & integer & ix,iy,iz,i & Loop counters \\
internal & integer & ixm1 & \texttt{ix-1}\\
internal & integer & iym1 & \texttt{iy-1}\\
internal & integer & izm1 & \texttt{iz-1}\\
internal & integer & ixp1 & \texttt{ix+1}\\
internal & integer & iyp1 & \texttt{iy+1}\\
internal & integer & izp1 & \texttt{iz+1}\\
\hline
\end{tabular}
\end{center}
%Starting from Maxwell equations,
%\[
%\nabla\times\mathbf{H}=\varepsilon_{0}
%\varepsilon(\mathbf{r})\frac{\partial \mathbf{E(r)}}{\partial t}
%\hspace{2cm}
%\nabla\cdot\varepsilon_{0}\varepsilon(\mathbf{r})\mathbf{E(r)}=\rho
%\]
%\[
%\nabla\times\mathbf{E}=-\mu_{0}\mu(\mathbf{r})\frac{\partial \mathbf{H(r)}}
%{\partial t}
%\hspace{2cm}
%\nabla\cdot\mu_{0}\mu(\mathbf{r})\mathbf{H(r)}=0
%\]
In fact our program actually
calculates $\nabla^{-}\cdot\hat{\varepsilon}(\mathbf{r})
\mathbf{\hat{E}(r)}$ and $\nabla^{+}\cdot\hat{\mu}(\mathbf{r})
\mathbf{\hat{H}(r)}$ but we can easily work out what those quantities
correspond to. Start from,
\[
\nabla\cdot\mathbf{D}=\rho
\]
So in integral form,
\[
\int_{\mathit{Vol}}\nabla\cdot\mathbf{D}\ dV=Q
\]
where $Q$ is the charge in volume $V$. Using the divergence theorem,
\[
\int_{\mathit{Surface}} \mathbf{D}\cdot dS=Q \\
\]
So integrating over the surface of the cell surrounding one of our
discrete mesh points\ldots
\begin{eqnarray*}
Q/\varepsilon_{0}&=&\varepsilon(\mathbf{r}) E^{3}(\mathbf{r})
\mathbf{u}_{3}\cdot 
Q_{1}Q_{2}(\mathbf{u_{1}\times u_{2}})
-\varepsilon(\mathbf{r-c}) E^{3}(\mathbf{r-c})\mathbf{u}_{3}\cdot Q_{1}Q_{2}
(\mathbf{u_{1}\times u_{2}}) \\
&+&\varepsilon(\mathbf{r}) E^{1}(\mathbf{r})\mathbf{u}_{1}\cdot Q_{2}Q_{3}
(\mathbf{u_{2}\times u_{3}})
-\varepsilon(\mathbf{r-a}) E^{1}(\mathbf{r-a})\mathbf{u}_{1}\cdot Q_{2}Q_{3}
(\mathbf{u_{2}\times u_{3}}) \\
&+&\varepsilon(\mathbf{r}) E^{2}(\mathbf{r})\mathbf{u}_{2}\cdot Q_{3}Q_{1}
(\mathbf{u_{3}\times u_{1}})
-\varepsilon(\mathbf{r-b}) E^{2}(\mathbf{r-b})\mathbf{u}_{2}\cdot Q_{3}Q_{1}
(\mathbf{u_{3}\times u_{1}}) \\
&=&Q_{0}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|
\left\{
\left[
\varepsilon(\mathbf{r})g^{3i}\hat{E}_{i}(\mathbf{r})
-\varepsilon(\mathbf{r-c})g^{3i}\hat{E}_{i}(\mathbf{r-c})\right]
\frac{Q_{1}Q_{2}Q_{3}}{Q_{0}Q_{i}Q_{3}} \right.\\