onyx.tex

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

&+&
\left.
\left[
\varepsilon(\mathbf{r})g^{1i}\hat{E}_{i}(\mathbf{r})
-\varepsilon(\mathbf{r-a})g^{1i}\hat{E}_{i}(\mathbf{r-a})\right]
\frac{Q_{1}Q_{2}Q_{3}}{Q_{0}Q_{i}Q_{1}} \right.\\
&+&
\left.
\left[
\varepsilon(\mathbf{r})g^{2i}\hat{E}_{i}(\mathbf{r})
-\varepsilon(\mathbf{r-b})g^{2i}\hat{E}_{i}(\mathbf{r-b})\right]
\frac{Q_{1}Q_{2}Q_{3}}{Q_{0}Q_{i}Q_{2}}
\right\} \\
Q/(Q_{0}\varepsilon_{0})&=&
 \hat{\varepsilon}^{3i}(\mathbf{r})\hat{E}_{i}(\mathbf{r})
-\hat{\varepsilon}^{3i}(\mathbf{r-c})\hat{E}_{i}(\mathbf{r-c}) \\
&+&\hat{\varepsilon}^{1i}(\mathbf{r})\hat{E}_{i}(\mathbf{r})
-\hat{\varepsilon}^{1i}(\mathbf{r-a})\hat{E}_{i}(\mathbf{r-a}) \\
&+&\hat{\varepsilon}^{2i}(\mathbf{r})\hat{E}_{i}(\mathbf{r})
-\hat{\varepsilon}^{2i}(\mathbf{r-b})\hat{E}_{i}(\mathbf{r-b}) \\
&=&\nabla^{-}\cdot\hat{\varepsilon}(\mathbf{r})
\mathbf{\hat{E}(r)}
\end{eqnarray*}
So instead of returning $Q$ the subroutine \texttt{calc\_div\_D} actually
returns $Q/(Q_{0}\varepsilon_{0})$ but that only differs from $Q$ by a
constant.
The same follows for calculating $\nabla\cdot\mathbf{B}$

\subsection{\texttt{subroutine calc\_energy\_density}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine 
\texttt{calc\_energy\_density}}\\
\hline
input & complex,pointer & e\_cur(:,:,:,:) 
& Electric Field at the current time\\
input & complex,pointer & e\_prev(:,:,:,:) 
& Electric Field at the previous time\\
input & complex,pointer & h\_prev(:,:,:,:) 
& Magnetic Field at the previous time\\
output & real,pointer & rho(:,:,:) & Energy density\\
input & complex,pointer & eps\_hat(:,:,:,:,:) & Renormalised permittivity\\
input & complex,pointer & mu\_hat(:,:,:,:,:) & Renormalised permeability\\
\hline
internal & integer & ix,iy,iz,i,j & Loop counters \\
internal & real & dtcq2 & $(\delta t\ c_0/Q_0)^2$\\
\hline
\end{tabular}
\end{center}
The correct discrete form for the energy density is\ldots
\[
U=\frac{1}{2}\left\{
\varepsilon_{0}\varepsilon(\mathbf{r})E_{\alpha}^{*}(t)
E^{\alpha}(t-\delta t)
+\mu_{0}\mu(\mathbf{r})H_{\alpha}^{*}(t-\delta t)H^{\alpha}(t-\delta t)
\right\}
\]
This has the correct form for the energy density in the limit $\delta t
\rightarrow 0$ and can be shown to be an exactly conserved quantity on
the lattice. Substituting the usual expressions for the generalised
co-ordinates you get\ldots
\[
U=\frac{U_{0}}{2}\left\{
\hat{\varepsilon}^{\alpha\beta}\hat{E}^{*}_{\alpha}(t)
\hat{E}_{\beta}(t-\delta t)
+(\frac{Q_{0}}{c_{0}\delta t})^{2}
\hat{\mu}^{\alpha\beta}\hat{H'}^{*}_{\alpha}(t-\delta t)
\hat{H'}_{\beta}(t-\delta t)
\right\}
\]
where $ U_{0}=Q_{0}\varepsilon_{0}/(Q_{1}Q_{2}Q_{3}|\mathbf{u_{1}\cdot
u_{2}\times u_{3}}|)$.

\subsection{\texttt{subroutine calc\_current}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{calc\_current}}\\
\hline
output & complex & J(3) & Energy current away from mesh point\\
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 & h\_prev(:,:,:,:) 
& Magnetic Field at the previous time\\
input & integer & ix,iy,iz & Position co-ordinates\\
\hline
\end{tabular}
\end{center}
The Poynting vector can be obtained fairly simply by considering
$\Delta_{t} U(t)=(U(t+\delta t)-U(t))/\delta t$. This gives\ldots
\begin{eqnarray*}
\Delta^{+}_{t}U(t)&=&\frac{U_{0}}{2}
\left\{\left[\hat{H}'^{*}_{2}(\mathbf{r-c},t)
-\hat{H}'^{*}_{3}(\mathbf{r-b},t)\right]\hat{E}_{1}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{H}'^{*}_{3}(\mathbf{r-a},t)
-\hat{H}'^{*}_{1}(\mathbf{r-c},t)\right]\hat{E}_{2}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{H}'^{*}_{1}(\mathbf{r-b},t)
-\hat{H}'^{*}_{2}(\mathbf{r-a},t)\right]\hat{E}_{3}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{H}'_{2}(\mathbf{r-c},t-\delta t)
-\hat{H}'_{3}(\mathbf{r-b},t-\delta t)\right]
\hat{E}^{*}_{1}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{H}'_{3}(\mathbf{r-a},t-\delta t)
-\hat{H}'_{1}(\mathbf{r-c},t-\delta t)\right]
\hat{E}^{*}_{2}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{H}'_{1}(\mathbf{r-b},t-\delta t)
-\hat{H}'_{2}(\mathbf{r-a},t-\delta t)\right]
\hat{E}^{*}_{3}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}_{2}(\mathbf{r+c},t)-\hat{E}_{3}(\mathbf{r+b},t)\right]
\hat{H}'^{*}_{1}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}_{3}(\mathbf{r+a},t)-\hat{E}_{1}(\mathbf{r+c},t)\right]
\hat{H}'^{*}_{2}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}_{1}(\mathbf{r+b},t)-\hat{E}_{2}(\mathbf{r+a},t)\right]
\hat{H}'^{*}_{3}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}^{*}_{2}(\mathbf{r+c},t)
-\hat{E}^{*}_{3}(\mathbf{r+b},t)\right]
\hat{H}'_{1}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}^{*}_{3}(\mathbf{r+a},t)
-\hat{E}^{*}_{1}(\mathbf{r+c},t)\right]
\hat{H}'_{2}(\mathbf{r},t) \right.\\
&+&\left.\left[\hat{E}^{*}_{1}(\mathbf{r+b},t)
-\hat{E}^{*}_{2}(\mathbf{r+a},t)\right]
\hat{H}'_{3}(\mathbf{r},t)
\right\}
\end{eqnarray*}
These terms can be divided up into those which carry energy \emph{into} 
a given mesh point and those which carry energy \emph{away} from it. 
The subroutine
\texttt{calc\_current} returns the three components of the energy flow 
\emph{away} from the given mesh point.

\section{Miscellaneous subroutines}

\subsection{\texttt{subroutine err}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{err}}\\
\hline
input & integer & ierr & Error number \\
\hline
\end{tabular}
\end{center}
This is a very simple subroutine which takes as its argument the
integer \texttt{ierr} and depending on its value writes a relevant
error message to the log file and then stops execution of the 
program. This allows other subroutines to halt the calculation when
they encounter a problem and leave the user with some clue as to
what has gone wrong. 

\subsection{\texttt{function matinv3}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for function \texttt{matinv3}}\\
\hline
output & complex & matinv3(3,3) & $A^{-1}$ \\
input & complex & A(3,3) & $3\times 3$ matrix\\
input & real & emach & Machine accuracy\\
output & logical & fail & True if $A$ has no inverse\\
\hline
internal & complex & det & $|A|$\\
internal & complex & B(3,3) & Work space\\
\hline
\end{tabular}
\end{center}
The function \texttt{matinv3} returns the inverse of the 
three by three matrix A. If A has no inverse then \texttt{fail} is
set to true else  \texttt{fail} is false.  \texttt{emach} is the machine
accuracy.

\subsection{\texttt{function cross\_prod}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for function \texttt{cross\_prod}}\\
\hline
output & real & cross\_prod(3) & $\mathbf{a\times b}$ \\
input & real & a(3),b(3) & Two 3-vectors \\
\hline
\end{tabular}
\end{center}
The function \texttt{cross\_prod} returns the cross-product of the two
three-vectors that the function takes as arguments. 

\subsection{\texttt{subroutine power\_spec} and \texttt{subroutine FFT}} 

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{power\_spec}}\\
\hline
input & complex,pointer & f(:,:) & Array of functions to be transformed\\
output & real,pointer & spectrum(:) & Power Spectrum\\
\hline
internal & integer & j,k & Loop counters\\
internal & integer & mm,m2,offset & Internal variables\\
internal & real & den,facm,facp,sumw & Internal variables\\
internal & real & window & Window function\\
internal & complex,pointer & w1(:,:),w2(:,:) & Work space\\
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{FFT}}\\
\hline
in/out & complex,pointer & f(:,:) & Array of functions to be transformed\\
input & integer & isign & +1 for transform, -1 for inverse transform\\
\hline
internal & complex & tmp,w,wp & Internal variables\\
internal & real & theta & Frequency\\
internal & integer & n & Number of functions to be transformed\\
internal & integer & nmax & Length of each function\\
internal & integer & nbits & Bit size of an integer\\
internal & integer & set\_bits & Counter for number of set bits\\
internal & integer & i,j,k,m,mmax,istep & Internal variables\\
\hline
\end{tabular}
\end{center}
The subroutines \texttt{FFT} and \texttt{power\_spec} are standard routines
for performing a fast Fourier transform and using that to obtain a power
spectrum. The code follows directly from Numerical Recipes pg 425ff.
The only substantial alteration is the inclusion of the loop after the 
fast Fourier transform has been completed,
\begin{verbatim}
do i=1,nmax/2
  w=2.0*pi*(i-1)/real(nmax)
  f(1,i)=f(1,i)*cexp(+isign*ci*w*0.50)
enddo
\end{verbatim}
This multiplies each frequency component by a factor $\exp{[i\omega/2]}$
which ensures that the Fourier transform we have calculated corresponds to,
\[
F(\omega)=\sum_{n=1}^{N_{t}} f(n\,\delta t) \exp{[i\omega(n-0.5)\delta t]}
\]
rather than, 
\[
F(\omega)=\sum_{n=1}^{N_{t}} f(n\,\delta t) \exp{[i\omega(n-1)\delta t]}
\]
which is what the standard Numerical Recipes fast Fourier transform
would give.
The half time step offset was found to be necessary in order to correctly
calculate the Green's function and avoid negatives in the density of states.

\subsection{\texttt{subroutine FT}}

\begin{center}
\begin{tabular}{||l|l|l|l||}
\hline
\multicolumn{4}{||c||}{Variables for subroutine \texttt{FT}}\\
\hline
in/out & complex,pointer & f(:,:) & Array of functions to be transformed\\
input & integer & isign & +1 for transform, -1 for inverse transform\\
\hline
internal & complex,pointer & ft\_f(:) 
& Temporary array to hold the transforms\\
internal & complex & tmp & Temporary variable\\
internal & real & w & Frequency\\
internal & integer & n & Number of functions to be transformed\\
internal & integer & nmax & Length of each function\\
internal & integer & iw,i,it & Loop counters\\
\hline
\end{tabular}
\end{center}
The subroutine \texttt{FT} performs a straight-forward slow Fourier transform
except for the half a time step offset in the exponential part. So what we
are calculating is,
\[
F(\omega)=\sum_{n=1}^{N_{t}} f(n\,\delta t) \exp{[i\omega(n-0.5)\delta t]}
\]

The array \texttt{f(n,nmax)} contains \texttt{n} separate functions of length
\texttt{nmax} each of which is Fourier transformed separately. The parameter
\texttt{isign} determines whether a inverse transform is performed.
If \texttt{isign}=1 it calculates the transform, 
if \texttt{isign}=-1 it calculates inverse transform, except for the
normalisation factor.

\section{Input / Output}
\label{sec:input}

\begin{center}
\begin{tabular}{||l|l|l||}
\hline
\multicolumn{3}{||c||}{Parameters read in from \texttt{infile.dat}}\\
\hline
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\\
real & dx & Mesh spacing in x direction in atomic units\\
real & dy & Mesh spacing in y direction in atomic units\\
real & dz & Mesh spacing in z direction in atomic units\\
real & dt & Size of the time step in atomic units\\
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 & 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\\
integer & n\_segment & Number of segments for the FFT (leave as 1)\\
logical & overlap & True if segments are to overlap (Leave as false)\\
real & damping & Controls amount fields are damped when stored\\
\hline
\e