onyx.tex

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

where,
\[
\hat{\sigma}^{ij}(\mathbf{r})=\sigma(\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{\sigma}_{m}^{ij}(\mathbf{r})=\sigma_m(\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,
\[
\Delta^{+}_{\tau}\hat{\mathbf{E}}(\mathbf{r},t)=
-\frac{\delta t \sigma}{\varepsilon_{0}\varepsilon(\mathbf{r})}
\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)=
-\frac{\delta t \sigma_{m}}{\mu_{0}\mu(\mathbf{r})}
\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)
\]
Writing 
\[
\hat{\sigma}=\frac{\delta t \sigma}{\varepsilon_{0}\varepsilon(\mathbf{r})}
\] 
and 
\[
\hat{\sigma}_{m}=\frac{\delta t \sigma_{m}}{\mu_{0}\mu(\mathbf{r})}
\]
Then,
\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*}
These are the update equations for the fields that we require.

\subsection{Stability Criterion}
These equations give a stable updating procedure for the fields if
the time step is kept sufficiently small. The criterion is easy to find.
Starting from the dispersion on the discrete mesh,
\[
\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\}
\]
The condition that the maximum value of the right hand side must
correspond to a real frequency gives,
\[
(\delta t)^2<\left(\frac{c_0^2}{Q_1^2}+\frac{c_0^2}{Q_2^2}
+\frac{c_0^2}{Q_3^2}\right)^{-1}
\]

\section{Program Overview}

\begin{center}
\begin{picture}(400,360)
%\put(0,0){\framebox(400,360){ }}
% Column 1
\put(40,200){\makebox(0,0){\fbox{Program ONYX}}}
\put(79,200){\line(1,0){11}}
\put(90,20){\line(0,1){300}}
\put(90,20){\vector(1,0){28}}
\put(90,50){\vector(1,0){24}}
\put(90,90){\vector(1,0){43}}
\put(90,180){\vector(1,0){28}}
\put(90,230){\vector(1,0){15}}
\put(90,260){\vector(1,0){40}}
\put(90,280){\vector(1,0){35}}
\put(90,300){\vector(1,0){42}}
\put(90,320){\vector(1,0){36}}
% Column 2
\put(150,20){\makebox(0,0){\fbox{finalproc\_dos}}}
\put(150,50){\makebox(0,0){\fbox{\shortstack{postproc\_band
\\postproc\_dos}}}}
\put(150,90){\makebox(0,0){\fbox{driver}}}
\put(150,180){\makebox(0,0){\fbox{\shortstack{initfields
\\initfields\_dos}}}}
\put(150,230){\makebox(0,0){\fbox{\shortstack{init\_store\_pts\_band
\\init\_store\_pts\_dos}}}}
\put(150,260){\makebox(0,0){\fbox{renorm}}}
\put(150,280){\makebox(0,0){\fbox{defmetric}}}
\put(150,300){\makebox(0,0){\fbox{defcell}}}
\put(150,320){\makebox(0,0){\fbox{setparam}}}
\put(174,280){\vector(1,0){69}}
\put(169,260){\vector(1,0){79}}
\put(182,180){\line(1,0){168}}
\put(240,180){\vector(0,-1){12}}
\put(350,180){\vector(0,-1){13}}
\put(167,90){\line(1,0){73}}
\put(240,90){\vector(0,1){42}}
\put(240,100){\vector(1,0){20}}
\put(240,120){\line(1,0){60}}
\put(300,120){\line(0,1){40}}
\put(300,160){\vector(1,0){41}}
\put(185,50){\vector(1,0){57}}
\put(210,50){\line(0,-1){20}}
\put(210,30){\vector(1,0){126}}
\put(315,100){\line(0,1){40}}
\put(315,140){\line(1,0){35}}
\put(350,140){\vector(0,1){12}}
% Column 3
\put(270,50){\makebox(0,0){\fbox{power\_spec}}}
\put(270,100){\makebox(0,0){\fbox{int}}}
\put(240,150){\makebox(0,0){\fbox{\shortstack{calc\_div\_B\\
calc\_div\_D\\calc\_energy\_density}}}}
\put(270,260){\makebox(0,0){\fbox{matinv3}}}
\put(270,280){\makebox(0,0){\fbox{cross\_prod}}}
\put(280,100){\line(1,0){35}}
\put(300,50){\line(1,0){50}}
\put(350,50){\vector(0,-1){8}}
% Column 4
\put(350,320){\makebox(0,0){\fbox{err}}}
\put(350,160){\makebox(0,0){\fbox{bc}}}
\put(350,30){\makebox(0,0){\fbox{\shortstack{FFT\\FT}}}}

\end{picture}
\end{center}

\section{The Main Program}
\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\\
\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.

\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