Hi,
I've tried to convert a simple example of
http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial/Material_dispersioninto
a C++ code.
Attached is my simple code. I was wondering, although I use same variables
as the example why I get wrong variables
For K=[0.3, 0, 0] ctl file gives:
freqs:, 1, 0.3, 0, 0, 0.180392569033513, 1.22101799986358
While C++ gives:
Too many frequencies!
P.S: The only difference I can assume is that I use Ez component for Harminv
peak detection, maybe in ctl something else is used!
I really appreciate if somebody can check my file
Thanks
Joe
/*****************************MY CODE
(CPP)***********************************/
#include <meep.hpp>
using namespace meep;
typedef complex<double> dcomp;
//Structure
const double xWidth = 1.0;
const double yWidth = 1.0;
const double zWidth = 1.0;
const double resolution = 20; // pixels per distance
//Source
const double freq = 1.0;
const double dfreq = 0.9; //fractional
const component srcComp = Ez;
const vec srcPosition = vec (0.5, 0.5, 0.5);
//Sampling
const component sampleComp = Ez;
const int maxNumOfBands = 10;
//Run Time
const double timeAfterSource = 200.0;
//functions
double eps(const vec &p) {
return 2.25;
}
double sigma1(const vec &p) {
return 0.5;
}
double sigma2(const vec &p) {
return 2e-5;
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
/** Structure */
grid_volume v = vol3d(xWidth, yWidth, zWidth, resolution);
structure s(v, eps);
/**Output*/
const char *dirname = make_output_directory("~/Desktop/run");
s.set_output_directory(dirname);
s.add_polarizability (sigma1, E_stuff, 1.1, 1e-5); //(sigma, field_type
E_stuff, omega_n, gamma_n);
s.add_polarizability (sigma2, E_stuff, 0.5, 0.1);
fields f(&s);
//Harminv variables
int NumOfSamples = int(timeAfterSource / f.dt + 0.5) + 1;
dcomp *samplePoint = new dcomp[NumOfSamples];
/** Source*/
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));
f.add_point_source(srcComp, src, srcPosition);
//
vec kVector = vec(0.3, 0.0, 0.0);
f.use_bloch(kVector);
/** Run*/
while (f.time() < f.last_source_time()) {
f.step();
}
int i = 0;
while (f.time() < f.last_source_time() + timeAfterSource) {
samplePoint[i] = f.get_field(sampleComp, srcPosition);
i++;
f.step();
}
/** Harmonic Inversion*/
dcomp amp[maxNumOfBands];
double freq_re[maxNumOfBands], freq_im[maxNumOfBands],
freq_err[maxNumOfBands];
int num1_freqs = do_harminv(samplePoint, i, f.dt, freq_min / 4.0,
freq_max * 4.0, maxNumOfBands,
amp, freq_re, freq_im, freq_err);
/** Results*/
for (int i = 0; i < num1_freqs; ++i) {
master_printf("freq[%d] = %g%+gi +/- %g (amp %g%+gi) Q = %g\n",i,
freq_re[i], freq_im[i], freq_err[i], real(amp[i]), imag(amp[i]), -0.5 *
freq_re[i] / freq_im[i]);
}
delete[] samplePoint;
return 0;
}
#include <meep.hpp>
using namespace meep;
typedef complex<double> dcomp;
//Structure
const double xWidth = 1.0;
const double yWidth = 1.0;
const double zWidth = 1.0;
const double resolution = 20; // pixels per distance
//Source
const double freq = 1.0;
const double dfreq = 0.9; //fractional
const component srcComp = Ez;
const vec srcPosition = vec (0.5, 0.5, 0.5);
//Sampling
const component sampleComp = Ez;
const int maxNumOfBands = 10;
//Run Time
const double timeAfterSource = 200.0;
//functions
double eps(const vec &p) {
return 2.25;
}
double sigma1(const vec &p) {
return 0.5;
}
double sigma2(const vec &p) {
return 2e-5;
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
/** Structure */
grid_volume v = vol3d(xWidth, yWidth, zWidth, resolution);
structure s(v, eps);
/**Output*/
const char *dirname = make_output_directory("~/Desktop/run");
s.set_output_directory(dirname);
s.add_polarizability (sigma1, E_stuff, 1.1, 1e-5); //(sigma, field_type E_stuff, omega_n, gamma_n);
s.add_polarizability (sigma2, E_stuff, 0.5, 0.1);
fields f(&s);
//Harminv variables
int NumOfSamples = int(timeAfterSource / f.dt + 0.5) + 1;
dcomp *samplePoint = new dcomp[NumOfSamples];
/** Source*/
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));
f.add_point_source(srcComp, src, srcPosition);
//
vec kVector = vec(0.3, 0.0, 0.0);
f.use_bloch(kVector);
/** Run*/
while (f.time() < f.last_source_time()) {
f.step();
}
int i = 0;
while (f.time() < f.last_source_time() + timeAfterSource) {
samplePoint[i] = f.get_field(sampleComp, srcPosition);
i++;
f.step();
}
/** Harmonic Inversion*/
dcomp amp[maxNumOfBands];
double freq_re[maxNumOfBands], freq_im[maxNumOfBands], freq_err[maxNumOfBands];
int num1_freqs = do_harminv(samplePoint, i, f.dt, freq_min / 4.0, freq_max * 4.0, maxNumOfBands,
amp, freq_re, freq_im, freq_err);
/** Results*/
for (int i = 0; i < num1_freqs; ++i) {
master_printf("freq[%d] = %g%+gi +/- %g (amp %g%+gi) Q = %g\n",i, freq_re[i], freq_im[i], freq_err[i], real(amp[i]), imag(amp[i]), -0.5 * freq_re[i] / freq_im[i]);
}
delete[] samplePoint;
return 0;
}
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
