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 [[email protected]]
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
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
