Hi Geeth, Attached the code below. I added some random comments on the meep source code from meep.scm that probably only I can understand. Hope it helps.
Dear All, This is implementing Transmission spectrum example in "Meep Tutorial/Band diagram, resonant modes, and transmission in a holey waveguide". It might be of some use for you guys. I spent quite some time to figure this out. =================================================== #include "meep.hpp" using namespace meep; const double eps_=13 ; ////dielectric constant of waveguide const double w =1.2 ; ////width of waveguide const double r= 0.36 ; //radius of holes const double d= 1.4 ; //defect spacing (ordinary spacing = 1 const int N= 3 ; //number of holes on either side of defect //The cell dimensions const double sy =6 ; //size of cell in y direction (perpendicular to wvg. const double pad= 2 ; //padding between last hole and PML edge const double dpml= 1.0; //PML thickness const double sx= 2*(pad+dpml+N)+d-1; const double up_bound=sy/2+w/2,low_bound=sy/2-w/2;//up_bound, low_bound of wvg double hole_centers[2*N] ; const double fcen= 0.25; // pulse center frequency const double df =0.2 ; //pulse width (in frequency) const double nfreq = 500; // number of frequencies at which to compute flux double eps(const vec &p){ const grid_volume vol=vol2d(sx,sy,20);//20 is the resolution per unit distance vec r_vec; if(p.y()<low_bound || p.y()>up_bound) return 1.0; for(int i=0;i<2*N;i++){ r_vec=p-vec(hole_centers[i],sy/2);//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 1.0;//if the length of r < 1, it is within the hole } return eps_; } complex<double> one(const vec &p) { return 1.0; } int main(int argc,char **argv){ initialize mpi(argc,argv); double wg_center=sx/2; //define centers of holes for(int i=0;i<N;i++){ hole_centers[N+i]=wg_center+d/2+i*1; hole_centers[N-1-i]=wg_center-d/2-i*1; } const grid_volume vol=vol2d(sx,sy,20);//20 is the resolution per unit distance const symmetry S = -mirror(Y, vol);//odd symmetry through/in X axis, Y is the normal direction structure s(vol,eps,pml(dpml),S); fields f(&s); gaussian_src_time src(fcen,df); // src center: (center (+ dpml (* -0.5 sx)) 0) (size 0 w)))) //origin (0,0) is at center of waveguide for Scheme //now origin is at the left down corner //original src is (dpml-0.5sx,-w/2) to (dpml-0.5sx,w/2) //now it is (dpml,low_bound) to (dpml, up_bound) volume src_plane(vec(dpml,low_bound),vec(dpml,up_bound)); f.add_volume_source(Ey,src,src_plane,1.0); //(define trans ; transmitted flux // (add-flux fcen df nfreq // (make flux-region // (center (- (* 0.5 sx) dpml 0.5) 0) (size 0 (* w 2))))) //flux region in scheme (0.5sx-dpml-0.5,-w) to (0.5sx-dpml-0.5,w) //in C++ (sx-dpml-0.5,0) to (sx-dpml-0.5,2w) volume flux_plane(vec(sx-dpml-0.5,low_bound),vec(sx-dpml-0.5,up_bound)); dft_flux ft = f.add_dft_flux_plane(flux_plane,fcen-df,fcen+df,nfreq); f.output_hdf5(Dielectric, vol.surroundings());//output waveguide /****************** (run-sources+ (stop-when-fields-decayed 50 Ey (vector3 (- (* 0.5 sx) dpml 0.5) 0) 1e-3) (at-beginning output-epsilon) (during-sources (in-volume (volume (center 0 0) (size sx 0)) (to-appended "hz-slice" (at-every 0.4 output-hfield-z))))) refer to meep.scm the definition of stop-when-fields-decayed ********************/ /*************************** (define (stop-when-fields-decayed dT c pt decay-by) (if (null? fields) (init-fields)) (let ((T0 (meep-round-time));this assignment is done once, outside any loop, later will be modified (max-abs (sqr (magnitude (meep-fields-get-field fields c pt))));outside any loop (cur-max 0));outside any loop (lambda () (let ((fabs (sqr (magnitude (meep-fields-get-field fields c pt)))));outside any loop, once value is calculated, the syntax sugar/value for substitution never changes (set! cur-max (max cur-max fabs));in loop 1 "<= old-cur (* max-abs decay-by)", cos when the cond is not satisfied, it goes to the beg of the loop and executes set! (if (<= (meep-round-time) (+ T0 dT)) false ; don't stop yet, here is loop 2 for "<= (meep-round-time) (+ T0 dT)" (let ((old-cur cur-max));after loop 2 (set! cur-max 0);after loop 2 (set! T0 (meep-round-time));after loop 2 (set! max-abs (max max-abs old-cur));after loop 2 (if (not (zero? max-abs)) (print "field decay(t = " (meep-time)"): " old-cur " / " max-abs " = " (/ old-cur max-abs) "\n")) (<= old-cur (* max-abs decay-by))))))));cond for loop 1 ******************************/ double T0=f.round_time(),dT=50;//in time steps of 50 units double dt=0.4; double max_abs=sqrt(norm(f.get_field(Ey,vec(sx-dpml-0.5,sy/2)))); master_printf("T0:%g dT:%g max_abs:%g\n",T0,dT,max_abs); double cur_max=0,old_cur=0; double decay_by=.001; double f_abs; //slice (volume (center 0 0) (size sx 0)) in scheme for output of flux volume check_plane(vec(0,sy/2),vec(sx,sy/2)); h5file *myfile=f.open_h5file("test4.h5",h5file::WRITE,"test4",true); double previous_dT=f.round_time(),previous_dt=previous_dT; bool flag=true; do { //master_printf("cur_max:%g f_abs:%g\n",cur_max,f_abs); f.step(); if (f.round_time()>=previous_dT+dT) { f_abs=sqrt(norm(f.get_field(Ey,vec(sx-dpml-0.5,sy/2)))); cur_max=max(cur_max,f_abs); previous_dT=f.round_time(); flag=false; old_cur=cur_max; cur_max=0; T0=f.round_time(); max_abs=max(max_abs,old_cur); if(max_abs!=0) { master_printf("field decay(t = %g): %0.12g / %0.12g = %0.12g \n",f.time(),old_cur,max_abs,old_cur/max_abs); } } //master_printf("time:%g T0:%g\n",f.round_time(),T0); if( flag && f.round_time()>previous_dt+dt)//every 0.4 time units, output h5. { previous_dt=f.round_time(); f.output_hdf5(Hz,check_plane,myfile,true,false,"test4"); } }while(old_cur>=max_abs*decay_by ); delete myfile; //while (f.time() < nfreq / df / 2) f.step(); double *flux = ft.flux(); //Shall I normalize??????????? double dfreq=df/(nfreq-1); master_printf("tranmission:, freq, T[i]\n"); for (int i=0; i<nfreq; ++i) master_printf("transmission:, %0.12g, %0.12g\n",fcen-df/2+dfreq*i,flux[i]); delete [] flux; master_printf("finished.\n"); return 0; } ================================================ Best Regards, Bin Huang (Bryan) MIT Graduate Student '10 Computation for Design and Optimization (+65)98947649 ________________________________________ From: Geethaka Devendra [gcdeven...@gmail.com] Sent: Tuesday, February 08, 2011 8:55 PM To: Bin Huang Subject: Meep c++ code for source Hi Bryan, Hope you remember this posting from last year about a source and how to code it in c++. I have the same issue and I am totally lost and dont know how to use the "fields" class to do this job. I would appreciate some help from you and if could show me how you implemented it in your c++ code, that would be great. I am new to c++ and find it sometimes confusing how to use all the arguments in the functions. Cheers, Geeth (run-sources+ (stop-when-fields-decayed 50 Ey (vector3 (- (* 0.5 sx) dpml 0.5) 0) 1e-3) (during-sources (in-volume (volume (center 0 0) (size sx 0)) (to-appended "hz-slice" (at-every 0.4 output-hfield-z))))) It seems like usually the while loop is as follows: --------------------------- while ( f.time() < f.last_source_time()+400) { f.step(); f.output_hdf5(Ey, vol.surroundings()); } --------------------------- How do I increment the time in steps of 50 time units till Ey has decayed to 1/1000 on the plane indicated by "(vector3 (- (* 0.5 sx) dpml 0.5) 0)"? How do I output hfield on that slice indicated by "(volume (center 0 0) (size sx 0))" every 0.4 time unit? It seems like f.output_hdf5() can only output the field data in the entire waveguide. _______________________________________________ meep-discuss mailing list meep-discuss@ab-initio.mit.edu http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss