\section{Nonlinear materials}

FIXME: Add a nice discussion of nonlinear materials here...

In this example, we will use a CW source and compute the amplitude of the
field at a given position and time as a function of the source amplitude.
The result will be linear as long as the material remains in the linear
regime.  Once we have departed from the linear regime, we get more
complicated behavior.  Yes, this is a stupid example...

\begin{figure}
\label{nonlinear_field_shape}
\caption{Field vs. source amplitude with nonlinear material}
\includegraphics[width=14cm,clip=true]{nonlinear-out/ez-000400-00}
\end{figure}

The system is shown in Figure~\ref{nonlinear_field_shape}.  It is a 2D
metallic waveguide with vacuum $\epsilon$ of 2.25.  In the center of the
cell is a small region that is linear which contains the source, and the
rest of the waveguide contains a nonlinear material.  Both ends of the
waveguide have PML absorbing boundary conditions.

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

#include <meep.hpp>
using namespace meep;

const double a = 10.0;
const int m = 0;
const double eps_value = 2.25;

const double pml_thickness = 2.0;
const double xmax = 60.0;
const double ymax = 1.0;

double eps(const vec &v) {
  return eps_value;
}
\end{comment}
\begin{verbatim}
const double alpha_value = 0.07;
double alpha(const vec &v) {
  if (fabs(v.x() - xmax/2.0) < .51) return 0.0;
  return alpha_value;
}
\end{verbatim}

\begin{comment}
int main(int argc, char **argv) {
  initialize mpi(argc, argv);
  deal_with_ctrl_c();
  const char *dirname = make_output_directory(__FILE__);
  volume v = vol2d(xmax, ymax, a);
  const symmetry S = mirror(X,v) + mirror(Y,v);
  structure s(v, eps, pml(pml_thickness, X), S);
  s.set_output_directory(dirname);
  grace g("field", dirname);
  g.output_point(0.0, 0.0);
  g.output_point(0.0, 0.0);
\end{comment}

The set\_chi3 method is used to set the chi3 coefficient.

\begin{verbatim}
  s.set_chi3(alpha);
\end{verbatim}

We use a CW source at a frequency of 0.4, which gives single-mode behavior
when the amplitude is small.  Also note that we use real fields, since
complex fields are incorrect for nonlinear materials.

\begin{verbatim}
  continuous_src_time my_source(0.4,  0.8);
  for (double amp = 0.05; amp <= 1.01 && !interrupt; amp += 0.01) {
    fields f(&s, m);
    f.use_real_fields();
    f.add_point_source(Ez, my_source, vec(xmax*0.5,ymax*0.5), amp);
\end{verbatim}

\begin{figure}
\label{nonlinear_field}
\caption{Field vs. source amplitude with nonlinear material}
\includegraphics[width=8.8cm,clip=true]{nonlinear-out/field}
\end{figure}

Time stepping, etc, is done as usual.  We monitor the field at one end of
the cell, which gives Figure~\ref{nonlinear_field}, which show the field
versus source amplitude.

\begin{comment}
    master_printf("Working with A=%g...\n", amp);
    double next_printtime = 400;
    while (f.time() < 10.1*xmax && !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());
	f.output_hdf5(Ez, f.total_volume());
      }
      f.step();
    }
    g.output_point(amp, real(f.get_field(Ez,vec(1.2*pml_thickness,ymax/2))));
  }
}
\end{comment}
