On Mon, 1 May 2006, Jun She wrote:
I wrote, f.initialize_with_n_tm(4). Then I compiled the code and run the
program. I saw, "meep: cannot require a hp compoint in a 2D grid"
This function is for cylindrical coordinates only, and initializes the
fields with a particular mode of a cylindrical metal waveguide. It's not
what you want for your problem. (It's one of several functions in the C++
interface that is fairly old problem-specific cruft and should arguably be
removed.)
In general, you don't need to initialize the field at all. To compute
band structures, just excite the field (initialized to zero, the default)
with a broadband pulse and then analyze the result with harminv.
I don't know why the field need to be intialized twice?
It's not being initialized twice. It's initializing two different
components of the field. And the fields:: function calls the
fields_chunk:: function to initialize the different "chunks" of the fields
(the grid is broken up into a set of chunks corresponding to different
CPUs, different update equations, etcetera).
Is require_componet(Ez) for 1d,2d,3d coordinates and require_component(Hp)
for cylindrical coordinates?
No. require_component just makes sure that those components are
allocated before they are initialized.
Why not require_component(Hr) here?
Because that component is not initialized here. (In any case, requiring
some components automatically implies that other related components are
allocated too.)
Could I ask how to compute a band structure using C++?
In general, I recommend using the libctl/Scheme interface, since its usage
is better documented. See the example bandstructure calculation in:
http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial/Band_diagram%2C_resonant_modes%2C_and_transmission_in_a_holey_waveguide
If you prefer to use C++, you'll want to use basically the same technique,
just with more manual coding. However, for your information I have
attached a C++ file that I once used to compute the band diagram of a
triangular lattice of metallic rods in Meep.
Steven
#include <math.h>
#include <getopt.h>
#include <stdlib.h>
#include <time.h>
#include <meep.hpp>
using namespace meep;
// parameters:
double eps_rod = -1e12; // defaults to "metal"
double r = 0.47; // radius of metal rods
/* we will make a triangular lattice using a 1 by sqrt(3) supercell */
double eps(const vec &v) {
double x = v.x(), y = v.y();
if (hypot(x, y) < r
|| hypot(x - 1, y) < r
|| hypot(x, y - sqrt(3.)) < r
|| hypot(x - 1, y - sqrt(3.)) < r
|| hypot(x - 0.5, y - 0.5 * sqrt(3.)) < r)
return eps_rod;
return 1.0;
}
double urand(double a, double b)
{
return a + rand() * (b - a) / RAND_MAX;
}
#define NUM_K 62
const double kpoints[NUM_K][2] = { /* from MPB (cartesian, w/o 2pi) */
{0.0, 0.0},
{0.0, 0.026243194054074},
{0.0, 0.0524863881081479},
{0.0, 0.0787295821622218},
{0.0, 0.104972776216296},
{0.0, 0.13121597027037},
{0.0, 0.157459164324444},
{0.0, 0.183702358378517},
{0.0, 0.209945552432591},
{0.0, 0.236188746486665},
{0.0, 0.262431940540739},
{0.0, 0.288675134594813},
{0.0, 0.314918328648887},
{0.0, 0.341161522702961},
{0.0, 0.367404716757035},
{0.0, 0.393647910811109},
{0.0, 0.419891104865182},
{0.0, 0.446134298919256},
{0.0, 0.47237749297333},
{0.0, 0.498620687027404},
{0.0, 0.524863881081478},
{0.0, 0.551107075135552},
{0.0, 0.577350269189626},
{-0.0256410256410256, 0.577350269189626},
{-0.0512820512820513, 0.577350269189626},
{-0.0769230769230769, 0.577350269189626},
{-0.102564102564103, 0.577350269189626},
{-0.128205128205128, 0.577350269189626},
{-0.153846153846154, 0.577350269189626},
{-0.179487179487179, 0.577350269189626},
{-0.205128205128205, 0.577350269189626},
{-0.230769230769231, 0.577350269189626},
{-0.256410256410256, 0.577350269189626},
{-0.282051282051282, 0.577350269189626},
{-0.307692307692308, 0.577350269189626},
{-0.333333333333333, 0.577350269189626},
{-0.320512820512821, 0.555144489605409},
{-0.307692307692308, 0.532938710021193},
{-0.294871794871795, 0.510732930436977},
{-0.282051282051282, 0.48852715085276},
{-0.269230769230769, 0.466321371268544},
{-0.256410256410256, 0.444115591684328},
{-0.243589743589744, 0.421909812100111},
{-0.230769230769231, 0.399704032515895},
{-0.217948717948718, 0.377498252931678},
{-0.205128205128205, 0.355292473347462},
{-0.192307692307692, 0.333086693763246},
{-0.17948717948718, 0.310880914179029},
{-0.166666666666667, 0.288675134594813},
{-0.153846153846154, 0.266469355010596},
{-0.141025641025641, 0.24426357542638},
{-0.128205128205128, 0.222057795842164},
{-0.115384615384615, 0.199852016257947},
{-0.102564102564103, 0.177646236673731},
{-0.0897435897435897, 0.155440457089515},
{-0.0769230769230769, 0.133234677505298},
{-0.0641025641025641, 0.111028897921082},
{-0.0512820512820513, 0.0888231183368655},
{-0.0384615384615385, 0.0666173387526491},
{-0.0256410256410256, 0.0444115591684328},
{-0.0128205128205128, 0.0222057795842164},
{0.0, 0.0}
};
int main(int argc, char **argv)
{
initialize mpi(argc, argv);
double resolution = 40; // pixels/a
bool hdf_output = false;
bool te = false;
double T = 100;
double freq = 0.3;
double dfreq = 0.5; // fractional
int which_k = -1; // default = all
srand(314159); // deterministic "random" nums
// srand(broadcast(0,time(NULL)));
int ch;
while ((ch = getopt(argc, argv, "a:hr:HAT:f:d:e:k:")) != -1)
switch (ch) {
case 'a':
resolution = atof(optarg);
break;
case 'T':
T = atof(optarg);
break;
case 'k':
which_k = atoi(optarg);
break;
case 'H':
hdf_output = true;
break;
case 'h':
te = true;
break;
case 'r':
r = atof(optarg);
break;
case 'e':
eps_rod = atof(optarg);
break;
case 'f':
freq = atof(optarg);
break;
case 'd':
dfreq = atof(optarg);
break;
default:
meep::abort("unrecognized argument %c", (char) ch);
}
volume v = vol2d(1.0, sqrt(3.), resolution);
structure s(v, eps);
if (hdf_output) {
trash_output_directory("trirods-out");
s.set_output_directory("trirods-out");
}
fields f(&s);
f.use_bloch(vec(0,0));
if (hdf_output) {
f.output_hdf5(Dielectric, geometric_volume(vec(-1.0,-sqrt(3.0)),
vec(2.0, 2*sqrt(3.0))));
master_printf("Output dielectric function.\n");
}
double freq_min = freq * (1 - dfreq);
double freq_max = freq * (1 + dfreq);
gaussian_src_time src((freq_min + freq_max) * 0.5,
0.5 / fabs(freq_max - freq_min),
0, 5 / fabs(freq_max - freq_min));
int Ncpt = int(T / f.dt + 0.5) + 1;
complex<double> *cpt = new complex<double>[Ncpt];
complex<double> *cpt2 = new complex<double>[Ncpt];
for (int kindex = 0; kindex < NUM_K; ++kindex) {
if (which_k > 0) kindex = which_k;
f.zero_fields();
f.remove_sources();
f.t = 0;
double kx = kpoints[kindex][0], ky = kpoints[kindex][1];
f.use_bloch(X, kx);
f.use_bloch(Y, ky);
// add sources, including Bloch-periodic replication of the primitive cell
double xsrc, ysrc; // random src location...make sure is in air!
do { xsrc = urand(0,1); ysrc = urand(0,sqrt(3.));
} while (eps(vec(xsrc,ysrc)) != 1.0);
double xsrc2 = xsrc, ysrc2 = ysrc; // equiv. src. location (by periodicity)
for (int i = -5; i < 5; ++i)
for (int j = -5; j < 5; ++j) {
double xsrc1 = xsrc + i * 1 + j * 0.5;
double ysrc1 = ysrc + i * 0 + j * 0.5 * sqrt(3.);
if ((i != 0 || j != 0)
&& xsrc1 >= 0 && xsrc1 <= 1.0
&& ysrc1 >= 0 && ysrc1 <= sqrt(3.))
f.add_point_source(Ez, src, vec(xsrc2 = xsrc1, ysrc2 = ysrc1),
polar(1.0, 2*pi*(kx * xsrc1 + ky * ysrc1)));
}
master_printf("src locations: (%g,%g) and (%g,%g)\n",
xsrc, ysrc, xsrc2, ysrc2);
int NT = int((f.last_source_time() + T) / f.dt + 0.5) + 1;
while (f.time() < f.last_source_time()) {
f.step();
if (f.t % 100 == 0)
master_printf("time step %d, %g%% done\n", f.t, f.t * 100.0 / NT);
}
int i = 0;
while (f.time() < f.last_source_time() + T) {
cpt[i] = f.get_field(Ez, vec(xsrc,ysrc));
cpt2[i] = f.get_field(Ez, vec(xsrc2,ysrc2));
i++;
f.step();
if (f.t % 100 == 0)
master_printf("time step %d, %g%% done\n", f.t, f.t * 100.0 / NT);
}
complex<double> amp[100];
double freq_re[100], freq_im[100], freq_err[100];
int num_freqs = do_harminv(cpt, i, f.dt, freq_min/4, freq_max*4, 100,
amp, freq_re, freq_im, freq_err);
complex<double> amp2[100];
double freq_re2[100], freq_im2[100], freq_err2[100];
int num_freqs2 = do_harminv(cpt2, i, f.dt, freq_min/4, freq_max*4, 100,
amp2, freq_re2, freq_im2, freq_err2);
for (int i = 0; i < num_freqs; ++i) {
master_printf("freq[%d] = %g%+gi +/- %g (amp %g%+gi)\n",
i, freq_re[i], freq_im[i], freq_err[i],
real(amp[i]), imag(amp[i]));
int i2 = 0;
for (int j = 0; j < num_freqs2; ++j)
if (fabs(freq_re[i] - freq_re2[j]) < fabs(freq_re[i] - freq_re2[i2]))
i2 = j;
if (fabs(freq_re[i] - freq_re2[i2]) < 1e-4) {
master_printf(" = freq2[%d] = %g%+gi +/- %g (amp %g%+gi)\n",
i2, freq_re2[i2], freq_im2[i2], freq_err2[i2],
real(amp2[i2]), imag(amp2[i2]));
complex<double> ph = amp2[i2] / amp[i]
/ polar(1.0, 2*pi * (kx * (xsrc2 - xsrc) + ky * (ysrc2 - ysrc)));
master_printf(" (rel. amp. = %g%+gi, phase err = %g)\n",
real(ph), imag(ph), abs(ph - 1.0));
if (abs(ph - 1.0) > 0.2)
freq_re[i] = 0;
}
else
freq_re[i] = 0;
}
master_printf("freqs:, %d, %g, %g", kindex, kx, ky);
for (int i = 0; i < num_freqs - 1; ++i) // sort by freq_re
for (int j = i + 1; j < num_freqs; ++j)
if (freq_re[i] > freq_re[j]) {
double foo = freq_re[i]; freq_re[i] = freq_re[j]; freq_re[j] = foo;
}
for (int i = 0; i < num_freqs; ++i)
if (freq_re[i] != 0.0)
master_printf(", %g", freq_re[i]);
master_printf("\n");
if (hdf_output && which_k > 0) {
for (int i = 0; i < 2.0 / (freq * f.dt); ++i) {
f.output_hdf5(Ez, geometric_volume(vec(-1.0,-sqrt(3.0)),
vec(2.0, 2*sqrt(3.0))));
f.step();
}
master_printf("Output Ez.\n");
}
if (which_k > 0) break;
}
delete[] cpt2;
delete[] cpt;
return 0;
}
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
