Transmission spectrum around a waveguide bend

#include <meep.hpp>
using namespace meep;

const double eps1 = 12.0;
const double eps2 = 1.0;
const double pml_thickness = 1.0;
const double x_center = 8.0;
const double y_center = 16.0;
const double sx = 2 * x_center;
const double sy = 2 * y_center;
const double pad = 4.0;
const double w = 1.0;
const double wvg_ycen = y_center - 0.5 * (sy - w - 2 * pad);
const double wvg_xcen = x_center + 0.5 * (sx - w - 2 * pad);
const double fcen = 0.15;
const double df = 0.1;
const double freq_max = fcen + 0.5 * df;
const double freq_min = fcen - 0.5 * df;
const double resolution = 10.0;
const double dT = 50.0;
const double decay_by = 0.001;
const int nfreq = 100;

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

double eps_no_bend (const vec &p)
{
  if (p.y() <= wvg_ycen + 0.5 * w && p.y() >= wvg_ycen - 0.5 * w) 
      return eps1;
  return eps2;
}

double eps_bend (const vec &p)
{
  if (p.y() <= wvg_ycen + 0.5 * w && p.y() >= wvg_ycen - 0.5 * w
      && p.x() >= 0.0 && p.x() <= sx - pad) 
      return eps1;
  if (p.x() <= wvg_xcen + 0.5 * w && p.x() >= wvg_xcen - 0.5 * w
      && p.y() <= sy && p.y() >= pad) 
      return eps1;
  return eps2;
}

void run_until_fields_decayed(fields &f, double dT, component c,
			      const vec &pt, double decay_by)
{
  double T0 = f.time();
  double max_abs = norm(f.get_field(c, pt));
  double cur_max = 0.0;
  double fa;
  bool stop = false;
  while (!stop) {
    while (f.time() < T0 + dT) {
      f.step();
      fa = norm(f.get_field(c, pt));
      cur_max = max(cur_max, fa);
    }
    T0 = f.time();
    max_abs = max(max_abs, cur_max);
    if (max_abs != 0.0) {
      master_printf("field decay(t = %g: %g / %g = %g\n",
		    f.time(), cur_max, max_abs, cur_max / max_abs);
    }
    if (cur_max <= max_abs * decay_by) stop = true;
    cur_max = 0.0;
  }
}
    

int main(int argc, char **argv)
{
  initialize mpi(argc, argv);
  const volume v = vol2d(2*x_center, 2*y_center, resolution);

  structure s0(v, eps_no_bend, pml(pml_thickness));
  s0.set_epsilon(eps_no_bend, true); // averaging
  fields f0(&s0);
  h5file *eps_no_bend = f0.open_h5file("eps_no_bend");
  f0.output_hdf5(Dielectric, v.surroundings(), eps_no_bend);
  delete eps_no_bend;
  
  structure s(v, eps_bend, pml(pml_thickness));
  s.set_epsilon(eps_bend, true); // averaging
  fields f(&s);
  h5file *eps_bend = f.open_h5file("eps_bend");
  f.output_hdf5(Dielectric, v.surroundings(), eps_bend);
  delete eps_bend;
  
  gaussian_src_time src(fcen, 1.0 / df, 0, 10.0 / df);
  geometric_volume src_plane(vec(1.0, wvg_ycen - 0.5*w),
			    vec(1.0, wvg_ycen + 0.5*w));
  f0.add_volume_source(Ez, src, src_plane);
  f.add_volume_source(Ez, src, src_plane);

  geometric_volume ft0_plane(vec(sx - 1.5, wvg_ycen - w),
			     vec(sx - 1.5, wvg_ycen + w));
  geometric_volume ft_plane(vec(wvg_xcen - w, sy - 1.5),
			    vec(wvg_xcen + w, sy - 1.5));
  geometric_volume fr_plane(vec(1.5, wvg_ycen - w),
			    vec(1.5, wvg_ycen + w));

  dft_flux ft0 = f0.add_dft_flux_plane(ft0_plane, freq_min, freq_max, nfreq);
  dft_flux fr0 = f0.add_dft_flux_plane(fr_plane, freq_min, freq_max, nfreq);
  dft_flux ft = f.add_dft_flux_plane(ft_plane, freq_min, freq_max, nfreq);
  dft_flux fr = f.add_dft_flux_plane(fr_plane, freq_min, freq_max, nfreq);

  run_until_fields_decayed(f0, dT, Ez, vec(sx-1.5, wvg_ycen), decay_by);
  fr -= fr0;
  run_until_fields_decayed(f, dT, Ez, vec(wvg_xcen, sy-1.5), decay_by);
  f.output_hdf5(Ez, v.surroundings());

  
  double *T = new double[nfreq];
  double *R = new double[nfreq];
  double *flux0 = ft0.flux();
  double *flux = ft.flux();
  for (int i = 0; i < nfreq; ++i)
    T[i] = flux[i] / flux0[i];
  flux = fr.flux();
  for (int i = 0; i < nfreq; ++i)
    R[i] = -flux[i] / flux0[i];

  for (int i = 0; i < nfreq; ++i)
    master_printf("transmission: %g, %g, %g, %g\n",
		  freq_min + i * df / (nfreq - 1), T[i], R[i],
		  1 - T[i] - R[i]);
  
  delete[] T;
  delete[] R;
  return 0;
}

pngファイルの作成

$ h5topng -S3 -Zc dkbluered -a yarg -A eps_bend.h5 -d ez.r ez-000300.00.h5

戻る