Fix points2d_rordered_frompoints last layer r;

dump lattices intests/ewalds.c


Former-commit-id: 1addca6a1c240310e86213c283df4162b2536b87

qpms/lattices2d.c

@@ -66,9 +66,14 @@ points2d_rordered_t points2d_rordered_annulus(const points2d_rordered_t *orig,
   ptrdiff_t imin, imax;
   imin = points2d_rordered_locate_r(orig, minr);
   imax = points2d_rordered_locate_r(orig, maxr);
+  // TODO check 
+  if(imax >= orig->nrs) --imax;
+  if(imax < 0) goto nothing;
+  // END TODO
   if (!inc_minr && (orig->rs[imin] <= minr)) ++imin;
   if (!inc_maxr && (orig->rs[imax] >= maxr)) --imax;
   if (imax < imin) { // it's empty
+nothing:
     p.nrs = 0;
     p.base = NULL;
     p.rs = NULL;
@@ -148,7 +153,7 @@ static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, c
   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
+      p->rs[ri] = rcursum / (double) 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;
@@ -158,6 +163,7 @@ static points2d_rordered_t *points2d_rordered_frompoints_c(point2d *orig_base, c
     rcursum += rcur;
     ++rcurcount;
   }
+  p->rs[ri] = rcursum / (double) rcurcount;
   p->r_offsets[rcount] = nmemb;
 
   return p;
@@ -178,7 +184,7 @@ points2d_rordered_t *points2d_rordered_shift(const points2d_rordered_t *orig, co
   point2d * shifted = malloc(n * sizeof(point2d));
 
   for(size_t i = 0; i < n; ++i) 
-    shifted[i] = cart2_add(orig->base[i], shift);
+    shifted[i] = cart2_add(orig->base[i+orig->r_offsets[0]], shift);
 
   return points2d_rordered_frompoints_c(shifted, n,rtol, atol, false);
 }

tests/ewalds.c

@@ -37,7 +37,8 @@ typedef struct ewaldtest_triang_results {
 
 
 ewaldtest_triang_params paramslist[] = {
-  { 3, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_HORIZONTAL},
+// lMax, beta, k, a, eta, maxR, maxK, csphase, orientation
+#if 0
   { 3, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
   { 3, {1.1, 0.23}, 2.3, 0.97, 0.5,  30, 30, -1., TRIANGULAR_VERTICAL},
   { 3, {1.1, 0.23}, 2.3, 0.97, 0.9,  30, 30, 1., TRIANGULAR_VERTICAL},
@@ -50,6 +51,17 @@ ewaldtest_triang_params paramslist[] = {
   { 6, {1.1, 0.23}, 2.3, 0.97, 4.5,  40, 40, 1., TRIANGULAR_VERTICAL},
   { 6, {1.1, 0.23}, 2.3, 0.97, 2.3,  100, 100, 1., TRIANGULAR_VERTICAL},
   { 6, {1.1, 0.23}, 2.3, 0.97, 2.9,  100, 100, 1., TRIANGULAR_VERTICAL},
+#endif
+  { 2, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 10, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 5, 1., TRIANGULAR_VERTICAL},
+  { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL},
+
 // end:
 //  { 0,  {0, 0}, 0, 0, 0, 0, 0, 0, 0}
 };
@@ -65,6 +77,18 @@ void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
   free(r);
 }
 
+
+void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
+  FILE *f = fopen(filename, "w");
+  for (size_t i = 0; i < ps->nrs; ++i) {
+    fprintf(f, "# r = %.16g\n", ps->rs[i]);
+    for (ptrdiff_t j = ps->r_offsets[i]; j < ps->r_offsets[i+1]; ++j)
+      fprintf(f, "%.16g %.16g\n", ps->base[j].x, ps->base[j].y);
+  }
+  fclose(f);
+}
+
+
 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
 
 int main() {
@@ -74,8 +98,8 @@ int main() {
     ewaldtest_triang_results *r = ewaldtest_triang(p);
     // TODO print per-test header here
     printf("===============================\n");
-    printf("Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), csphase = %g\n",
+    printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), csphase = %g\n",
-        p.maxK, p.maxR, p.lMax, p.eta, p.k, p.beta.x, p.beta.y, p.csphase);
+        p.a, 4*M_PI/sqrt(3)/p.a, p.maxK, p.maxR, p.lMax, p.eta, p.k, p.beta.x, p.beta.y, p.csphase);
     printf("sigma0: %.16g%+.16gj\n", creal(r->sigma0), cimag(r->sigma0));    
     for (qpms_l_t n = 0; n <= p.lMax; ++n) {
       for (qpms_m_t m = -n; m <= n; ++m){
@@ -101,6 +125,9 @@ int main() {
 }
 
 
+int ewaldtest_counter = 0;
+
+
 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   const double a = p.a; //const double a = p.h * sqrt(3);
  
@@ -134,6 +161,10 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
       8*DBL_EPSILON, 8*DBL_EPSILON);
 
+  char filename[BUFSIZ];
+  sprintf(filename, "betalattice_%d.out", ewaldtest_counter);
+  dump_points2d_rordered(Kpoints_plus_beta, filename);
+
   point2d particle_shift = {0,0}; // TODO make this a parameter
 
   if (0!=ewald32_sigma_long_shiftedpoints(results->sigmas_long,
@@ -192,5 +223,6 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   qpms_ewald32_constants_free(c);
   triangular_lattice_gen_free(Klg);
   triangular_lattice_gen_free(Rlg);
+  ++ewaldtest_counter;
   return results;
 }