diff --git a/src/meep.hpp b/src/meep.hpp index 416c63ab6..c1a66c0be 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -737,6 +737,24 @@ class gaussian_src_time : public src_time { double freq, width, peak_time, cutoff; }; + +// 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 user_cutoff); + virtual ~band_src_time() {} + + virtual std::complex 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 std::complex frequency() const { return freq; } + virtual void set_frequency(std::complex f) { freq = real(f); } + + private: + double freq, width, peak_time, cutoff; + //complex test; +}; // Continuous (CW) source with (optional) slow turn-on and/or turn-off. class continuous_src_time : public src_time { public: diff --git a/src/sources.cpp b/src/sources.cpp index 311f6805c..bec6bde0e 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -23,6 +23,13 @@ #include "meep.hpp" #include "meep_internals.hpp" +// For sine-integral function used in the band-source +#define HAVE_LIBGSL_EXPINT + +#ifdef HAVE_LIBGSL_EXPINT +# include +#endif + using namespace std; namespace meep { @@ -109,6 +116,70 @@ bool gaussian_src_time::is_equal(const src_time &t) const else return 0; } +/// Band source +band_src_time::band_src_time(double f, double fwidth, double user_cutoff) +{ + freq = f; + width = 1.0 / fwidth; + peak_time = user_cutoff/2; + cutoff = user_cutoff; +#ifdef HAVE_LIBGSL_EXPINT + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tExperimental: using GSL si() function for flat-top spectrum\n"); +#else + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tWarning: not compiled with GSL, the source spectrum will not be flat-top\n", peak_time, cutoff); +#endif + cutoff = float(cutoff); // don't make cutoff sensitive to roundoff error +} + +complex band_src_time::dipole(double time) const +{ + double tt = time - peak_time; + + // The function that introduces a rectangular band in the spectrum (centered around zero frequency) +#ifdef HAVE_LIBGSL_EXPINT + // The emitted field is the derivative of the dipole, so if possible, we use the sine integral: + // Si(x) = \int_0^x dt sin(t)/t + complex func = gsl_sf_Si(tt*2*pi*(freq-.5/width)) - gsl_sf_Si(tt*2*pi*(freq+.5/width)); +#else + // If GSL not available, a reasonable approximation is the sinc(t) function (but has not flat top) + complex func = sin(tt*2*pi / width/2)/(tt*2*pi * width) * polar(1.0, -2*pi*freq*tt); +#endif + + // The envelope that suppresses ringing and side lobes of the rectangle + double wnd_BlackmanN= (0.3635819 + 0.4891775*cos(tt/(cutoff)*pi*2) + + 0.1365995*cos(tt/(cutoff)*pi*4)+ 0.0106411*cos(tt/(cutoff)*pi*6)); + //double wnd_Hann = (.5 + .5*cos(tt/(cutoff)*pi*2)); + //double wnd_Hamming = (0.53836+0.46164*cos(tt/(cutoff)*pi*2)) ; + //double wnd_Blackman = (0.42659 + 0.49656*cos(tt/(cutoff)*pi*2) + 0.076849*cos(tt/(cutoff)*pi*4)) ; + //double wnd_Gauss = exp(-((tt)/(cutoff)*3)**2) // Gaussian window, char. width = 1/3; + //double wnd_ContRect = // should crop sinc when passing zero: rect(tt, 0, (trunc((cutoff)/8/width))*8*width) ; + + // correction factor so that current amplitude (= d(dipole)/dt) is + // ~ 1 near the peak of the band. + complex amp = 1.0 / complex(0,-2*pi*freq); // TODO + + // The needed time-domain source amplitude is the sinc function constrained by window and shifted + // in frequency by complex exponential + if (abs(tt) < cutoff/2) + //return sinc; + return func * wnd_BlackmanN; + //return func * wnd_BlackmanN * amp; + else + return 0; + +} + +bool band_src_time::is_equal(const src_time &t) const +{ + const band_src_time *tp = dynamic_cast(&t); + if (tp) + return(tp->freq == freq && tp->width == width && + tp->peak_time == peak_time && tp->cutoff == cutoff); + else + return 0; +} complex continuous_src_time::dipole(double time) const { diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index c3c1a9184..7aac77b7f 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -111,6 +111,21 @@ void *lorentzian_susceptibility::new_internal_data( return (void*) d; } + +/* Return true if the discretized Lorentzian ODE is intrinsically unstable, + i.e. if it corresponds to a filter with a pole z outside the unit circle. + Note that the pole satisfies the quadratic equation: + (z + 1/z - 2)/dt^2 + g*(z - 1/z)/(2*dt) + w^2 = 0 + where w = 2*pi*omega_0 and g = 2*pi*gamma. It is just a little + algebra from this to get the condition for a root with |z| > 1. */ +static bool lorentzian_unstable(double omega_0, double gamma, double dt) { + double w = 2*pi*omega_0, g = 2*pi*gamma; + double g2 = g*dt/2, w2 = (w*dt)*(w*dt); + double b = (1 - w2/2) / (1 + g2), c = (1 - g2) / (1 + g2); + return b*b > c && 2*b*b - c + 2*fabs(b)*sqrt(b*b - c) > 1; +} + + void lorentzian_susceptibility::init_internal_data( realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *data) const { @@ -128,6 +143,11 @@ void lorentzian_susceptibility::init_internal_data( P += 2*ntot; P_prev += 2*ntot; } + + // TODO - ensure this is printed only once for each oscillator, not for every chunk + if (!no_omega_0_denominator && gamma >= 0 + && lorentzian_unstable(omega_0, gamma, dt)) + master_printf("Warning: Lorentzian may be unstable according to the Von Neumann stability analysis (omega_0=%g, gamma=%g, dt=%g). Proceed with caution.\n", omega_0, gamma, dt); } void *lorentzian_susceptibility::copy_internal_data(void *data) const { @@ -147,18 +167,6 @@ void *lorentzian_susceptibility::copy_internal_data(void *data) const { return (void*) dnew; } -/* Return true if the discretized Lorentzian ODE is intrinsically unstable, - i.e. if it corresponds to a filter with a pole z outside the unit circle. - Note that the pole satisfies the quadratic equation: - (z + 1/z - 2)/dt^2 + g*(z - 1/z)/(2*dt) + w^2 = 0 - where w = 2*pi*omega_0 and g = 2*pi*gamma. It is just a little - algebra from this to get the condition for a root with |z| > 1. */ -static bool lorentzian_unstable(double omega_0, double gamma, double dt) { - double w = 2*pi*omega_0, g = 2*pi*gamma; - double g2 = g*dt/2, w2 = (w*dt)*(w*dt); - double b = (1 - w2/2) / (1 + g2), c = (1 - g2) / (1 + g2); - return b*b > c && 2*b*b - c + 2*fabs(b)*sqrt(b*b - c) > 1; -} #define SWAP(t,a,b) { t SWAP_temp = a; a = b; b = SWAP_temp; } @@ -177,10 +185,6 @@ void lorentzian_susceptibility::update_P const double omega0dtsqr_denom = no_omega_0_denominator ? 0 : omega0dtsqr; (void) W_prev; // unused; - if (!no_omega_0_denominator && gamma >= 0 - && lorentzian_unstable(omega_0, gamma, dt)) - abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); - FOR_COMPONENTS(c) DOCMP2 if (d->P[c][cmp]) { const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)]; if (w && s) {