Mode frequency calculations for empty lattices.

Former-commit-id: ca14f8d765399c5bf92bb64cc7f1dc6722d53ebb

qpms/latticegens.c

@@ -585,3 +585,18 @@ const PGenClassInfo PGen_xyWeb = {
 };
 
 
+size_t PGen_xyWeb_sizecap(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
+    double minR, bool inc_minR, double maxR, bool inc_maxR) 
+{
+  l2d_reduceBasis(b1, b2, &b1, &b2);
+  LatticeFlags lf = l2d_detectRightAngles(b1, b2, rtol);
+  double layer_min_height = l2d_hexWebInCircleRadius(b1, b2);
+  long layer = ceil(minR / layer_min_height);
+  if(!inc_minR && (layer * layer_min_height)  <= minR)
+    ++layer; 
+  long last_layer = floor(maxR / layer_min_height);
+  if(!inc_maxR && (last_layer * layer_min_height) >= maxR)
+    --(last_layer);
+  // TODO less crude estimate (this one should be safe, however)
+  return ((lf & ORTHOGONAL_01) ? 4 : 6) * (last_layer - layer + 1);
+}

qpms/lattices.h

@@ -542,7 +542,7 @@ static inline PGenCart3ReturnData PGen_next_cart3(PGen *g) {
 	}
 }
 
-/// Ǧenerate a point in a 2D real space (cartesian coordinates).
+/// Generate a point in a 2D real space (cartesian coordinates).
 static inline PGenCart2ReturnData PGen_next_cart2(PGen *g) {
 	if (g->c->next_cart2) 
 		return g->c->next_cart2(g);
@@ -618,11 +618,14 @@ PGen PGen_1D_new_minMaxR(double period, ///< Distance between points.
 	PGen_1D_incrementDirection incdir ///< Order of generated points.
 );
 
-
 extern const PGenClassInfo PGen_xyWeb;
 PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, 
 		double minR, bool inc_minR, double maxR, bool inc_maxR);
 
+/// Returns a number larger or equal than the number of all the points generated by a PGen_xyWeb.
+size_t PGen_xyWeb_sizecap(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
+		double minR, bool inc_minR, double maxR, bool inc_maxR);
+
 
 /*
  * THE NICE PART (adaptation of lattices2d.py)
@@ -796,6 +799,58 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig, d
 		double maxr, bool maxr_inc);
 
 
+#if 0
+/// Gives the frequency of \a n-th empty lattice mode at a given wave vector \a k.
+double qpms_emptylattice2_mode_nth(
+                cart2_t b1_rec, ///< First reciprocal lattice base vector
+                cart2_t b2_rec, ///< Second reciprocal lattice base vector
+                double rtol, ///< Relative tolerance to detect right angles
+                cart2_t k, ///< The wave vector
+                double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index).
+                size_t N ///< Index of the mode (note that degenerate modes are counted multiple times).
+                );
+
+/// Gives the first `maxindex` frequencies of empty lattice modes at a given wave vector \a k.
+void qpms_emptylattice2_modes_maxindex(
+                double target_freqs[], ///< Target array of size maxindex.
+                cart2_t b1_rec, ///< First reciprocal lattice base vector
+                cart2_t b2_rec, ///< Second reciprocal lattice base vector
+                double rtol, ///< Relative tolerance to detect right angles
+                cart2_t k, ///< The wave vector
+                double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index).
+                size_t maxindex ///< Number of the frequencies generated.
+                );
+#endif
+
+/// Gives the frequencies of empty lattice modes at a given wave vector \a k up to \a maxfreq and one more.
+/**
+ * The frequencies are saved to a newly allocated array *target_freqs (to be deallocated
+ * using free() by the caller).
+ *
+ * \returns Number of found mode frequencies lower or equal than \a maxfreq plus one.
+ */
+size_t qpms_emptylattice2_modes_maxfreq(
+                double **target_freqs,
+                cart2_t b1_rec, ///< First reciprocal lattice base vector
+                cart2_t b2_rec, ///< Second reciprocal lattice base vector
+                double rtol, ///< Relative tolerance to detect right angles
+                cart2_t k, ///< The wave vector
+                double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index).
+                double maxfreq ///< The maximum frequency.
+                );
+
+/// Gives the frequencies of two empty lattice modes nearest to \a omega at a given wave vector \a k.
+void qpms_emptylattice2_modes_nearest(
+                double target[2], ///< Target array with lower ([0]) and upper ([1]) frequency.
+                cart2_t b1_rec, ///< First reciprocal lattice base vector
+                cart2_t b2_rec, ///< Second reciprocal lattice base vector
+                double rtol, ///< Relative tolerance to detect right angles
+                cart2_t k, ///< The wave vector
+                double wave_speed, ///< Wave speed in a given medium (i.e. vacuum speed / refractive index).
+                double omega ///< The frequency around which the frequencies are searched.
+                );
+
+
 
 /*
  * THE UGLY PART

qpms/lattices2d.c

@@ -866,4 +866,83 @@ double l2d_unitcell_area(cart2_t b1, cart2_t b2) {
 }
   
 
+#if 0
+double qpms_emptylattice2_mode_nth(
+                cart2_t b1_rec,
+                cart2_t b2_rec,
+                double rtol,
+                cart2_t k,
+                double c,
+                size_t N 
+                );
+
+void qpms_emptylattice2_modes_maxindex(
+                double target_freqs[], 
+                cart2_t b1_rec,
+                cart2_t b2_rec,
+                double rtol,
+                cart2_t k,
+                double c,
+                size_t maxindex 
+                );
+#endif
+
+static int dblcmp(const void *p1, const void *p2) {
+  const double *x1 = (double *) p1, *x2 = (double *) p2;
+  double dif = *x1 - *x2;
+  if(dif > 0) return 1;
+  else if (dif < 0) return -1;
+  else return 0;
+}
+
+size_t qpms_emptylattice2_modes_maxfreq(double **target_freqs, 
+  cart2_t b1, cart2_t b2, double rtol, cart2_t k, 
+  double c, double maxfreq)
+{
+  QPMS_UNTESTED;
+  l2d_reduceBasis(b1, b2, &b1, &b2);
+  double maxk_safe = cart2norm(k) + maxfreq / c + cart2norm(b1) + cart2norm(b2); // This is an overkill
+  size_t capacity = PGen_xyWeb_sizecap(b1, b2, rtol, k, 0, true, maxk_safe, true);
+  cart2_t *Kpoints;
+  QPMS_CRASHING_MALLOC(Kpoints, sizeof(cart2_t) * capacity);
+  PGen Kgen = PGen_xyWeb_new(b1, b2, rtol, k, 0, true, maxk_safe, true);
+
+  size_t generated;
+  PGenReturnDataBulk rd;
+  while((rd = PGen_fetch_cart2(&Kgen, capacity - generated, Kpoints + generated)).flags & PGEN_NOTDONE)
+    generated += rd.generated;
+
+  double *thefreqs;
+  QPMS_CRASHING_MALLOC(thefreqs, generated * sizeof(double));
+  for(size_t i = 0; i < generated; ++i) thefreqs[i] = cart2norm(Kpoints[i]) * c;
+
+  free(Kpoints);
+
+  qsort(thefreqs, generated, sizeof(double), dblcmp);
+  size_t count;
+  bool hitmax = false;
+  for(count = 0; count < generated; ++count)
+    if(thefreqs[count] > maxfreq) {
+      if(hitmax)
+        break;
+      else 
+        hitmax = true;
+    }
+
+  *target_freqs = thefreqs;
+  return count;
+}
+
+void qpms_emptylattice2_modes_nearest(double target[2],
+    cart2_t b1, cart2_t b2, double rtol,
+    cart2_t k, double c, double omega) 
+{
+  QPMS_UNTESTED;
+  double *freqlist;
+  size_t n = qpms_emptylattice2_modes_maxfreq(&freqlist,
+      b1, b2, rtol, k, c, omega);
+  target[0] = freqlist[n-2];
+  target[1] = freqlist[n-1];
+  free(freqlist);
+}