WIP generate 3D lattices

Former-commit-id: 16c6f7779db68a3169dbc77a6adb9ffdd6ad374a

qpms/latticegens.c

@@ -585,3 +585,170 @@ const PGenClassInfo PGen_xyWeb = {
 };
 
 
+//==== PGen_xyzWeb ====
+// 3D lattice generator in the "spiderweb" style, generated in the "perimetre" order,
+// not strictly ordered (or limited) by distance from origin.
+// The minR and maxR here refer to the TODO WWHAT
+
+extern const PGenClassInfo PGen_xyzWeb; // forward declaration needed by constructor (may be placed in header file instead)
+
+// Internal state structure
+typedef struct PGen_xyzWeb_StateData {
+  long i[6]; // real size given by nsteps
+  unsigned short phase; // 0 to 5
+  long layer;
+  long last_layer; // generation stops when layer > last_layer
+  double layer_min_height; // this * layer is what minR and maxR are compared to
+  double minR, maxR;
+  bool inc_minR, inc_maxR;
+  cart3_t b[3]; // lattice vectors (TODO needed?)
+  cart3_t steps[6]; // generated by l3d_shortestBase356
+  int nsteps; // real size of steps (3, 5, or 6), depending on lattice type
+  cart3_t offset; // offset of the zeroth lattice point from origin (TODO will be normalised to the WS cell)
+  // TODO type rectangular vs. triangular
+  LatticeFlags lf;
+} PGen_xyzWeb_StateData;
+
+// Constructor
+PGen PGen_xyzWeb_new(const cart3_t b[3], double rtol, cart2_t offset, double minR, bool inc_minR, double maxR, bool inc_maxR) {
+  PGen_xyzWeb_StateData *s =  calloc(1, sizeof(PGen_xyzWeb_StateData));
+  s->minR = minR; s->maxR = maxR;
+  s->inc_minR = inc_minR;
+  s->inc_maxR = inc_maxR;
+  l3d_reduceBasis(b, &(s->b));
+  s->nsteps = l3d_shortestBase356(s->b, s->steps);
+  s->offset = offset; // TODO shorten into the WS cell ?
+  s->lf = l3d_detectRightAnles(s->b, rtol);
+  s->layer_min_height = l3d_inSphereRadius(s->b);
+  s->layer = ceil(s->minR/s->layer_min_height);
+  if(!inc_minR && (s->layer * s->layer_min_height)  <= minR)
+    ++(s->layer); 
+  s->i = s->layer; s->j = 0; s->phase = 0; // init indices
+  s->last_layer = floor(s->maxR/s->layer_min_height);
+  if(!inc_maxR && (s->last_layer * s->layer_min_height) >= maxR)
+    --(s->last_layer);
+  PGen g = {&PGen_xyzWeb, (void *) s};
+  return g;
+}
+
+// Destructor
+void PGen_xyzWeb_destructor(PGen *g) {
+  free(g->stateData);
+  g->stateData = NULL;
+}
+
+// Extractor (3D cartesian, native)
+PGenCart3ReturnData PGen_xyzWeb_next_cart3(PGen *g) {
+  if (g->stateData == NULL) // already destroyed
+    return PGenCart2DoneVal;
+  else {
+    PGen_xyzWeb_StateData * const s = (PGen_xyzWeb_StateData *) g->stateData;
+    assert(s->layer >= 0);
+    if (s->layer <= s->last_layer) {
+      const cart2_t thePoint = cart2_add(s->offset,
+          cart2_add(cart2_scale(s->i, s->b1), cart2_scale(s->j, s->b2)));
+      if(s->layer == 0) { // origin is unique, proceed with next layer
+        ++s->layer;
+        s->phase = 0;
+        s->i = s->layer;
+        s->j = 0;
+      }
+      else if(s->lt & (SQUARE | RECTANGULAR)) {
+        // rectangular or square lattice, four perimeters
+        switch(s->phase) {
+          case 0: // initial i = l, j = 0
+            --s->i;
+            ++s->j;
+            if(s->i <= 0) ++s->phase;
+            break;
+          case 1: // initial i = 0, j = l
+            --s->i;
+            --s->j;
+            if(s->j <= 0) ++s->phase;
+            break;
+          case 2: // initial i = -l, j = 0
+            ++s->i;
+            --s->j;
+            if(s->i >= 0) ++s->phase;
+            break;
+          case 3: // initial i = 0, j = -l
+            ++s->i;
+            ++s->j;
+            if(s->j >= 0) ++s->phase;
+            break;
+          default:
+            abort();
+        }
+        if(s->phase == 4) { // phase overflow, start new layer
+          ++s->layer;
+          s->phase = 0;
+          s->i = s->layer;
+          s->j = 0;
+        }
+      } else { // non-rectangular lattice, six perimeters
+        switch(s->phase) {
+          case 0:
+            --s->i;
+            ++s->j;
+            if(s->i <= 0) ++s->phase;
+            break;
+          case 1:
+            --s->i;
+            if(s->i + s->j <= 0) ++s->phase;
+            break;
+          case 2:
+            --s->j;
+            if(s->j <= 0) ++s->phase;
+            break;
+          case 3:
+            ++s->i;
+            --s->j;
+            if(s->i >= 0) ++s->phase;
+            break;
+          case 4:
+            ++s->i;
+            if(s->i + s->j >= 0) ++s->phase;
+            break;
+          case 5:
+            ++s->j;
+            if(s->j >= 0) ++s->phase;
+            break;
+          default:
+            abort();
+        }
+        if(s->phase == 6) { // phase overflow, start next layer
+          ++s->layer;
+          s->phase = 0;
+          s->i = s->layer;
+          s->j = 0;
+        }
+      }
+      PGenCart2ReturnData retval = {(PGEN_NOTDONE | PGEN_AT_XY | PGEN_COORDS_CART2), thePoint};
+      return retval;
+    } else {
+      PGen_destroy(g);
+      return PGenCart2DoneVal;
+    }
+  }
+}
+
+// Class metadata structure; TODO maybe this can rather be done by macro.
+const PGenClassInfo PGen_xyzWeb = {
+  "PGen_xyzWeb",
+  3,
+  PGEN_COORDS_CART3,
+  NULL,//PGen_xyzWeb_next,
+  NULL,//PGen_xyzWeb_next_z,
+  NULL,//PGen_xyzWeb_next_pol,
+  PGen_next_sph_from_cart3, //NULL,//PGen_xyzWeb_next_sph,
+  NULL,//PGen_xyzWeb_next_cart2,
+  PGen_xyzWeb_next_cart3, // native
+  NULL,//PGen_xyzWeb_fetch,
+  NULL,//PGen_xyzWeb_fetch_z,
+  NULL,//PGen_xyzWeb_fetch_pol,
+  NULL,//PGen_xyzWeb_fetch_sph,
+  NULL,//PGen_xyzWeb_fetch_cart2,
+  NULL,//PGen_xyzWeb_fetch_cart3,
+  PGen_xyzWeb_destructor
+};
+

qpms/lattices.h

@@ -620,7 +620,11 @@ PGen PGen_1D_new_minMaxR(double period, ///< Distance between points.
 
 
 extern const PGenClassInfo PGen_xyWeb;
-PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, 
+PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double orthogonality_rtol, cart2_t offset, 
+		double minR, bool inc_minR, double maxR, bool inc_maxR);
+
+extern const PGenClassInfo PGen_xyzWeb;
+PGen PGen_xyzWeb_new(const cart3_t basis[3], double orthogonality_rtol, cart3_t offset,
 		double minR, bool inc_minR, double maxR, bool inc_maxR);
 
 
@@ -691,6 +695,8 @@ int l2d_shortestBase46(cart2_t i1, cart2_t i2,  cart2_t *o1, cart2_t *o2, cart2_
 // variant
 int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
 
+int l3d_shortestBase356(const cart3_t basis[3], cart3_t shortest356[6], double rtol);
+
 // Determines whether angle between inputs is obtuse
 bool l2d_is_obtuse_r(cart2_t i1, cart2_t i2, double rtol);
 bool l2d_is_obtuse(cart2_t i1, cart2_t i2);
@@ -729,6 +735,8 @@ double l2d_unitcell_area(cart2_t b1, cart2_t b2);
 // returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple
 double l2d_hexWebInCircleRadius(cart2_t b1, cart2_t b2);
 
+double l3d_inSphereRadius(const cart3_t b[3]);
+
 /*
  * THE MORE OR LESS OK PART
  * ========================