## File: epsilon_polariton_1d.tex

package info (click to toggle)
meep 0.10-2.1
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 \begin{comment} #include #include #include #include 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 *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 *epsilon = new complex[numl]; for (int i=0;i 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