Unfinished generation of equilateral triangular lattices (dudom).

Former-commit-id: 4cc4da2e9755c279aba0d14c296f76de2cbebd04

qpms/lattices.h

@@ -0,0 +1,88 @@
+#ifndef LATTICES_H
+#define LATTICES_H
+#include <math.h>
+#include <stdbool.h>
+
+#define M_SQRT3 1.7320508075688772935274463415058724
+
+// This might be reduced to x, y only; not sure yet
+typedef struct {
+	double key; // distance key in a given lattice
+	double x, y, r, phi;
+} point2d;
+
+
+static inline point2d point2d_fromxy(const double x, const double y) {
+	point2d p;
+	p.x = x;
+	p.y = y;
+	p.r = sqrt(x*x+y*y);
+	p.phi = atan2(y, x);
+	return p;
+}
+
+/* 
+ * General set of points ordered by the r-coordinate.
+ * Typically, this will include all lattice inside a certain circle.
+ * This structure is internally used by the "lattice generators" below.
+ * It does not have its memory management of its own, as it is handled
+ * by the "generators". For everything except the generators,
+ * this structure shall be read-only.
+ */
+typedef struct {
+	size_t nrs; // number of different radii
+	double *rs; // the radii; of length nrs (largest contained radius == rs[nrs-1])
+	point2d **points_at_r; // of length nrs+1
+	/* // redundand (therefore removed) members
+         * point2d *points; // redundant as it is the same as points_at_r[0]
+	 * size_t npoints;  // redundant as it is the same as points_at_r[nrs]-points_at_r[0]
+	 */ 
+} points2d_rordered_t;
+
+
+
+/* 
+ * EQUILATERAL TRIANGULAR LATTICE
+ */
+
+typedef enum {
+	TRIANGULAR_VERTICAL, // there is a lattice base vector parallel to the y-axis
+  	TRIANGULAR_HORIZONTAL // there is a lattice base vector parallel to the x-axis
+} TriangularLatticeOrientation;
+
+// implementation data structures; not needed in the header file
+typedef struct triangular_lattice_gen_privstuff_t triangular_lattice_gen_privstuff_t;
+
+typedef struct {
+	// public:
+	points2d_rordered_t ps;
+	TriangularLatticeOrientation orientation;
+	double a; // lattice vector length
+
+	// not sure if needed:
+	bool includes_origin;
+
+	// private:
+	triangular_lattice_gen_privstuff_t *priv;
+
+} triangular_lattice_gen_t;
+
+triangular_lattice_gen_t *triangular_lattice_gen_init(double a, TriangularLatticeOrientation ori, bool include_origin); 
+const points2d_reordered_t * triangular_lattice_gen_getpoints(const triangular lattice_generator_t *g);
+int triangular_lattice_gen_extend_to_r(triangular_lattice_generator_t *g, double r);
+int triangular_lattice_gen_extend_to_steps(triangular_lattice_generator_t *g, int maxsteps);
+void triangular_lattice_gen_free(triangular_lattice_generator_t *g);
+
+#if 0
+
+/*
+ * HONEYCOMB LATTICE
+ */
+
+typedef struct {
+
+} honeycomb_lattice_generator_t;
+
+#endif
+
+#endif // LATTICES_H

qpms/lattices2d.c

@@ -0,0 +1,113 @@
+#include "lattices.h"
+#include <assert.h>
+
+typedef struct {
+  int i, j;
+} intcoord2d;
+
+static inline int sqi(int x) { return x*x; }
+
+
+/*
+ * EQUILATERAL TRIANGULAR LATTICE
+ */
+
+
+/*
+ * N. B. the possible radii (distances from origin) of the lattice points can be described as
+ *
+ *     r**2 / a**2 == i**2 + j**2 + i*j ,
+ *
+ * where i, j are integer indices describing steps along two basis vectors (which have
+ * 60 degree angle between them).
+ *
+ * The plane can be divided into six sextants, characterized as:
+ *
+ * 0) i >= 0 && j >= 0,
+ *    [a] i > 0,
+ *    [b] j > 0,
+ * 1) i <= 0 && {j >= 0} && i + j >= 0,
+ *    [a] i + j > 0,
+ *    [b] i < 0,
+ * 2) {i <= 0} && j >= 0 && i + j <= 0,
+ *    [a] j > 0,
+ *    [b] i + j < 0,
+ * 3) i <= 0 && j <= 0,
+ *    [a] i < 0,
+ *    [b] j < 0,
+ * 4) i >= 0 && {j <= 0} && i + j <= 0,
+ *    [a] i + j < 0,
+ *    [b] i > 0,
+ * 5) {i >= 0} && j <= 0 && i + j >= 0,
+ *    [a] j < 0,
+ *    [b] i + j > 0.
+ *
+ * The [a], [b] are two variants that uniquely assign the points at the sextant boundaries.
+ * The {conditions} in braces are actually redundant.
+ *
+ * In each sextant, the "minimum steps from the origin" value is calculated as:
+ * 0) i + j,
+ * 1) j
+ * 2) -i
+ * 3) -i - j,
+ * 4) -j,
+ * 5) i.
+ *
+ * The "spider web" generation for s steps from the origin (s-th layer) goes as following (variant [a]):
+ * 0) for (i = s, j = 0; i > 0; --i, ++j)
+ * 1) for (i = 0, j = s; i + j > 0; --i)
+ * 2) for (i = -s, j = s; j > 0; --j)
+ * 3) for (i = -s, j = 0; i < 0; ++i, --j)
+ * 4) for (i = 0, j = -s; i + j < 0; ++i)
+ * 5) for (i = s, j = -s; j < 0; ++j)
+ *
+ *
+ * Length of the s-th layer is 6*s for s >= 1. Size (number of lattice points) of the whole s-layer "spider web"
+ * is therefore 3*s*(s+1), excluding origin.
+ * The real area inside the web is (a*s)**2 * 3 * sqrt(3) / 2.
+ * Area of a unit cell is a**2 * sqrt(3)/2.
+ * Inside the web, but excluding the circumscribed circle, there is no more
+ * than 3/4.*s*(s+1) lattice cells (FIXME pretty stupid but safe estimate).
+ *
+ * s-th layer circumscribes a circle of radius a * s * sqrt(3)/2.
+ *
+ */
+
+static inline int trilat_r2_ij(const int i, const int j) {
+  return sqi(i) + sqi(j) + i*j;
+}
+
+// Classify points into sextants (variant [a])
+static int trilat_sextant_ij_a(const int i, const int j) {
+  const int w = i + j;
+  if (i >  0 && j >= 0) return  0;
+  if (i <= 0 && w >  0) return  1;
+  if (w <= 0 && j >  0) return  2;
+  if (i <  0 && j <= 0) return  3;
+  if (i >= 0 && w <  0) return  4;
+  if (w >= 0 && j <  0) return  5;
+  if (i == 0 && j == 0) return -1; // origin
+  assert(0); // other options should be impossible
+}
+
+
+
+typedef struct {
+  TODO;
+} triangular_lattice_gen_t_privstuff_t;
+
+triangular_lattice_gen_t * triangular_lattice_gen_init(double a, TriangularLatticeOrientation ori, bool include_origin)
+{
+  triangular_lattice_gen_t *g = malloc(sizeof(triangular_latice_gen_t));
+  g->orientation = ori;
+  g->includes_origin = include_origin;
+  g->ps.nrs = 0;
+  g->ps.rs = NULL;
+  g->ps.points_at_r = NULL;
+  g->priv = malloc(sizeof(triangular_lattice_gen_privstuff_t));
+  TODO;
+
+  return g;
+}
+
+