Dear MEEP users, let me share a little code snippet with you. Sometimes it happens that one wants to run a FDTD simulation with a broadband source, which however avoids some frequencies. A single gaussian source is not useful in such case, as its spectrum has always broad wings. Composing the source from many narrow-band gaussians is not practical.
For such a case I prepared a band_src_time() class, which is used similarly as the gaussian source except for that it transmits only in a spectral window given by parameters f, fwidth. Technically the source is a sinc(t) function multiplied by Blackman-Nutall window. Such a waveform has proven to provide a nice rectangle with highly suppressed side lobes and therefore spectral leakage under -25 dB with default settings of a short source duration (much better values available with longer source time). To apply, put the following snippet in the sources.cpp in MEEP 1.2 sources: ---------------------------------------------------------------------------------------------------------- /// Band source band_src_time::band_src_time(double f, double fwidth, double s) { freq = f; width = 1.0 / fwidth; peak_time = width * s; cutoff = width * s * 2; master_printf("Initializing band source for time_domain (f fwidth s)"); // this is to make last_source_time as small as possible while (exp(-cutoff*cutoff / (2*width*width)) < 1e-100) cutoff *= 0.9; cutoff = float(cutoff); // don't make cutoff sensitive to roundoff error } complex<double> band_src_time::dipole(double time) const { // e_i2pif = exp(-1j*2*pi*freq*tt) // wnd_Hann = cos(-(tt/(cutoff*2)*pi))**2 // wnd_Hamming= (0.53836+0.46164*cos(tt/(cutoff*2)*pi*2)) // wnd_Blackman= (0.42659 + 0.49656*cos(tt/(cutoff*2)*pi*2) + 0.076849*cos(tt/(cutoff*2)*pi*4)) // wnd_Gauss = exp(-((tt)/(cutoff*2)*3)**2) // Gaussian window, char. width = 1/3 // wnd_ContRect =rect(tt, 0, (trunc((cutoff*2)/8/width))*8*width) double tt = time - peak_time; double sinc = sin(tt*2*pi / width/2)/(tt*2*pi * width); double wnd_BlackmanN= (0.3635819 + 0.4891775*cos(tt/(cutoff*2)*pi*2) + 0.1365995*cos(tt/(cutoff*2)*pi*4)+ 0.0106411*cos(tt/(cutoff*2)*pi*6)); //if (float(fabs(tt)) > cutoff) //return 0.0; // correction factor so that current amplitude (= d(dipole)/dt) is // ~ 1 near the peak of the band. complex<double> amp = 1.0 / complex<double>(0,-2*pi*freq); return sinc * wnd_BlackmanN * polar(1.0, -2*pi*freq*tt) * amp; //if (abs(tt)>1e-9) //return exp(-tt*tt / (2*width*width)) * sin(tt)/tt * polar(1.0, -2*pi*freq*tt) * amp; //else //return 1; } bool band_src_time::is_equal(const src_time &t) const { const band_src_time *tp = dynamic_cast<const band_src_time*>(&t); if (tp) return(tp->freq == freq && tp->width == width && tp->peak_time == peak_time && tp->cutoff == cutoff); else return 0; } ---------------------------------------------------------------------------------------------------------- And add the following into meep.hpp: ---------------------------------------------------------------------------------------------------------- // band source with given minimum and maximum frequency and duration class band_src_time : public src_time { public: band_src_time(double f, double fwidth, double s = 5.0); virtual ~band_src_time() {} virtual complex<double> dipole(double time) const; virtual double last_time() const { return float(peak_time + cutoff); }; virtual src_time *clone() const { return new band_src_time(*this); } virtual bool is_equal(const src_time &t) const; virtual complex<double> frequency() const { return freq; } virtual void set_frequency(complex<double> f) { freq = real(f); } private: double freq, width, peak_time, cutoff; }; ---------------------------------------------------------------------------------------------------------- _______________________________________________ meep-discuss mailing list meep-discuss@ab-initio.mit.edu http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss