Material dispersion

#include <meep.hpp>
using namespace meep;

const double eps_inf = 2.25;
const double sigma1 = 0.5;
const double omega1 = 1.1;
const double gamma1 = 1e-5;
const double sigma2 = 2e-5;
const double omega2 = 0.5;
const double gamma2 = 0.1;

const double fcen = 1.0;
const double df = 2.0;
const double freq_max = fcen + 0.5 * df;
const double freq_min = fcen - 0.5 * df;
const double resolution = 20.0;

const int maxbands = 5;
const double harminv_time = 200.0;

const double kmin = 0.3;
const double kmax = 2.2;
const int k_num = 101;
const double dk = (kmax - kmin) / (k_num - 1);

inline double max(double a, double b) { return (a > b ? a : b); }
inline double sqr(double a) { return a*a; }

double feps (const vec &p)
{
  return eps_inf;
}

double fsigma1 (const vec &p)
{
  return sigma1;
}

double fsigma2 (const vec &p)
{
  return sigma2;
}

int main(int argc, char **argv)
{
  initialize mpi(argc, argv);
  const volume v = vol3d(0.0, 0.0, 0.0, resolution);
  structure s(v, feps);
  s.set_epsilon(feps, true); // averaging
  s.add_polarizability(fsigma1, omega1, gamma1);
  s.add_polarizability(fsigma2, omega2, gamma2); 
 
  complex<double> *amps = new complex<double>[maxbands];
  double *freq_re = new double[maxbands];
  double *freq_im = new double[maxbands];
  double *errors = new double[maxbands];
  for (int j = 0; j < k_num; j++) {
    double k = kmin + j * dk;
    fields f(&s);
    f.use_bloch(vec(k,0.0,0.0));
    gaussian_src_time src(fcen, 1.0 / df, 0.0, 10.0 / df);
    f.add_point_source(Ez, src, v.center());
    while (f.time() < f.last_source_time()) f.step();
    int ttot = int(harminv_time / f.dt)+1;
    complex<double> *mon_data = new complex<double>[ttot];
    int i = 0;
    while (f.time() < f.last_source_time() + harminv_time) {
      f.step();
      mon_data[i++] = f.get_field(Ez, v.center());
    }
    int nfreq = do_harminv(mon_data, ttot, f.dt, freq_min, freq_max,
			   maxbands, amps, freq_re, freq_im, errors); 
    for (int i = 0; i < nfreq; ++i) {
      complex<double> freq = complex<double>(freq_re[i], freq_im[i]);
      complex<double> e = pow(k / freq, 2);
      master_printf("eps:, %g, %g, %g, %g, %g\n",
		    k, freq_re[i], freq_im[i], real(e), imag(e));
    }
    delete[] mon_data;
  }

  delete[] errors;
  delete[] freq_im;
  delete[] freq_re;
  delete[] amps;
  return 0;
}