oldtests/ewaldshift_3g.c

   1 // c99 -o ewaldshift_3g -ggdb -Wall -I ../ ewaldshift_3g.c ../qpms/ewald.c ../qpms/ewaldsf.c  ../qpms/lattices2d.c ../qpms/latticegens.c -lgsl -lm -lblas
   2
   3 // implementation of the [LT(4.16)] test
   4 #include <math.h>
   5 #define M_SQRTPI 1.7724538509055160272981674833411452
   6 #define M_SQRT3 1.7320508075688772935274463415058724
   7 #include <qpms/ewald.h>
   8 #include <qpms/tiny_inlines.h>
   9 #include <qpms/indexing.h>
  10 #include <stdlib.h>
  11 #include <stdio.h>
  12 #include <float.h>
  13 #include <gsl/gsl_sf_legendre.h>
  14 typedef struct ewaldtest_triang_params {
  15   qpms_l_t lMax;
  16   point2d beta;
  17   point2d particle_shift;
  18   double k;
  19   double a;
  20   double eta;
  21   double maxR;
  22   double maxK;
  23   double csphase;
  24   TriangularLatticeOrientation orientation;
  25 } ewaldtest_triang_params;
  26
  27 typedef struct ewaldtest_triang_results {
  28   ewaldtest_triang_params p;
  29   complex double *sigmas_short,
  30           *sigmas_long,
  31           sigma0,
  32           *sigmas_total;
  33   double *err_sigmas_short,
  34     *err_sigmas_long,
  35     err_sigma0,
  36     *err_sigmas_total;
  37   complex double *regsigmas_416;
  38 } ewaldtest_triang_results;
  39
  40
  41 /*
  42 const double a = 582e-9;
  43 const double inv_a = 4*M_PI/a/M_SQRT3;
  44 const double Klen = 4*M_PI/a/3;
  45 */
  46
  47 #define AA (582.e-9)
  48 #define INV_A (4*M_PI/AA/M_SQRT3)
  49 #define KLEN (4*M_PI/AA/3)
  50
  51 ewaldtest_triang_params paramslist[] = {
  52 // lMax, beta, shift, k, a, eta, maxR, maxK, csphase, orientation
  53 /*
  54  { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  55   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  56   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  57   { 2, {2.7, 1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  58
  59   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  60   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  61   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  62   { 2, {1.1, 1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  63
  64   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  65   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  66   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  67   { 2, {1.1, 1}, {0.5,0.}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  68 */
  69
  70 //
  71   { 3, {0,4198609.6394310603}, {0,0}, 11255786.828366444, 9.9766126515967311e-07, 29088820.866572164, 20*9.9766126515967311e-07, 20*7272205.21664304, 1., TRIANGULAR_VERTICAL},
  72   { 3, {0,4198609.6394310603}, {0,0}, 11255786.828366444, 9.9766126515967311e-07, 0.18*29088820.866572164, 20*9.9766126515967311e-07, 20*7272205.21664304, 1., TRIANGULAR_VERTICAL},
  73   { 3, {0,4198609.6394310603}, {0,0}, 11255786.828366444, 9.9766126515967311e-07, 0.13*29088820.866572164, 20*9.9766126515967311e-07, 20*7272205.21664304, 1., TRIANGULAR_VERTICAL},
  74   { 3, {0,4198609.6394310603}, {0,0}, 11255786.828366444, 9.9766126515967311e-07, 0.07*29088820.866572164, 20*9.9766126515967311e-07, 20*7272205.21664304, 1., TRIANGULAR_VERTICAL},
  75   { 3, {0,4198609.6394310603}, {0,0}, 11255786.828366444, 9.9766126515967311e-07, 0.03*29088820.866572164, 20*9.9766126515967311e-07, 20*7272205.21664304, 1., TRIANGULAR_VERTICAL},
  76
  77 //  { 3, {0,KLEN}, {0,0}, 2.62 * 4 * M_PI/3/AA, AA, 0.5 / AA, 20*AA, 160/AA, 1., TRIANGULAR_VERTICAL},
  78   { 3, {0,KLEN}, {0,0}, 2.62 * 4 * M_PI/3/AA, AA, 2.5 / AA, 20*AA, 160/AA, 1., TRIANGULAR_VERTICAL},
  79   { 3, {0,KLEN}, {0,0}, 2.62 * 4 * M_PI/3/AA, AA, 4.5 / AA, 20*AA, 160/AA, 1., TRIANGULAR_VERTICAL},
  80   { 3, {0,KLEN}, {0,0}, 2.62 * 4 * M_PI/3/AA, AA, 6.5 / AA, 20*AA, 160/AA, 1., TRIANGULAR_VERTICAL},
  81   { 3, {0,KLEN}, {0,0}, 2.62 * 4 * M_PI/3/AA, AA, 8.5 / AA, 20*AA, 160/AA, 1., TRIANGULAR_VERTICAL},
  82 /*
  83   { 3, {0,Klen}, {0,0}, 2.62 * 4 * M_PI/3/a, a, 0.5 / a, 20*a, 2*M_PI*160/a, 1., TRIANGULAR_VERTICAL},
  84   { 3, {0,Klen}, {0,0}, 2.62 * 4 * M_PI/3/a, a, 2.5 / a, 20*a, 2*M_PI*160/a, 1., TRIANGULAR_VERTICAL},
  85   { 3, {0,Klen}, {0,0}, 2.62 * 4 * M_PI/3/a, a, 4.5 / a, 20*a, 2*M_PI*160/a, 1., TRIANGULAR_VERTICAL},
  86   { 3, {0,Klen}, {0,0}, 2.62 * 4 * M_PI/3/a, a, 6.5 / a, 20*a, 2*M_PI*160/a, 1., TRIANGULAR_VERTICAL},
  87   { 3, {0,Klen}, {0,0}, 2.62 * 4 * M_PI/3/a, a, 8.5 / a, 20*a, 2*M_PI*160/a, 1., TRIANGULAR_VERTICAL},
  88 */
  89
  90
  91   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  92   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  93   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  94   { 2, {1.1, 2.1}, {0.5,0.1325}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  95 /*
  96   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  97   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  98   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  99   { 2, {0, 3.1}, {0.5,0}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 100
 101   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 102   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 103   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 104   { 2, {0, 1.1}, {0.5,0}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 105
 106   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 107   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 108   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 109   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 110
 111   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 112   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 113   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 114   { 2, {1.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 115
 116   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 117   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 118   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 119   { 2, {3.1,0}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 120 */
 121   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 122   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 123   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 124   { 2, {3.1*0.5,-3.1*0.8}, {0.8,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 125
 126   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 127   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 128   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 129   { 2, {1.1*0.5,-1.1*0.8}, {0.8,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 130
 131   // Poloviční posun oproti přodchozímu
 132   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 133   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 134   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 135   { 2, {3.1*0.5,-3.1*0.8}, {0.4,0.25}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 136
 137   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 138   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 139   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 140   { 2, {1.1*0.5,-1.1*0.8}, {0.4,0.25}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 141
 142   // miniposun
 143   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 144   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 145   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 146   { 3, {3.1*0.5,-3.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 147
 148   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 149   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 150   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 151   { 3, {1.1*0.5,-1.1*0.8}, {0.004,0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 152
 153   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 154   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 155   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 156   { 3, {3.1*0.5,-3.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 157
 158   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 159   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 160   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 161   { 3, {1.1*0.5,-1.1*0.8}, {-0.004,-0.0025}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 162
 163
 164
 165
 166 /*
 167   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 168   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 169   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 170   { 2, {0, 3.1}, {0,0.5}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
 171 */
 172
 173
 174
 175
 176 // end:
 177 //  { 0,  {0, 0}, 0, 0, 0, 0, 0, 0, 0}
 178 };
 179
 180 void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
 181   free(r->sigmas_short);
 182   free(r->sigmas_long);
 183   free(r->sigmas_total);
 184   free(r->err_sigmas_long);
 185   free(r->err_sigmas_total);
 186   free(r->err_sigmas_short);
 187   free(r->regsigmas_416);
 188   free(r);
 189 }
 190
 191 static inline double san(double x) {
 192   return fabs(x) < 1e-13 ? 0 : x;
 193 }
 194
 195 ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p);
 196
 197 int main() {
 198   gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
 199   for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) {
 200     ewaldtest_triang_params p = paramslist[i];
 201     ewaldtest_triang_results *r = ewaldtest_triang_3g(p);
 202     // TODO print per-test header here
 203     printf("===============================\n");
 204     printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), ps = (%g,%g), csphase = %g\n",
 205         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);
 206     printf("sigma0: %.16g%+.16gj\n", creal(r->sigma0), cimag(r->sigma0));
 207     for (qpms_l_t n = 0; n <= p.lMax; ++n) {
 208       for (qpms_m_t m = -n; m <= n; ++m){
 209         if ((m+n)%2) continue;
 210         qpms_y_t y = qpms_mn2y_sc(m,n);
 211         qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
 212         // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
 213         printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n"
 214             //"| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n"
 215             ,
 216             y, n, m, creal(san(r->sigmas_total[y])), san(cimag(r->sigmas_total[y])),
 217             r->err_sigmas_total[y],
 218             san(creal(r->sigmas_long[y])), san(cimag(r->sigmas_long[y])),
 219             r->err_sigmas_long[y],
 220             san(creal(r->sigmas_short[y])), san(cimag(r->sigmas_short[y])),
 221             r->err_sigmas_short[y]
 222             //san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])),
 223             //san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])),
 224             //san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj]))
 225         );
 226       }
 227     }
 228     ewaldtest_triang_results_free(r);
 229   }
 230   return 0;
 231 }
 232
 233
 234 int ewaldtest_counter = 0;
 235
 236
 237 ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) {
 238   const double a = p.a; //const double a = p.h * sqrt(3);
 239   cart3_t beta3 = cart22cart3xy(p.beta);
 240   cart3_t particle_shift3 = cart22cart3xy(p.particle_shift);
 241
 242   cart2_t b1, b2, rb1, rb2;
 243   if (p.orientation == TRIANGULAR_VERTICAL) {
 244     b1.x = 0;
 245     b1.y = a;
 246     b2.x = a * M_SQRT3 * .5;
 247     b2.y = a * .5;
 248   } else {
 249     b1.x = a;
 250     b1.y = 0;
 251     b2.x = a * .5;
 252     b2.y = a * M_SQRT3 * .5;
 253   }
 254   if (QPMS_SUCCESS != l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2))
 255     abort();
 256
 257   const double A = l2d_unitcell_area(b1, b2); // sqrt(3) * a * a / 2.; // unit cell size
 258   const double K_len = cart2norm(rb1); //4*M_PI/a/sqrt(3); // reciprocal vector length
 259
 260
 261   ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results));
 262   results->p = p;
 263
 264   //triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, p.orientation, true, 0); // N.B. orig is included
 265   //triangular_lattice_gen_extend_to_r(Rlg, p.maxR + a);
 266   //triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0);
 267   //triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len);
 268
 269   //point2d *Rpoints = Rlg->ps.base;
 270   //size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs];
 271
 272   /*if (fabs(p.particle_shift.x) ==0 && fabs(p.particle_shift.y) == 0) {
 273     points2d_rordered_t Rpos = points2d_rordered_annulus(&(Rlg->ps), 0, false, INFINITY, false);
 274     Rpoints = Rpos.base + Rpos.r_offsets[0];
 275     nR = Rpos.r_offsets[Rpos.nrs] - Rpos.r_offsets[0];
 276   }*/
 277
 278   //point2d  *Kpoints = Klg->ps.base;
 279   //size_t nK = Klg->ps.r_offsets[Klg->ps.nrs];
 280
 281   /*
 282   point2d particle_shift = p.particle_shift;
 283   point2d minus_ps = {-particle_shift.x, -particle_shift.y};
 284   point2d Rpoints_plus_shift[nR];
 285   for(size_t i = 0; i < nR; ++i){
 286     Rpoints_plus_shift[i].x = Rpoints[i].x - particle_shift.x;
 287     Rpoints_plus_shift[i].y = Rpoints[i].y - particle_shift.y;
 288   }
 289   */
 290
 291
 292   // skip zeroth point if it coincides with origin
 293   bool include_origin = !(fabs(p.particle_shift.x) == 0
 294       && fabs(p.particle_shift.y) == 0);
 295
 296   PGen Rlgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, CART2_ZERO, 0, include_origin, p.maxR + a, false);
 297   //PGen Rlgen_plus_shift = PGen_xyWeb_new(b1, b2, BASIS_RTOL, cart2_scale(-1 /* CHECKSIGN */, particle_shift2), 0, include_origin, p.maxR + a, false);
 298   PGen Klgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, CART2_ZERO, 0, true, p.maxK + K_len, false);
 299   //PGen Klgen_plus_beta = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, beta2, 0, true, p.maxK + K_len, false);
 300
 301   qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax);
 302
 303   results->sigmas_short = malloc(sizeof(complex double)*nelem_sc);
 304   results->sigmas_long = malloc(sizeof(complex double)*nelem_sc);
 305   results->sigmas_total = malloc(sizeof(complex double)*nelem_sc);
 306   results->err_sigmas_short = malloc(sizeof(double)*nelem_sc);
 307   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
 308   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 309
 310   qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 311
 312   //points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
 313   //    8*DBL_EPSILON, 8*DBL_EPSILON);
 314
 315   //char filename[BUFSIZ];
 316   //sprintf(filename, "betalattice_%d.out", ewaldtest_counter);
 317   //dump_points2d_rordered(Kpoints_plus_beta, filename);
 318
 319
 320   if (0!=ewald3_sigma_long(results->sigmas_long,
 321         results->err_sigmas_long, c, p.eta, p.k, A,
 322         LAT_2D_IN_3D_XYONLY, &Klgen, false, beta3, particle_shift3))
 323       abort();
 324 #if 0
 325   if (0!=ewald32_sigma_long_points_and_shift(results->sigmas_long,
 326         results->err_sigmas_long, c, p.eta, p.k, A,
 327         nK, Kpoints,
 328         p.beta,
 329         particle_shift /*minus_ps*/ ))
 330     abort();
 331 #endif
 332   if (0!=ewald3_sigma_short(
 333     results->sigmas_short, results->err_sigmas_short, c,
 334     p.eta, p.k, LAT_2D_IN_3D_XYONLY, &Rlgen, false,  beta3, particle_shift3))
 335     abort();
 336   /*if (0!=ewald32_sigma_short_points_and_shift(
 337         results->sigmas_short, results->err_sigmas_short, c,
 338         p.eta, p.k,
 339         nR, Rpoints, p.beta, particle_shift))
 340     abort();*/
 341   //if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
 342   if (0!=ewald3_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
 343     abort();
 344   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 345     results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y];
 346     results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y];
 347   }
 348   results->sigmas_total[0] += results->sigma0;
 349   results->err_sigmas_total[0] += results->err_sigma0;
 350
 351   // Now calculate the reference values [LT(4.16)]
 352   results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double));
 353   results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)];
 354
 355 #if 0 // not yet implemented for the new API
 356   {
 357     double legendres[gsl_sf_legendre_array_n(p.lMax)];
 358     points2d_rordered_t sel =
 359       points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false);
 360     if (0 != sel.nrs)
 361     {
 362       point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0];
 363       size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0];
 364       for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) {
 365         point2d beta_pq = beta_pq_lessthan_k[i];
 366         double rbeta_pq = cart2norm(beta_pq);
 367         double arg_pq = atan2(beta_pq.y, beta_pq.x);
 368         double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq);
 369         if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
 370               p.lMax, denom/p.k, p.csphase, legendres) != 0)
 371           abort();
 372         for (qpms_y_t y = 0; y < nelem_sc; ++y) {
 373           qpms_l_t n; qpms_m_t m;
 374           qpms_y2mn_sc_p(y, &m, &n);
 375           if ((m+n)%2 != 0)
 376             continue;
 377           complex double eimf = cexp(I*m*arg_pq);
 378           results->regsigmas_416[y] +=
 379             4*M_PI*ipow(n)/p.k/A
 380             * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m)
 381             / denom;
 382         }
 383       }
 384     }
 385   }
 386 #else
 387   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 388           qpms_l_t n; qpms_m_t m;
 389           qpms_y2mn_sc_p(y, &m, &n);
 390           if ((m+n)%2 != 0)
 391             continue;
 392           results->regsigmas_416[y] = NAN;
 393   }
 394 #endif
 395
 396
 397   //points2d_rordered_free(Kpoints_plus_beta);
 398   qpms_ewald3_constants_free(c);
 399   //triangular_lattice_gen_free(Klg);
 400   //triangular_lattice_gen_free(Rlg);
 401   ++ewaldtest_counter;
 402   return results;
 403 }