Dear Prof Johnson, (also Meepers) I am using the c++ interface of Meep to simulate a Bragg grating. The grating parameters : wavelength 1550nm, unit length a=1.0 micron, res=10, pixel width 0.1 micron, grating has 50 periods (to start with), Bragg period is 0.531496 micron. When I calculate the flux for transmission, I observe a shift in the Bragg wavelength, which is at 1570 nm. There is a 20nm shift. I have checked every aspect of the code and also tested the how the grating is drawn by Meep. The index of the grating is calculated using a sine function, which returns the index when the position is fed. When I plot the index of the grating after Meep draws it against length and also plot the sine function on top, I see an averaging. The index value Meep used is slightly different from the sine curve. I do not understand what Meep does here?
I have increased the sampling points, without increasing the resolution, of the grating in the hope that Meep will follow the sine and I will be able to see the index profile similar to the sine curve. But it did not do the job, When I increase the resolution, for example to 40, the Bragg transmission minimum comes back to 1550nm and it looks OK and I see the Bragg wavelength at 1550nm as expected for this structure. When the resolution is at 40, a pixel width is about 0.025 micron. when the resolution is at 10, a pixel width is 0.1 micron. Can this cause the problem, if YES, what is the relationship here? I have spent now about three months to figure out why this is there and how I could rectify this, but in vain! There should be a way to fix this. Can you help me and give me some advice to understand what is going on here. I have attached my c++ code here. Since you dont have the fdm code, just set the "multi_src" at the beginning of the code to "false", this will launch a gaussian point source. The code should work. I have also attached a plot of one period of the grating where you can see how Meep draws the period. I hope you will give me some advice, which will help me to get over this hump. Many thanks Geeth
<<attachment: index_x_a1.0_res10_x1_x10_x20_sinecurve_3.png>>
#include <iostream>
#include <fstream>
#include <meep.hpp>
#include <stdlib.h>
#include <iomanip>
using namespace meep;
const double a = 1.0 ; // unit length taken as 1 micron, so that all lengths are given as lengths of a
int res = 10 ; // sets number of pixels per unit length
int extra_t = 200; // simulation runtime after source turned off in meep units
// source details
double wlen = 1.55/a; // 1.55 micron
double src_df = 0.05; // delta_lambda = nm
component c = Ez; // Field component of the launching source
bool fluxoutput = false; // true outputs frequent flux data at count=10
bool multi_src = true;
// geometry parameters
double Neff = 1.458147488349648; // effective index obtained from mode solver
const double bragg_period = 1.550/Neff/2.0/a;
double dpml = 1.0/a, blank = 20.0/a;
const double period = (2.0*pi)/bragg_period; // period is the period of sine function of index
int n_periods = 50 ;
double grt_length = bragg_period * n_periods;
// setting waveguide and grating parameters
const double core_thick = 3.0/a;
double contrast = 0.2 ; // contrast clad-core
double clad_n = 1.45;
double clad_eps = clad_n * clad_n; // index 1.45
double core_n = clad_n * 1.01;
double core_eps = core_n * core_n;
double core_grt_n = (core_n * ((100.0+contrast)/100));
double clad_grt_n = (clad_n * ((100.0+contrast)/100));
double amp_core = core_grt_n - core_n; // index difference grating half period and core
double amp_clad = clad_grt_n - clad_n;
//double amp_core = 0.002929;
//double amp_clad = 0.0029;
// lattice dimensions
double lx = grt_length + 2.0*blank + 2.0*dpml;
double ly = core_thick * 10.0;
// grating location
double grt_start = (lx*0.5) - (grt_length*0.5);
double grt_end = (lx*0.5) + (grt_length*0.5);
// Source lcoation
//double src_x = 0.5* (lx*0.5 - grt_start) ;
double src_x = lx/2.0 - grt_length/2.0 - 18 ;
double src_y = ly * 0.5;
// geometry definition WG without grating - reference waveguide
double eps0(const vec &p) {
if (p.y() <= ly*0.5+core_thick*0.5 && p.y() >= ly*0.5-core_thick*0.5)
return core_eps;
else
return clad_eps;
}
// geometry definition WG with grating
double eps(const vec &p) {
// f = a * sin (b * x - c) + d; c=0
// double f_eps_core = amp_core * sin(period * p.x()) + core_eps;
// double f_eps_clad = amp_clad * sin(period * p.x()) + clad_eps;
double f_n_core = amp_core * sin(period * p.x()) + core_n;
double f_n_clad = amp_clad * sin(period * p.x()) + clad_n;
double f_eps_core = f_n_core * f_n_core;
double f_eps_clad = f_n_clad * f_n_clad;
if (p.y() < ly*0.5-core_thick*0.5 && p.x() > grt_start && p.x() < grt_end)
return f_eps_clad;
if (p.y() < ly*0.5-core_thick*0.5 && p.x() < grt_start && p.x() > 0)
return clad_eps;
if (p.y() < ly*0.5-core_thick*0.5 && p.x() > grt_end && p.x() < lx)
return clad_eps;
if (p.y() >= ly*0.5-core_thick*0.5 && p.y() <= ly*0.5+core_thick*0.5 &&
p.x() > grt_start && p.x() < grt_end )
return f_eps_core;
if (p.y() >= ly*0.5-core_thick*0.5 && p.y() <= ly*0.5+core_thick*0.5 &&
p.x() < grt_start && p.x() > 0 )
return core_eps;
if (p.y() >= ly*0.5-core_thick*0.5 && p.y() <= ly*0.5+core_thick*0.5 &&
p.x() > grt_end && p.x() < lx)
return core_eps;
if (p.y() > ly*0.5+core_thick*0.5 && p.x() > grt_start && p.x() < grt_end)
return f_eps_clad;
if (p.y() > ly*0.5+core_thick*0.5 && p.x() < grt_start && p.x() > 0)
return clad_eps;
return clad_eps;
}
// gaussian envelope function whic multiplies the gaussian source points\
// This is done to compare with the mode profile as source and the gaussian profile as source
double my_amp(double y) {
double L = 0.0809649, M = 3.73774;
double f = (L * exp (- (2 * ((y - ly/2)*(y - ly/2))/ (M*M))));
return f;
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
double resolution = res; // pixels per unit distance, a
grid_volume v = vol2d(lx,ly, resolution);
cout << "\n clad_n:" << sqrt(clad_eps) <<'\t' << "core_n :" << sqrt(core_eps) << "\n"
<< "core_grt_n :"<< core_grt_n << '\t' << "clad_grt_n :" << clad_grt_n <<
'\n' << "Bragg period :" << bragg_period << '\t' << "lx :" << lx <<"\n\n";
cout << "amp_core:" << amp_core << '\t' << "amp_clad: " << amp_clad << "\n\n" ;
// create two separate structures with (s) and without (s0) grating
structure s0(v, eps0, pml(dpml));
structure s(v, eps, pml(dpml));
fields my_field0(&s0); // field for structure without grating
fields my_field(&s); // field for structure with grating
h5file *epsfile_ref = my_field0.open_h5file("data/eps_ref", h5file::WRITE);
h5file *epsfile_grat = my_field.open_h5file("data/eps_grat", h5file::WRITE);
my_field0.output_hdf5(Dielectric, v.surroundings(), epsfile_ref );
my_field.output_hdf5(Dielectric, v.surroundings(), epsfile_grat);
// output the index along the x direction in the core
double *eps_core_x;
double *eps_clad_x;
double index_tot_cld = 0, index_tot_core =0 ;
double index_avg_core = 0, index_avg_cld =0;
int x_points = lx * res * 1 ; // number of pixels
cout << " lx : " << lx << '\t' << "x_points : " << x_points << '\n' ;
cout << "grt_length :" << grt_length << '\n' ;
cout << "grt_start :" << grt_start << '\t' << "grt_end :" << grt_end << '\n';
eps_core_x = new double [x_points+1];
eps_clad_x = new double [x_points+1];
int pix_count =0;
std::ofstream wgindexfile_x("./wgindex_x.dat");
std::ofstream grat_index_file("./grating.dat");
for (int i=1; i <= x_points; i++){
double delta_x = lx / (x_points);
double x_cord = delta_x * i;
//double x_cord = (delta_x* i) - delta_x*0.5;
if ((x_cord > grt_start) && (x_cord < grt_end))
{
eps_core_x[i] = s.get_eps(vec(x_cord, src_y)); // eps of core
eps_clad_x[i] = s.get_eps(vec(x_cord, (ly * 0.3)));
grat_index_file << x_cord << "\t" << sqrt(eps_core_x[i]) << "\t" << sqrt(eps_clad_x[i]) << "\n";
index_tot_core += sqrt(eps_core_x[i]);
index_tot_cld += sqrt(eps_clad_x[i]);
pix_count ++;
}
wgindexfile_x << x_cord << "\t" << sqrt(s.get_eps(vec(x_cord, src_y))) << "\t" << sqrt(s.get_eps(vec(x_cord, (ly * 0.3)))) << "\n";
index_avg_core = index_tot_core / pix_count;
index_avg_cld = index_tot_cld / pix_count ;
}
grat_index_file.close();
wgindexfile_x.close();
cout << "\nindex_avg_core :" << index_avg_core << '\t' << "index_avg_cld :" << index_avg_cld <<'\n';
cout << "pix_count :" << pix_count << '\n';
//************************************************************************************************
// insertion of mode solver as source
// For the mode solver to calcualte the supported fundamental mode of the waveguide, the location in a Y slice is
// extracted along with its's corresponding index value
// this information is then fed into a FDM mode solver to work out the mode, which is fed back to Meep as the
// source
ofstream wgindexfile_y;
wgindexfile_y.open("./wgindex_y.dat");
ofstream my_indexfile;
my_indexfile.open("./waveguide.index");
double fcen = 1.0 / wlen; double fwidth = src_df;
gaussian_src_time src(fcen, fwidth);
if (multi_src == true ) {
double *eps_y;
int y_points = ly * res * 100 ; //number of points along y direction
eps_y = new double [y_points+1];
double delta_y = ly / (y_points) ;
//cout << "src_gap :" << delta_y << '\n' ;
for (int i=1; i <= y_points; i++){
double src_y_cord = i * delta_y;
eps_y[i] = s.get_eps(vec(src_x, (src_y_cord))); // eps of ref structure (s0)
my_indexfile << delta_y << '\t' << sqrt(eps_y[i]) << '\n' ;
wgindexfile_y << src_y_cord << '\t' << sqrt(eps_y[i]) << '\n' ;
}
my_indexfile.close();
wgindexfile_y.close();
// exit(0);
master_printf("\n...guided mode launched as source...\n\n");
//FD program input parameters
ofstream para_file("./fd_params.dat");
para_file << "E\n" //TE mode
<< "1.55\n" //wavelength
<< y_points + 1 << "\n"
<< 1 << "\n" //search for 1 mode
<< core_n << "\n"; //guess
para_file.close();
system("1D_fd_waveguide waveguide.index < fd_params.dat");
//Now get the mode profile
ifstream mode_file("./eigenvectors.out");
double y, mode, dummy;
ofstream out_file("mode.dat");
for (int i=1; i <= y_points; i++){
mode_file >> y >> mode >> dummy;
y*=1e6;
my_field0.add_point_source(c, src, vec(src_x, y), mode);
my_field.add_point_source(c, src, vec(src_x, y), mode);
//my_field0.add_point_source(c, src, vec(src_x, y), my_amp(y));
//my_field.add_point_source(c, src, vec(src_x, y), my_amp(y));
out_file << y << "\t" << mode << "\t" << my_amp(y) << "\n";
}
out_file.close();
exit (0);
}
else {
// only a single source point launched in the waveguide
my_field0.add_point_source(c, src, vec(src_x, src_y));
my_field.add_point_source(c, src, vec(src_x, src_y));
master_printf("\n...single point source added...\n\n");
}
//**************************************************************************************************
// define coordinates for flux planes. Ref flux plane and Reflection flux plane
// both at same X-location. They are run successively, so wont be a problem
//double rflux_x = lx/2 - grt_length/2 - blank/2;
//double trflux_x = lx/2.0 + grt_length/2.0 + blank*0.75 ;
double trflux_x = lx/2.0 + grt_length/2.0 + blank * 0.5 ;
double rflux_x = lx/2.0 - grt_length/2.0 - blank * 0.5 ;
double flux_y1 = ly/2.0 + core_thick/2.0 ;
double flux_y2 = ly/2.0 - core_thick/2.0 ;
// definition of flux planes & my objects for calling function in class "dft_flux"
int nfreq = 500;
double freq_min = fcen - fwidth*0.5;
double freq_max = fcen + fwidth*0.5;
volume ref_flux_plane(vec(rflux_x, flux_y1), vec(rflux_x, flux_y2));
volume tr_ref_flux_plane(vec(trflux_x, flux_y1), vec(trflux_x, flux_y2));
volume trans_flux_plane(vec(trflux_x, flux_y1), vec(trflux_x, flux_y2));
volume refl_flux_plane(vec(rflux_x, flux_y1), vec(rflux_x, flux_y2));
dft_flux my_ref_flux = my_field0.add_dft_flux_plane(ref_flux_plane, freq_min, freq_max, nfreq);
dft_flux my_tr_ref_flux = my_field0.add_dft_flux_plane(tr_ref_flux_plane, freq_min, freq_max, nfreq);
dft_flux my_trans_flux = my_field.add_dft_flux_plane(trans_flux_plane, freq_min, freq_max, nfreq);
dft_flux my_refl_flux = my_field.add_dft_flux_plane(refl_flux_plane, freq_min, freq_max, nfreq);
// simulate WG structure without grating, time stepping
master_printf("\n...simulating waveguide without grating...\n\n");
while (my_field0.time() < my_field0.last_source_time()+ extra_t)
{
my_field0.step();
}
h5file *fieldfile_ref = my_field0.open_h5file("data/h5files/ez_ref", h5file::WRITE);
my_field0.output_hdf5(c, v.surroundings(), fieldfile_ref);
double *ref_flux_data = my_ref_flux.flux();
double *tr_ref_flux_data = my_tr_ref_flux.flux();
// simulate WG structure with grating, time stepping
master_printf("\n...simulating waveguide with grating structure...\n\n");
//h5file *fieldfile_grat = my_field.open_h5file("/data/ez_grat", h5file::WRITE);
if (fluxoutput == true) {
int count = 0;
int i = 1;
ostringstream ss;
while (my_field.time() < my_field.last_source_time()+ extra_t)
{
my_field.step();
if (count == 10)
{
// field filename generating with a count every 10 steps
ss << "/data/h5files/ez_grat_" << setw(5) << setfill('0') << int(i);
string s = ss.str();
int len = s.length();
char *filename = new char [len+1];
s.copy(filename, len , 0);
filename[len]=0;
// output field components
h5file *fieldfile_grat = my_field.open_h5file(filename, h5file::WRITE);
my_field.output_hdf5(c, v.surroundings(), fieldfile_grat);
count = 0;
i ++;
ss.str("");
}
count ++;
}}
else
{
while (my_field.time() < my_field.last_source_time()+ extra_t)
{
my_field.step();
}
h5file *fieldfile_grat = my_field.open_h5file("/data/h5files/ez_grat_once", h5file::WRITE);
my_field.output_hdf5(c, v.surroundings(), fieldfile_grat);
}
// record flux at flux planes defined before
my_refl_flux -= my_ref_flux; // deduct the incident field from source before normalising
double *refl_flux_data = my_refl_flux.flux();
double *trans_flux_data = my_trans_flux.flux();
double dfreq = src_df / (nfreq -1);
// normalise flux data and output to file
double *trans, *refl, *test_ref;
trans = new double [nfreq];
refl = new double [nfreq];
// test_ref = new double [nfreq];
// write the flux data into a file
ofstream my_datafile;
my_datafile.open("./bragg_wg.dat");
for (int i=0; i < nfreq; i++) {
refl[i] = refl_flux_data[i] / ref_flux_data[i];
trans[i] = trans_flux_data[i] / tr_ref_flux_data[i];
my_datafile << my_refl_flux.freq_min + my_refl_flux.dfreq*i << '\t' << ref_flux_data[i] <<
'\t'<< refl[i] << '\t' << trans[i] << '\t' << refl_flux_data[i] << '\t' <<
trans_flux_data[i] << '\n';
}
my_datafile.close();
master_printf("\n...flux data written to file...\n");
// deallocating memory for these
delete [] trans_flux_data;
delete [] refl_flux_data;
delete [] ref_flux_data;
delete [] tr_ref_flux_data;
return 0;
}
_______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
