A stupid shift function for points2d_rordered_t

Former-commit-id: bbba5408081965d8550d78fed5b7569c5cff95b8

qpms/lattices.h

@@ -105,7 +105,7 @@ int l2d_cellCornersWS_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
 
 // Reciprocal bases
 void l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
-void l2D_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
+void l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
 
 
 /*
@@ -133,10 +133,26 @@ typedef struct {
 	 */ 
 } points2d_rordered_t;
 
-// returns a copy but scaled by a factor
-points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig, double factor);
 
-void points2d_rordered_free(points2d_rordered_t *); // use only for result of points2d_rordered_scale
+// sorts arbitrary points and creates points2d_rordered_t
+points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base,
+	       	size_t nmemb, double rtol, double atol);
+
+// returns a copy but shifted by a constant (actually in a stupid way, but whatever)
+points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig,
+		point2d shift, double rtol, double atol);
+
+// returns a copy but scaled by a factor
+points2d_rordered_t *points2d_rordered_scale(const points2d_rordered_t *orig,
+       	double factor);
+
+
+/* The destructor: use only for results of
+ *  - points2D_rordered_frompoints,
+ *  - points2d_rordered_shift,
+ *  - points2d_rordered_scale.
+ */
+void points2d_rordered_free(points2d_rordered_t *);
 
 static inline point2d points2d_rordered_get_point(const points2d_rordered_t *ps, int r_order, int i) {
 	assert(i >= 0);

qpms/lattices2d.c

@@ -1,6 +1,7 @@
 #include "lattices.h"
 #include <assert.h>
 #include <stdlib.h>
+#include <string.h>
 
 typedef struct {
   int i, j;
@@ -8,6 +9,7 @@ typedef struct {
 
 static inline int sqi(int x) { return x*x; }
 
+static inline double sqf(double x) { return x*x; }
 
 void points2d_rordered_free(points2d_rordered_t *p) {
   free(p->rs);
@@ -80,6 +82,107 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig,
   return p;
 }
 
+
+static inline double pr2(const point2d p) {
+  return sqf(p.x) + sqf(p.y);
+}
+
+static inline double prn(const point2d p) {
+  return sqrt(pr2(p));
+}
+
+static int point2d_cmp_by_r2(const void *p1, const void *p2) {
+  const point2d *z1 = (point2d *) p1, *z2 = (point2d *) p2; 
+  double dif = pr2(*z1) - pr2(*z2);
+  if(dif > 0) return 1;
+  else if(dif < 0) return -1;
+  else return 0;
+}
+
+static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, const size_t nmemb,
+    const double rtol, const double atol, bool copybase)
+{
+  // TODO should the rtol and atol relate to |r| or r**2? (Currently: |r|)
+  assert(rtol >= 0);
+  assert(atol >= 0);
+  points2d_rordered_t *p = malloc(sizeof(points2d_rordered_t));
+  if(nmemb == 0) {
+    p->nrs = 0;
+    p->rs = NULL;
+    p->base = NULL;
+    p->r_offsets = NULL;
+    return p;
+  }
+
+  if (copybase) {
+    p->base = malloc(nmemb * sizeof(point2d));
+    memcpy(p->base, orig_base, nmemb * sizeof(point2d));
+  } else 
+    p->base = orig_base;
+
+  qsort(p->base, nmemb, sizeof(point2d), point2d_cmp_by_r2);
+
+  // first pass: determine the number of "different" r's.
+  size_t rcount = 0;
+  double rcur = -INFINITY;
+  double rcurmax = -INFINITY;
+  for (size_t i = 0; i < nmemb; ++i)
+    if ((rcur = prn(p->base[i])) > rcurmax) {
+      ++rcount;
+      rcurmax = rcur * (1 + rtol) + atol;
+    }
+
+  p->nrs = rcount;
+  // TODO check malloc return values
+  p->rs = malloc(rcount * sizeof(double));
+  p->r_offsets = malloc((rcount+1) * sizeof(ptrdiff_t));
+
+
+  // second pass: fill teh rs;
+  size_t ri = 0;
+  size_t rcurcount = 0;
+  rcur = prn(p->base[0]);
+  rcurmax = rcur * (1 + rtol) + atol;
+  double rcursum = 0;
+  p->r_offsets[0] = 0;
+  for (size_t i = 0; i < nmemb; ++i) {
+    rcur = prn(p->base[i]);
+    if (rcur > rcurmax) {
+      p->rs[ri] = rcursum / rcurcount; // average of the accrued r's within tolerance
+      ++ri;
+      p->r_offsets[ri] = i; //r_offsets[ri-1] + rcurcount (is the same)
+      rcurcount = 0;
+      rcursum = 0;
+      rcurmax = rcur * (1 + rtol) + atol;
+    } 
+    rcursum += rcur;
+    ++rcurcount;
+  }
+  p->r_offsets[rcount] = nmemb;
+
+  return p;
+}
+
+points2d_rordered_t *points2d_rordered_frompoints(const point2d *orig_base, const size_t nmemb,
+    const double rtol, const double atol) 
+{
+  return points2d_rordered_frompoints_c((point2d *)orig_base, nmemb, rtol, atol, true);
+}
+
+points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig, const  point2d shift,
+    double rtol, double atol)
+{
+  size_t n = (orig->nrs > 0) ?
+      orig->r_offsets[orig->nrs] - orig->r_offsets[0] : 0;
+
+  point2d * shifted = malloc(n * sizeof(point2d));
+
+  for(size_t i = 0; i < n; ++i) 
+    shifted[i] = cart2_add(orig->base[i], shift);
+
+  return points2d_rordered_frompoints_c(shifted, n,rtol, atol, false);
+}
+
 /*
  * EQUILATERAL TRIANGULAR LATTICE
  */

tests/lattice/2d_triangular_lattice.c

@@ -1,5 +1,6 @@
 #include <qpms/lattices.h>
 #include <stdio.h>
+#include <float.h>
 
 void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
   FILE *f = fopen(filename, "w");
@@ -65,6 +66,10 @@ int main() {
   dump_points2d_rordered(&(h->ps), "hex_h_empty.out");
   honeycomb_lattice_gen_extend_to_steps(h, 7);
   dump_points2d_rordered(&(h->ps), "hex_h_s7.out");
+  point2d shift = {0, h->h};
+  p = points2d_rordered_shift(&(h->ps), shift, DBL_EPSILON, h->h * DBL_EPSILON);
+  dump_points2d_rordered(p, "hex_h_s7_shifted.out");
+  points2d_rordered_free(p);
   honeycomb_lattice_gen_extend_to_steps(h, 120);
   dump_points2d_rordered(&(h->ps), "hex_h_s120.out");
   honeycomb_lattice_gen_free(h);