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