The following code provides a work around using 3d such that Ey and Ex behave as per the 1d Ex case for frequency and amplitude, albeit with a slightly lower quality factor. The Ez case is also run to demonstrate that it completes execution.
#include <meep.hpp> using namespace meep; //with changes from http://ab-initio.mit.edu/wiki/index.php/Meep_C-plus-plus_Tutorial#Calculatin g_the_flux_through_a_close_packed_opal_monolayer const double eps1 = 12.0; // epsilon of layer 1 const double eps2 = 1.0; // epsilon of layer 2 const double grating_periodicity = 1.0; const double d1 = sqrt(eps2) / (sqrt(eps1)+sqrt(eps2)); // quarter wave stack dimensions (0.2240) const double d2 = sqrt(eps1) / (sqrt(eps1)+sqrt(eps2)); // (0.7760) const double half_cavity_width = d2; const int N = 5; const double pml_thickness = 1.0; const double z_center = half_cavity_width + N*grating_periodicity + pml_thickness; // (6.7760) complex<double> one(const vec &p) { return 1.0; } double eps(const vec &p) { vec r = p - vec(z_center); if (abs(r.z()) < half_cavity_width) return 1.0; else { double dz = abs(r.z()) - half_cavity_width; while (dz > grating_periodicity) dz -= grating_periodicity; return (dz < d1) ? eps1 : eps2; } } void excitefilm(component c) { const double amicron = 10.0; const volume vol = vol3d(0.6,1.0,2*z_center,amicron); structure s(vol,eps,pml(pml_thickness,Z)); fields f(&s); f.use_bloch(vec(0,0,0)); const double w_midgap = (sqrt(eps1)+sqrt(eps2))/(4*sqrt(eps1)*sqrt(eps2)); //0.3222 gaussian_src_time src(w_midgap, 1.0, 0.0, 5.0); geometric_volume src_plane(vec(0,0,z_center),vec(0.6,1.0,z_center)); f.add_volume_source(c,src,src_plane,one,1.0); f.output_hdf5(Dielectric, vol.surroundings()); while (f.time() < f.last_source_time()) f.step(); int maxbands = 5; int ttot = int(400/f.dt)+1; complex<double> *p = new complex<double>[ttot]; complex<double> *amps = new complex<double>[maxbands]; double *freq_re = new double[maxbands]; double *freq_im = new double[maxbands]; int i=0; while (f.time() < f.last_source_time() + 400) { f.step(); p[i++] = f.get_field(c,vol.center()); } int num = do_harminv(p, ttot, f.dt, 0.8*w_midgap, 1.2*w_midgap, maxbands, amps, freq_re, freq_im); master_printf("frequency,amplitude,quality factor\n"); for (int i=0; i<num; ++i) master_printf("%0.6g%+0.6gi,%0.5g%+0.5gi,%g\n",freq_re[i],freq_im[i],real(am ps[i]),imag(amps[i]),-0.5*freq_re[i]/freq_im[i]); double curr_time = f.time(); while (f.time() < curr_time + 1/w_midgap) { f.output_hdf5(c, vol.surroundings()); f.step(); } } int main(int argc, char **argv) { initialize mpi(argc,argv); master_printf("myfileresonator3d.cpp starting.\n"); excitefilm(Ex); excitefilm(Ey); excitefilm(Ez); master_printf("myfileresonator3d.cpp finished.\n"); return 0; } -----Original Message----- From: [email protected] [mailto:[email protected]] On Behalf Of Kristian Medri Sent: February 8, 2010 6:50 PM To: [email protected] Subject: [Meep-discuss] 1D Fabry-Perot Dear Meep community, I would appreciate some input as to why the 1D Fabry-Perot only works with Ex. The code below is for the example from: http://ab-initio.mit.edu/wiki/index.php/Meep_C-plus-plus_Tutorial The tutorial states that "For a 1D system, Meep considers a computational cell along the z coordinate." Thus why does running a simulation with Ey behave differently from running it with Ex? ------------ #include <meep.hpp> using namespace meep; //if Ey instead of Ex - meep: cannot require a ey component in a 1D grid //if Ez instead of Ex - meep: cannot require a ez component in a 1D grid const double eps1 = 12.0; // epsilon of layer 1 const double eps2 = 1.0; // epsilon of layer 2 const double grating_periodicity = 1.0; const double d1 = sqrt(eps2) / (sqrt(eps1)+sqrt(eps2)); // quarter wave stack dimensions const double d2 = sqrt(eps1) / (sqrt(eps1)+sqrt(eps2)); const double half_cavity_width = d2; const int N = 5; const double pml_thickness = 1.0; const double z_center = half_cavity_width + N*grating_periodicity + pml_thickness; double eps(const vec &p) { vec r = p - vec(z_center); if (abs(r.z()) < half_cavity_width) return 1.0; else { double dz = abs(r.z()) - half_cavity_width; while (dz > grating_periodicity) dz -= grating_periodicity; return (dz < d1) ? eps1 : eps2; } } int main(int argc, char **argv){ initialize mpi(argc,argv); const double amicron = 10.0; const volume vol = vol1d(2*z_center,amicron); structure s(vol,eps,pml(pml_thickness)); fields f(&s); const double w_midgap = (sqrt(eps1)+sqrt(eps2))/(4*sqrt(eps1)*sqrt(eps2)); gaussian_src_time src(w_midgap, 1.0, 0.0, 5.0); f.add_point_source(Ex, src, vol.center()); f.output_hdf5(Dielectric, vol.surroundings()); while (f.time() < f.last_source_time()) f.step(); int maxbands = 5; int ttot = int(400/f.dt)+1; complex<double> *p = new complex<double>[ttot]; complex<double> *amps = new complex<double>[maxbands]; double *freq_re = new double[maxbands]; double *freq_im = new double[maxbands]; int i=0; while (f.time() < f.last_source_time() + 400) { f.step(); p[i++] = f.get_field(Ex,vol.center()); } int num = do_harminv(p, ttot, f.dt, 0.8*w_midgap, 1.2*w_midgap, maxbands, amps, freq_re, freq_im); master_printf("frequency,amplitude,quality factor\n"); for (int i=0; i<num; ++i) master_printf("%0.6g%+0.6gi,%0.5g%+0.5gi,%g\n",freq_re[i],freq_im[i],real(am ps[i]),imag(amps[i]),-0.5*freq_re[i]/freq_im[i]); double curr_time = f.time(); while (f.time() < curr_time + 1/w_midgap) { f.output_hdf5(Ex, vol.surroundings()); f.step(); } return 0; } _______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss _______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
