epsilon_polariton_1d.tex

麻省理工的计算光子晶体的程序

\begin{comment}
#include <stdio.h>
#include <stdlib.h>
#include <signal.h>

#include <meep.hpp>
using namespace meep;
\end{comment}

\section{Epsilon of a polaritonic material in one dimension}

In this example, we compute epsilon as a function of frequency for a simple
polaritonic material.  This example is done in one dimension for speed
purposes.

One thing to be aware of when using polaritonic materials, is that
generally you will be needing a rather higher grid resolution than you may
be used to in order to properly model the material.  Here I am using an $a$
of 40.

\begin{figure}
\label{epsilon_polariton}
\caption{Epsilon of a polaritonic material.}
\includegraphics[width=8.8cm,clip=true]{epsilon_polariton_1d-out/eps}
\end{figure}


\begin{comment}
const double a = 10;
const double pml_thickness = 1.0;
const double middlesize = 1.0;
const double zsize = middlesize + 2*pml_thickness;
\end{verbatim}

For our example polaritonic material, we'll use an $\epsilon(0)$ of 13.4. %'

\begin{verbatim}
double eps(const vec &) { return 13.4; }
\end{verbatim}
\begin{comment}
double one(const vec &p) {
  return 1;
}

int main(int argc, char **argv) {
  initialize mpi(argc, argv);
  deal_with_ctrl_c();
  const double ttot =2000.0;
  const volume v = volone(zsize, a);
  const symmetry S = mirror(Z, v);
  structure s(v, eps, pml(pml_thickness), S);
  const char *dirname = make_output_directory(__FILE__);
  s.set_output_directory(dirname);
\end{comment}
Altough in this calculation the polaritonic material will not be within the
PML, it is all right to have polaritonic material within PML regions.
\begin{verbatim}
  s.add_polarizability(one, 0.4, 0.01, 27.63);
\end{verbatim}
\begin{comment}
  fields f(&s);
  double sourceloc = pml_thickness+1.0/(double)a;
\end{comment}
We use a single rather high frequency (and very broad) point source, to
cover a broad frequency range.
\begin{verbatim}
  f.add_point_source(Ex, 0.9, 0.8, 0.0, 8.0, vec(sourceloc));
\end{verbatim}
We use a couple of monitor points to determine epsilon.
\begin{verbatim}
  monitor_point *left = NULL, *right = NULL, *middle = NULL;
\end{verbatim}
\begin{comment}
  double next_printtime = 100;
  while (f.time() <= ttot && !interrupt) {
    if (f.time() >= next_printtime) {
      next_printtime += 100;
      master_printf("Working on time %g...  ", f.time());
      master_printf("energy is %g\n", f.field_energy());
    }
\end{comment}
The monitor points are located one grid spacing from one another.  The
\verb*|get_new_point| method appends the fields at a given time to a
monitor point linked list.
\begin{verbatim}
    left  = f.get_new_point(vec(sourceloc+1.0/a), left );
    middle = f.get_new_point(vec(sourceloc+2.0/a), middle);
    right = f.get_new_point(vec(sourceloc+3.0/a), right);
\end{verbatim}
\begin{comment}
    f.step();
  }
  grace g("eps", dirname);
  complex<double> *al, *ar, *am, *freqs;
  int numl, numr;
  master_printf("Working on left fourier transform...\n");
\end{comment}
When the time stepping is over, we take a fourier transform of the fields
at the two monitor points.
\begin{verbatim}
  left->fourier_transform(Ex, &al, &freqs, &numl, 0.301, 0.5, 300);
\end{verbatim}
\begin{comment}
  delete[] freqs;
  master_printf("Working on middle fourier transform...\n");
  middle->fourier_transform(Ex, &am, &freqs, &numr, 0.301, 0.5, 300);
  delete[] freqs;
  master_printf("Working on right fourier transform...\n");
  right->fourier_transform(Ex, &ar, &freqs, &numr, 0.301, 0.5, 300);
  if (numl != numr)
     master_printf("Aaack you need both nums to be the same!\n");
  g.new_set();
  g.set_legend("\\x\\e\\s1\\N");
\end{comment}
Finally we calculate epsilon from the second derivative of the field using
\begin{equation*}
-k^2 H_z(\omega) = \nabla^2 H_z(\omega) = \nabla^2
\end{equation*}
\begin{verbatim}
  complex<double> *epsilon = new complex<double>[numl];
  for (int i=0;i<numl;i++) {
    complex<double> ksqr = -(ar[i]+al[i]-2.0*am[i])*a*a/am[i];
    epsilon[i] = ksqr/freqs[i]/freqs[i]/(2*pi*2*pi);
  } 
  for (int i=0;i<numl;i++)
    g.output_point(real(freqs[i]), real(epsilon[i]));
  g.new_set();
  g.set_legend("\\x\\e\\s2\\N");
  for (int i=0;i<numl;i++)
    g.output_point(real(freqs[i]), imag(epsilon[i]));
\end{verbatim}
\begin{comment}
  g.new_set();
  g.set_legend("analytic \\x\\e\\s1\\N");
  for (int i=0;i<numl;i++)
    g.output_point(real(freqs[i]),
                   real(f.analytic_epsilon(real(freqs[i]),vec(sourceloc+3.0/a))));
  g.new_set();
  g.set_legend("analytic \\x\\e\\s2\\N");
  for (int i=0;i<numl;i++)
    g.output_point(real(freqs[i]),
                   imag(f.analytic_epsilon(real(freqs[i]),vec(sourceloc+3.0/a))));
}
\end{comment}