Latticegen xyWeb, slightly tested as well.

Former-commit-id: 1a51ed4080f5a17b797b27fb2425052beacfb687

qpms/latticegens.c

@@ -432,8 +432,12 @@ PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) {
     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
+      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) {

qpms/lattices2d.c

@@ -733,6 +733,39 @@ int l2d_shortestBase46(const cart2_t i1, const cart2_t i2,  cart2_t *o1, cart2_t
   }
 }
 
+/*
+ * Given two basis vectors, returns 2D Bravais lattice type.
+ */
+LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol)
+{
+  l2d_reduceBasis(b1, b2, &b1, &b2);
+  cart2_t b3 = cart2_substract(b2, b1);
+  double b1s = cart2_normsq(b1), b2s = cart2_normsq(b2), b3s = cart2_normsq(b3);
+  double eps = rtol * (b2s + b1s); //FIXME what should eps be?
+  // avoid obtuse angle between b1 and b2. TODO this should be yet tested
+  // TODO use is_obtuse here?
+  if (b3s - b2s - b1s > eps) {
+    b3 = b2;
+    b2 = cart2_add(b2, b1);
+    // N.B. now the assumption |b3| >= |b2| is no longer valid
+    // b3 = cart2_substract(b2, b1)
+    b2s = cart2_normsq(b2); 
+    b3s = cart2_normsq(b3);
+  }
+  if (fabs(b2s-b1s) < eps || fabs(b2s - b3s) < eps) { // isoscele
+    if (fabs(b3s-b1s) <  eps)
+      return EQUILATERAL_TRIANGULAR;
+    else if (fabs(b3s - 2*b1s))
+      return SQUARE;
+    else
+      return RHOMBIC;
+  } else if (fabs(b3s-b2s-b1s) < eps)
+    return RECTANGULAR;
+  else
+    return OBLIQUE;
+}
+
+
 # if 0
 // variant
 int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
@@ -740,10 +773,6 @@ int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
 // Determines whether angle between inputs is obtuse
 bool l2d_is_obtuse_r(cart2_t i1, cart2_t i2, double rtol);
 
-/*
- * Given two basis vectors, returns 2D Bravais lattice type.
- */
-LatticeType l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol);
 
 // Other functions in lattices2d.py: TODO?
 // range2D()

tests/test_latticegenxyweb.c

@@ -0,0 +1,54 @@
+//c99 -o test_latticegenxyweb -ggdb -Wall -I ../ test_latticegenxyweb.c ../qpms/latticegens.c -lm
+#define QPMS_VECTORS_NICE_TRANSFORMATIONS
+#include <qpms/lattices.h>
+#include <stdio.h>
+
+void fprint_PGenCart2ReturnData(FILE *f, PGenCart2ReturnData d) {
+  cart2_t c = d.point_cart2;
+  fprintf(f, "%g\t%g\tflags: %#xd\n",
+      c.x, c.y, d.flags);
+}
+
+void dump_PGenCart2(char *filename, PGen *g) {
+  FILE *out = fopen(filename, "w");
+  PGenCart2ReturnData d;
+  while ((d = PGen_next_cart2(g)).flags & PGEN_NOTDONE) {
+    fprint_PGenCart2ReturnData(out, d);
+  } 
+}
+
+#if 0
+void print_PGenSphReturnData(PGenSphReturnData d) {
+  sph_t s = d.point_sph;
+  cart3_t c = sph2cart(s);
+  printf("sph: (%g, %g π, %g π), cart: (%g, %g, %g), flags: %#xd\n",
+      s.r, s.theta / M_PI, s.phi / M_PI, c.x, c.y, c.z, d.flags);
+}
+
+void dump_PGenSph(PGen *g) {
+  PGenSphReturnData d;
+  do {
+    d = PGen_next_sph(g);
+    print_PGenSphReturnData(d);
+  } while (d.flags & PGEN_NOTDONE);
+}
+#endif
+
+#define DO_AND_PRINT(label, s) printf(#label ":\n" #s "\n"); s ;
+
+int main(int argc, char **argv) {
+  PGen g;
+  cart2_t b1 = {1.1,0}, b2 = {0,-1.3}, offset = {0.1, -0.1};
+  g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 0, true, 5, false);
+  dump_PGenCart2("rect1.xydump", &g);
+  g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 5, true, 8, false);
+  dump_PGenCart2("rect2.xydump", &g);
+  b1.x = 1.342, b1.y = 4.3121; b2.x = -1.83, b2.y = 1.4;
+  g = PGen_xyWeb_new(b1, b2, 1e-13, offset, 0, true, 8, false);
+  dump_PGenCart2("oblique1.xydump", &g);
+
+  return 0;
+
+
+}
+