File: epsilon_polariton_1d.tex

package info (click to toggle)
meep 0.10-2.1
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 4,380 kB
  • ctags: 5,469
  • sloc: cpp: 50,653; sh: 8,380; haskell: 744; makefile: 367; perl: 10
file content (147 lines) | stat: -rw-r--r-- 4,702 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
\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}