Band source

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<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 std::complex<double> frequency() const { return freq; }
+  virtual void set_frequency(std::complex<double> f) { freq = real(f); }
+
+ private:
+  double freq, width, peak_time, cutoff;
+  //complex<double> test;
+};
 // Continuous (CW) source with (optional) slow turn-on and/or turn-off.
 class continuous_src_time : public src_time {
  public:

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 <gsl/gsl_sf_expint.h>
+#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<double> 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<double> 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<double> 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<double> amp = 1.0 / complex<double>(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<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;
+}
 
 complex<double> continuous_src_time::dipole(double time) const
 {

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) {