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