Dear meep users
I am trying to reproduce transmission in fig(4)of chapter 10 in book
Photonic crystals molding the flow of light second edition.But I am
getting peak transmission 0.09.I have used plane wave source and line
detector for the simulation and run until field decayed(0.001). I have
attached the c++ code also.Please indicate what could be the error in the
code?
#include "meep.hpp"
#include <iostream>
#include <fstream>
using namespace meep;
const double eps_=11.4 ; ////dielectric constant of waveguide
const double w =11.0; ////width of waveguide
const double r= 0.20 ; //radius of holes
const double d= 1.0 ; //defect spacing (ordinary spacing = 1
const int N1= 17 ; //number of holes on either side of defect
const int N2= 11 ; //number of holes on either side of defect
const double pad= 2.0 ; //padding between last hole and PML edge
const double dpml= 5.0; //PML thickness
//The cell dimensions
const double sy =w; ; //size of cell in y direction (perpendicular to wvg.
const double sx= 2*(pad+dpml)+17;
const double sx0= 2*(pad+dpml)+17;
const double up_bound=sy/2+w/2,low_bound=sy/2-w/2;//up_bound, low_bound of
wvg
double hole_centerx[N1] ;
double hole_centery[N2] ;
const double f_max=0.48;
const double f_min=0.28;
const double fcen=(f_max+f_min)*0.5;
const double df=(f_max-f_min);
const int nfreq=2001;
const double dT = 50.0;
const double decay_by = 0.001;
inline double max(double a, double b) { return (a > b ? a : b); }
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;
}
}
double eps(const vec &p){
for(int i=0;i<N1;i++)
hole_centerx[i]=dpml+pad+0.5+1*i;
for(int i=0;i<N2;i++)
hole_centery[i]=0.5+1*i;
vec r_vec;
if(p.y()<low_bound || p.y()>up_bound)
return 1.0;
for(int j=0;j<N2;j++)
{
for(int i=0;i<N1;i++)
{
if(j!=5)
{
r_vec=p-vec(hole_centerx[i],hole_centery[j]);//compare to one hole
center
//master_printf("x:%g y:%g r:%g",r.x(),r.y(),r.r());
if (r_vec.x()*r_vec.x()+r_vec.y()*r_vec.y()<=r*r)
return eps_;//if the length of r < 1, it is within the hole
}
if(j==5)
{
if(i==6 || i==7 || i==9 || i==10)
{
vec r_vec1;
r_vec1=p-vec(hole_centerx[i],hole_centery[j]);//compare to one hole center
//master_printf("x:%g y:%g r:%g",r.x(),r.y(),r.r());
if (r_vec1.x()*r_vec1.x()+r_vec1.y()*r_vec1.y()<=r*r)
return eps_;//if the length of r < 1, it is within the hole
}
}
}
}
return 1.0 ;
}
double eps0(const vec &){
return 1.0 ;
}
complex<double> one(const vec &p) {
return 1.0;
}
int main(int argc,char **argv)
{
initialize mpi(argc,argv);
const grid_volume v0=vol2d(sx0,sy,20);
const grid_volume v=vol2d(sx,sy,20);
const symmetry S = mirror(Y, v);//odd symmetry through/in X axis, Y is the
normal direction
const symmetry S0 = mirror(Y, v0);//odd symmetry through/in X axis, Y is
the normal direction
structure s0(v0,eps0,pml(dpml,X),S0);
structure s(v,eps,pml(dpml,X),S);
fields f0(&s0);
f0.use_bloch(Y,0);
//f0.use_real_fields();
//f0.output_hdf5(Dielectric,v.surroundings());
fields f(&s);
f.use_bloch(Y,0);
//f.use_real_fields();
f.output_hdf5(Dielectric,v.surroundings());
gaussian_src_time src(fcen,df);
volume src_plane(vec(dpml+0.5,0),vec(dpml+0.5,sy));
f0.add_volume_source(Ez,src,src_plane,one,1.0);
f.add_volume_source(Ez,src,src_plane,one,1.0);
volume rflux_plane(vec(dpml+0.5,0),vec(dpml+0.5,sy));
volume tflux_plane(vec(sx-dpml-1.5,0),vec(sx-dpml-1.5,sy));
volume t0flux_plane(vec(sx0-dpml-1.5,0),vec(sx0-dpml-1.5,sy));
dft_flux ft0 = f0.add_dft_flux_plane(t0flux_plane,f_min,f_max,nfreq);
dft_flux ft = f.add_dft_flux_plane(tflux_plane,f_min,f_max,nfreq);
dft_flux fr0 = f0.add_dft_flux_plane(rflux_plane,f_min,f_max,nfreq);
dft_flux fr = f.add_dft_flux_plane(rflux_plane,f_min,f_max,nfreq);
run_until_fields_decayed(f0, dT, Ez, vec(sx0-dpml-1.5, sy/2), decay_by);
fr -= fr0;
run_until_fields_decayed(f, dT, Ez, vec(sx-dpml-1.5, sy/2), decay_by);
double dfreq=df/(nfreq-1);
double *flux = ft.flux();
double *flux0 = ft0.flux();
double *T,*R;
T=new double[nfreq];
R=new double[nfreq];
for (int i = 0; i < nfreq; ++i)
T[i] = flux[i] / flux0[i];
delete[] flux;
flux = fr.flux();
for (int i = 0; i < nfreq; ++i)
R[i] = -flux[i] / flux0[i];
ofstream myfile;
myfile.open ("test.dat");
for (int i=0; i<nfreq; ++i)
myfile << "transmission:,"<<f_min+dfreq*i<<","<<T[i]<<","<<R[i]<<"\n";
myfile.close();
delete[] flux;
delete[] flux0;
delete[] T;
delete[] R;
f.output_hdf5(Ez, v.surroundings());
return 0;
}
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
