oldtests/ewaldshift.c

   1 // c99 -ggdb -Wall -I ../ ewaldshift.c ../qpms/ewald.c ../qpms/ewaldsf.c  ../qpms/lattices2d.c -lgsl -lm -lblas
   2
   3 // implementation of the [LT(4.16)] test
   4 #include <math.h>
   5 #define M_SQRTPI 1.7724538509055160272981674833411452
   6 #include <qpms/ewald.h>
   7 #include <qpms/tiny_inlines.h>
   8 #include <qpms/indexing.h>
   9 #include <stdlib.h>
  10 #include <stdio.h>
  11 #include <float.h>
  12 #include <gsl/gsl_sf_legendre.h>
  13 #include <gsl/gsl_const_mksa.h>
  14
  15
  16 static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
  17 static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  18 static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  19
  20
  21 typedef struct ewaldtest_triang_params {
  22   qpms_l_t lMax;
  23   point2d beta;
  24   point2d particle_shift;
  25   double k;
  26   double a;
  27   double eta;
  28   double maxR;
  29   double maxK;
  30   double csphase;
  31   TriangularLatticeOrientation orientation;
  32 } ewaldtest_triang_params;
  33
  34 typedef struct ewaldtest_triang_results {
  35   ewaldtest_triang_params p;
  36   complex double *sigmas_short,
  37           *sigmas_long,
  38           sigma0,
  39           *sigmas_total;
  40   double *err_sigmas_short,
  41     *err_sigmas_long,
  42     err_sigma0,
  43     *err_sigmas_total;
  44   complex double *regsigmas_416;
  45 } ewaldtest_triang_results;
  46
  47
  48 ewaldtest_triang_params paramslist[] = {
  49 // lMax, beta, shift, k, a, eta, maxR, maxK, csphase, orientation
  50   {3, {2.7/576e-9, 1/576e-9}, {576e-9/2,0}, 4.1*M_PI/576e-9, 576e-9, 1/576e-9, 40*576e-9, 40/576e-9, 1., TRIANGULAR_VERTICAL},
  51   {3, {2.7/576e-9, 1/576e-9}, {576e-9/2,0}, 4.1*M_PI/576e-9, 576e-9, 2/576e-9, 40*576e-9, 40/576e-9, 1., TRIANGULAR_VERTICAL},
  52   {3, {2.7/576e-9, 1/576e-9}, {576e-9/2,0}, 4.1*M_PI/576e-9, 576e-9, 3/576e-9, 40*576e-9, 40/576e-9, 1., TRIANGULAR_VERTICAL},
  53   {3, {2.7/576e-9, 1/576e-9}, {576e-9/2,0}, 4.1*M_PI/576e-9, 576e-9, 4/576e-9, 40*576e-9, 40/576e-9, 1., TRIANGULAR_VERTICAL},
  54 /*
  55  { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  56   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  57   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  58   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  59
  60   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  61   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  62   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  63   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  64
  65   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  66   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  67   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  68   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  69 */
  70   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  71   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  72   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  73   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  74
  75   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  76   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  77   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  78   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  79 /*
  80   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  81   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  82   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  83   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  84
  85   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  86   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  87   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  88   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  89
  90   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  91   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  92   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  93   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  94
  95   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  96   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  97   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  98   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  99
 100   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 101   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 102   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 103   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 104 */
 105   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 106   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 107   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 108   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 109
 110   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 111   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 112   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 113   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 114
 115   // Poloviční posun oproti přodchozímu
 116   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 117   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 118   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 119   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 120
 121   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 122   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 123   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 124   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 125
 126   // miniposun
 127   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 128   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 129   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 130   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 131
 132   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 133   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 134   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 135   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 136
 137   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 138   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 139   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 140   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 141
 142   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 143   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 144   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 145   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 146
 147
 148
 149
 150 /*
 151   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 152   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 153   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 154   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 155 */
 156
 157
 158
 159
 160 // end:
 161 //  { 0,  {0, 0}, 0, 0, 0, 0, 0, 0, 0}
 162 };
 163
 164 void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
 165   free(r->sigmas_short);
 166   free(r->sigmas_long);
 167   free(r->sigmas_total);
 168   free(r->err_sigmas_long);
 169   free(r->err_sigmas_total);
 170   free(r->err_sigmas_short);
 171   free(r->regsigmas_416);
 172   free(r);
 173 }
 174
 175
 176 void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
 177   FILE *f = fopen(filename, "w");
 178   for (size_t i = 0; i < ps->nrs; ++i) {
 179     fprintf(f, "# r = %.16g\n", ps->rs[i]);
 180     for (ptrdiff_t j = ps->r_offsets[i]; j < ps->r_offsets[i+1]; ++j)
 181       fprintf(f, "%.16g %.16g\n", ps->base[j].x, ps->base[j].y);
 182   }
 183   fclose(f);
 184 }
 185
 186
 187 static inline double san(double x) {
 188   return fabs(x) < 1e-13 ? 0 : x;
 189 }
 190
 191 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
 192
 193 int main() {
 194   gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
 195   for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) {
 196     ewaldtest_triang_params p = paramslist[i];
 197     ewaldtest_triang_results *r = ewaldtest_triang(p);
 198     // TODO print per-test header here
 199     printf("===============================\n");
 200     printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), ps = (%g,%g), csphase = %g\n",
 201         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.particle_shift.x, p.particle_shift.y, p.csphase);
 202     printf("sigma0: %.16g%+.16gj\n", creal(r->sigma0), cimag(r->sigma0));
 203     for (qpms_l_t n = 0; n <= p.lMax; ++n) {
 204       for (qpms_m_t m = -n; m <= n; ++m){
 205         if ((m+n)%2) continue;
 206         qpms_y_t y = qpms_mn2y_sc(m,n);
 207         qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
 208         // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
 209         printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n"
 210             //"| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n"
 211             ,
 212             y, n, m, creal(san(r->sigmas_total[y])), san(cimag(r->sigmas_total[y])),
 213             r->err_sigmas_total[y],
 214             san(creal(r->sigmas_long[y])), san(cimag(r->sigmas_long[y])),
 215             r->err_sigmas_long[y],
 216             san(creal(r->sigmas_short[y])), san(cimag(r->sigmas_short[y])),
 217             r->err_sigmas_short[y]
 218             //san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])),
 219             //san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])),
 220             //san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj]))
 221         );
 222       }
 223     }
 224     ewaldtest_triang_results_free(r);
 225   }
 226   return 0;
 227 }
 228
 229
 230 int ewaldtest_counter = 0;
 231
 232
 233 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
 234   const double a = p.a; //const double a = p.h * sqrt(3);
 235
 236   const double A = sqrt(3) * a * a / 2.; // unit cell size
 237   const double K_len = 4*M_PI/a/sqrt(3); // reciprocal vector length
 238
 239
 240   ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results));
 241   results->p = p;
 242
 243   triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, p.orientation, true, 0); // N.B. orig is included
 244   triangular_lattice_gen_extend_to_r(Rlg, p.maxR + a);
 245   triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0);
 246   triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len);
 247
 248   point2d *Rpoints = Rlg->ps.base; //point2d  *Kpoints = Klg->ps.base;
 249   size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs],
 250          nK = Klg->ps.r_offsets[Klg->ps.nrs];
 251
 252   point2d particle_shift = p.particle_shift;
 253   point2d minus_ps = {-particle_shift.x, -particle_shift.y};
 254   point2d Rpoints_plus_shift[nR];
 255   for(size_t i = 0; i < nR; ++i){
 256     Rpoints_plus_shift[i].x = Rpoints[i].x - particle_shift.x;
 257     Rpoints_plus_shift[i].y = Rpoints[i].y - particle_shift.y;
 258   }
 259
 260   qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax);
 261
 262   results->sigmas_short = malloc(sizeof(complex double)*nelem_sc);
 263   results->sigmas_long = malloc(sizeof(complex double)*nelem_sc);
 264   results->sigmas_total = malloc(sizeof(complex double)*nelem_sc);
 265   results->err_sigmas_short = malloc(sizeof(double)*nelem_sc);
 266   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
 267   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 268
 269   qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 270
 271   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
 272       8*DBL_EPSILON, 8*DBL_EPSILON);
 273
 274   char filename[BUFSIZ];
 275   sprintf(filename, "betalattice_%d.out", ewaldtest_counter);
 276   dump_points2d_rordered(Kpoints_plus_beta, filename);
 277
 278
 279   if (0!=ewald32_sigma_long_shiftedpoints(results->sigmas_long,
 280         results->err_sigmas_long, c, p.eta, p.k, A,
 281         nK, Kpoints_plus_beta->base,
 282         p.beta,
 283         particle_shift /*minus_ps*/ ))
 284     abort();
 285   if (0!=ewald32_sigma_short_shiftedpoints(
 286         results->sigmas_short, results->err_sigmas_short, c,
 287         p.eta, p.k,
 288         nR, Rpoints_plus_shift, p.beta, particle_shift))
 289     abort();
 290   if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
 291     abort();
 292   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 293     results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y];
 294     results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y];
 295   }
 296   results->sigmas_total[0] += results->sigma0;
 297   results->err_sigmas_total[0] += results->err_sigma0;
 298
 299   // Now calculate the reference values [LT(4.16)]
 300   results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double));
 301   results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)];
 302
 303   {
 304     double legendres[gsl_sf_legendre_array_n(p.lMax)];
 305     points2d_rordered_t sel =
 306       points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false);
 307     if (0 != sel.nrs)
 308     {
 309       point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0];
 310       size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0];
 311       for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) {
 312         point2d beta_pq = beta_pq_lessthan_k[i];
 313         double rbeta_pq = cart2norm(beta_pq);
 314         double arg_pq = atan2(beta_pq.y, beta_pq.x);
 315         double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq);
 316         if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
 317               p.lMax, denom/p.k, p.csphase, legendres) != 0)
 318           abort();
 319         for (qpms_y_t y = 0; y < nelem_sc; ++y) {
 320           qpms_l_t n; qpms_m_t m;
 321           qpms_y2mn_sc_p(y, &m, &n);
 322           if ((m+n)%2 != 0)
 323             continue;
 324           complex double eimf = cexp(I*m*arg_pq);
 325           results->regsigmas_416[y] +=
 326             4*M_PI*ipow(n)/p.k/A
 327             * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m)
 328             / denom;
 329         }
 330       }
 331     }
 332   }
 333
 334   points2d_rordered_free(Kpoints_plus_beta);
 335   qpms_ewald3_constants_free(c);
 336   triangular_lattice_gen_free(Klg);
 337   triangular_lattice_gen_free(Rlg);
 338   ++ewaldtest_counter;
 339   return results;
 340 }