diff --git a/qpms/ewald.c b/qpms/ewald.c
index d4d6a29..46e1c60 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -57,14 +57,14 @@ typedef enum {
                               * spherical harmonic. See notes/ewald.lyx.
                               */
 
-} ewald32_constants_option;
+} ewald3_constants_option;
 
-static const ewald32_constants_option type = EWALD32_CONSTANTS_AGNOSTIC;
+static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC;
 
-qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, const ewald32_constants_option type */,
+qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, const ewald3_constants_option type */,
     const int csphase)
 {
-  qpms_ewald32_constants_t *c = malloc(sizeof(qpms_ewald32_constants_t));
+  qpms_ewald3_constants_t *c = malloc(sizeof(qpms_ewald3_constants_t));
   //if (c == NULL) return NULL; // Do I really want to do this?
   c->lMax = lMax;
   c->nelem_sc = qpms_lMax2nelem_sc(lMax);
@@ -160,7 +160,7 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
   return c;
 }
 
-void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) {
+void qpms_ewald3_constants_free(qpms_ewald3_constants_t *c) {
   free(c->legendre0);
   free(c->legendre_plus1);
   free(c->legendre_minus1);
@@ -175,7 +175,7 @@ void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) {
 
 
 int ewald3_sigma0(complex double *result, double *err,
-    const qpms_ewald32_constants_t *c,
+    const qpms_ewald3_constants_t *c,
     const double eta, const complex double k)
 {
   qpms_csf_result gam;
@@ -189,198 +189,10 @@ int ewald3_sigma0(complex double *result, double *err,
   return 0;
 }
 
-int ewald32_sigma0(complex double *result, double *err,
-    const qpms_ewald32_constants_t *c,
-    const double eta, const double k)
-{
-  return ewald3_sigma0(result, err, c, eta, k);
-}
-
-
-
-int ewald32_sigma_long_shiftedpoints (
-    complex double *target, // must be c->nelem_sc long
-    double *err,
-    const qpms_ewald32_constants_t *c,
-    const double eta, const double k, const double unitcell_area,
-    const size_t npoints, const point2d *Kpoints_plus_beta,
-    const point2d beta, // not needed
-    const point2d particle_shift // target - src
-    ) 
-{
-  const qpms_y_t nelem_sc = c->nelem_sc;
-  const qpms_l_t lMax = c->lMax;
-  
-  // Manual init of the ewald summation targets
-  complex double *target_c = calloc(nelem_sc, sizeof(complex double));
-  memset(target, 0, nelem_sc * sizeof(complex double));
-  double *err_c = NULL;
-  if (err) {
-    err_c = calloc(nelem_sc, sizeof(double));
-    memset(err, 0, nelem_sc * sizeof(double));
-  }
-
-  const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC)
-  assert(commonfac > 0);
-
-  // space for Gamma_pq[j]'s
-  qpms_csf_result Gamma_pq[lMax/2+1];
-
-  // CHOOSE POINT BEGIN
-  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
-    const point2d beta_pq = Kpoints_plus_beta[i];
-    const point2d K_pq = {beta_pq.x - beta.x, beta_pq.y - beta.y};
-    const double rbeta_pq = cart2norm(beta_pq);
-  // CHOOSE POINT END
-
-    const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!!
-    const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT
-
-    // R-DEPENDENT BEGIN
-    const complex double gamma_pq = lilgamma(rbeta_pq/k);
-    const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
-    for(qpms_l_t j = 0; j <= lMax/2; ++j) {
-      int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
-      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
-      if(creal(z) < 0) 
-        Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
-      if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
-    }
-    // R-DEPENDENT END
-    
-    // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
-    // and just fetched for each n, m pair
-    for(qpms_l_t n = 0; n <= lMax; ++n)
-      for(qpms_m_t m = -n; m <= n; ++m) {
-        if((m+n) % 2 != 0) // odd coefficients are zero.
-          continue;
-        const qpms_y_t y = qpms_mn2y_sc(m, n);
-        const complex double e_imalpha_pq = cexp(I*m*arg_pq);
-        complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
-        double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
-        assert((n-abs(m))/2 == c->s1_jMaxes[y]);
-        for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
-          complex double summand = pow(rbeta_pq/k, n-2*j) 
-            * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
-            * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
-            * c->s1_constfacs[y][j];
-          if(err) {
-            // FIXME include also other errors than Gamma_pq's relative error
-             kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand));
-          }
-          summand *= Gamma_pq[j].val; // GGG
-          ckahanadd(&jsum, &jsum_c, summand);
-        }
-        jsum *= phasefac; // PFC
-        ckahanadd(target + y, target_c + y, jsum);
-        if(err) kahanadd(err + y, err_c + y, jsum_err);
-      }
-  } // END POINT LOOP
- 
-  free(err_c);
-  free(target_c);
-  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
-    target[y] *= commonfac;
-  if(err)
-    for(qpms_y_t y = 0; y < nelem_sc; ++y)
-      err[y] *= commonfac;
-  return 0;
-}
-
-int ewald32_sigma_long_points_and_shift (
-    complex double *target, // must be c->nelem_sc long
-    double *err,
-    const qpms_ewald32_constants_t *c,
-    const double eta, const double k, const double unitcell_area,
-    const size_t npoints, const point2d *Kpoints,
-    const point2d beta, 
-    const point2d particle_shift // target - src
-    ) 
-{
-  const qpms_y_t nelem_sc = c->nelem_sc;
-  const qpms_l_t lMax = c->lMax;
-  
-  // Manual init of the ewald summation targets
-  complex double *target_c = calloc(nelem_sc, sizeof(complex double));
-  memset(target, 0, nelem_sc * sizeof(complex double));
-  double *err_c = NULL;
-  if (err) {
-    err_c = calloc(nelem_sc, sizeof(double));
-    memset(err, 0, nelem_sc * sizeof(double));
-  }
-
-  const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC)
-  assert(commonfac > 0);
-
-  // space for Gamma_pq[j]'s
-  qpms_csf_result Gamma_pq[lMax/2+1];
-
-  // CHOOSE POINT BEGIN
-  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
-    // Only these following two lines differ from ewald32_sigma_long_points_and_shift()!!!  WTFCOMMENT?!
-    const point2d K_pq = Kpoints[i]; 
-    const point2d beta_pq = {K_pq.x + beta.x, K_pq.y + beta.y};
-    const double rbeta_pq = cart2norm(beta_pq);
-  // CHOOSE POINT END
-
-    const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!!
-    const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT
-
-    // R-DEPENDENT BEGIN
-    const complex double gamma_pq = lilgamma(rbeta_pq/k);
-    const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
-    for(qpms_l_t j = 0; j <= lMax/2; ++j) {
-      int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
-      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
-      if(creal(z) < 0) 
-        Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
-      if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
-    }
-    // R-DEPENDENT END
-    
-    // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
-    // and just fetched for each n, m pair
-    for(qpms_l_t n = 0; n <= lMax; ++n)
-      for(qpms_m_t m = -n; m <= n; ++m) {
-        if((m+n) % 2 != 0) // odd coefficients are zero.
-          continue;
-        const qpms_y_t y = qpms_mn2y_sc(m, n);
-        const complex double e_imalpha_pq = cexp(I*m*arg_pq);
-        complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
-        double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
-        assert((n-abs(m))/2 == c->s1_jMaxes[y]);
-        for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
-          complex double summand = pow(rbeta_pq/k, n-2*j) 
-            * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
-            * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
-            * c->s1_constfacs[y][j];
-          if(err) {
-            // FIXME include also other errors than Gamma_pq's relative error
-             kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand));
-          }
-          summand *= Gamma_pq[j].val; // GGG
-          ckahanadd(&jsum, &jsum_c, summand);
-        }
-        jsum *= phasefac; // PFC
-        ckahanadd(target + y, target_c + y, jsum);
-        if(err) kahanadd(err + y, err_c + y, jsum_err);
-      }
-  } // END POINT LOOP
- 
-  free(err_c);
-  free(target_c);
-  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
-    target[y] *= commonfac;
-  if(err)
-    for(qpms_y_t y = 0; y < nelem_sc; ++y)
-      err[y] *= commonfac;
-  return 0;
-}
-
 int ewald3_21_xy_sigma_long (
     complex double *target, // must be c->nelem_sc long
     double *err,
-    const qpms_ewald32_constants_t *c,
+    const qpms_ewald3_constants_t *c,
     const double eta, const complex double k,
     const double unitcell_volume /* with the corresponding lattice dimensionality */,
     const LatticeDimensionality latdim,
@@ -517,7 +329,7 @@ int ewald3_21_xy_sigma_long (
 int ewald3_1_z_sigma_long (
     complex double *target, // must be c->nelem_sc long
     double *err,
-    const qpms_ewald32_constants_t *c,
+    const qpms_ewald3_constants_t *c,
     const double eta, const complex double k,
     const double unitcell_volume /* length (periodicity) in this case */,
     const LatticeDimensionality latdim,
@@ -633,7 +445,7 @@ int ewald3_1_z_sigma_long (
 int ewald3_sigma_long (
     complex double *target, // must be c->nelem_sc long
     double *err,
-    const qpms_ewald32_constants_t *c,
+    const qpms_ewald3_constants_t *c,
     const double eta, const complex double k,
     const double unitcell_volume /* with the corresponding lattice dimensionality */,
     const LatticeDimensionality latdim,
@@ -749,139 +561,10 @@ static int ewald32_sr_integral_ck(double r, complex double k, int n, double eta,
   return retval;
 }
 
-int ewald32_sigma_short_shiftedpoints(
-    complex double *target, // must be c->nelem_sc long
-    double *err,
-    const qpms_ewald32_constants_t *c, // N.B. not too useful here
-    const double eta, const double k,
-    const size_t npoints, const point2d *Rpoints_plus_particle_shift,
-    const point2d beta,
-    const point2d particle_shift           // used only in the very end to multiply it by the phase
-    )
-{
-  const qpms_y_t nelem_sc = c->nelem_sc;
-  const qpms_l_t lMax = c->lMax;
-  gsl_integration_workspace *workspace = 
-    gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
-  
-  // Manual init of the ewald summation targets
-  complex double * const target_c = calloc(nelem_sc, sizeof(complex double));
-  memset(target, 0, nelem_sc * sizeof(complex double));
-  double *err_c = NULL;
-  if (err) {
-    err_c = calloc(nelem_sc, sizeof(double));
-    memset(err, 0, nelem_sc * sizeof(double));
-  }
-  
-
-// CHOOSE POINT BEGIN
-  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
-    const point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
-    const double r_pq_shifted = cart2norm(Rpq_shifted);
-// CHOOSE POINT END
-
-    const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT
-    const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
-    
-    for(qpms_l_t n = 0; n <= lMax; ++n) {
-      const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
-      const double R_pq_pown = pow(r_pq_shifted, n);
-      // TODO the integral here
-      double intres, interr;
-      int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta,
-          &intres, &interr, workspace);
-      if (retval) abort();
-      for (qpms_m_t m = -n; m <= n; ++m){
-        if((m+n) % 2 != 0) // odd coefficients are zero.
-          continue; // nothing needed, already done by memset
-        const complex double e_imf = cexp(I*m*Rpq_shifted_arg);
-        const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
-        const qpms_y_t y = qpms_mn2y_sc(m,n);
-        if(err)
-          kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
-              * interr)); // TODO include also other errors
-        ckahanadd(target + y, target_c + y,
-            prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m));
-      }
-
-    }
-  }
-
-  gsl_integration_workspace_free(workspace);
-  if(err) free(err_c);
-  free(target_c);
-  return 0;
-}
-
-int ewald32_sigma_short_points_and_shift(
-    complex double *target, // must be c->nelem_sc long
-    double *err,
-    const qpms_ewald32_constants_t *c, // N.B. not too useful here
-    const double eta, const double k,
-    const size_t npoints, const point2d *Rpoints,
-    const point2d beta,
-    const point2d particle_shift           // used only in the very end to multiply it by the phase
-    )
-{
-  const qpms_y_t nelem_sc = c->nelem_sc;
-  const qpms_l_t lMax = c->lMax;
-  gsl_integration_workspace *workspace = 
-    gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
-  
-  // Manual init of the ewald summation targets
-  complex double * const target_c = calloc(nelem_sc, sizeof(complex double));
-  memset(target, 0, nelem_sc * sizeof(complex double));
-  double *err_c = NULL;
-  if (err) {
-    err_c = calloc(nelem_sc, sizeof(double));
-    memset(err, 0, nelem_sc * sizeof(double));
-  }
-  
-
-// CHOOSE POINT BEGIN
-  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
-    //const point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
-    const point2d Rpq_shifted = cart2_add(Rpoints[i], cart2_scale(-1,particle_shift)); // CHECKSIGN!!!!
-    const double r_pq_shifted = cart2norm(Rpq_shifted);
-// CHOOSE POINT END
-
-    const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT
-    const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
-    
-    for(qpms_l_t n = 0; n <= lMax; ++n) {
-      const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
-      const double R_pq_pown = pow(r_pq_shifted, n);
-      // TODO the integral here
-      double intres, interr;
-      int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta,
-          &intres, &interr, workspace);
-      if (retval) abort();
-      for (qpms_m_t m = -n; m <= n; ++m){
-        if((m+n) % 2 != 0) // odd coefficients are zero.
-          continue; // nothing needed, already done by memset
-        const complex double e_imf = cexp(I*m*Rpq_shifted_arg);
-        const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
-        const qpms_y_t y = qpms_mn2y_sc(m,n);
-        if(err)
-          kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
-              * interr)); // TODO include also other errors
-        ckahanadd(target + y, target_c + y,
-            prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m));
-      }
-
-    }
-  }
-
-  gsl_integration_workspace_free(workspace);
-  if(err) free(err_c);
-  free(target_c);
-  return 0;
-}
-
 int ewald3_sigma_short( 
                 complex double *target, // must be c->nelem_sc long
                 double *err, // must be c->nelem_sc long or NULL
-                const qpms_ewald32_constants_t *c,
+                const qpms_ewald3_constants_t *c,
                 const double eta, const complex double k /* TODO COMPLEX */,
                 const LatticeDimensionality latdim, // apart from asserts and possible optimisations ignored, as the SR formula stays the same
                 PGen *pgen_R, const bool pgen_generates_shifted_points
@@ -1027,23 +710,3 @@ int ewald3_sigma_short(
   return 0;
 }
 
-
- 
-#if 0
-
-int ewald32_sigma_long_shiftedpoints_rordered(
-    complex double *target_sigmalr_y, // must be c->nelem_sc long
-    const qpms_ewald32_constants_t *c,
-    double eta, double k, double unitcell_area,
-    const points2d_rordered_t *Kpoints_plus_beta_rordered,
-    point2d particle_shift
-    );
-int ewald32_sigma_short_points_rordered(
-    complex double *target_sigmasr_y, // must be c->nelem_sc long
-    const qpms_ewald32_constants_t *c, // N.B. not too useful here
-    double eta, double k,
-    const points2d_rordered_t *Rpoints_plus_particle_shift_rordered,
-    point2d particle_shift    // used only in the very end to multiply it by the phase
-    );
-
-#endif
diff --git a/qpms/ewald.h b/qpms/ewald.h
index a4a5a56..ca794b3 100644
--- a/qpms/ewald.h
+++ b/qpms/ewald.h
@@ -76,10 +76,10 @@ typedef struct {
 					This is because I dont't actually consider this fixed in
 					translations.c */
 
-} qpms_ewald32_constants_t;
+} qpms_ewald3_constants_t;
 
-qpms_ewald32_constants_t *qpms_ewald32_constants_init(qpms_l_t lMax, int csphase);
-void qpms_ewald32_constants_free(qpms_ewald32_constants_t *);
+qpms_ewald3_constants_t *qpms_ewald3_constants_init(qpms_l_t lMax, int csphase);
+void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
 
 
 typedef struct { // as gsl_sf_result, but with complex val
@@ -144,14 +144,14 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
 // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space
 
 int ewald3_sigma0(complex double *result, double *err,
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, complex double k
 );
 
 int ewald3_sigma_short(
 		complex double *target_sigmasr_y, // must be c->nelem_sc long
 		double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		const double eta, const complex double k,
 		const LatticeDimensionality latdim, // apart from asserts and possible optimisations ignored, as the SR formula stays the same
 		PGen *pgen_R, const bool pgen_generates_shifted_points 
@@ -168,7 +168,7 @@ int ewald3_sigma_short(
 int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
 		complex double *target_sigmalr_y, // must be c->nelem_sc long
 		double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		const double eta, const complex double k,
 	        const double unitcell_volume /* with the corresponding lattice dimensionality */,
 		const LatticeDimensionality latdim,
@@ -182,11 +182,10 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
 		const cart3_t particle_shift
 		);
 		
-/// !!!!!!!!!!!!!!! ZDE JSEM SKONČIL !!!!!!!!!!!!!!!!!!!!!!.
-
+#ifdef EWALD_LEGACY // moved to ewald_legacy.c, not even everything implemented
 
 int ewald32_sigma0(complex double *result, double *err, // actually, this should be only alias for ewald3_sigma0
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k
 );
 
@@ -196,7 +195,7 @@ int ewald32_sigma0(complex double *result, double *err, // actually, this should
 int ewald32_sigma_long_shiftedpoints ( 
 		complex double *target_sigmalr_y, // must be c->nelem_sc long
 		double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k, double unitcell_area,
 		size_t npoints, const point2d *Kpoints_plus_beta,
 		point2d beta,
@@ -205,7 +204,7 @@ int ewald32_sigma_long_shiftedpoints (
 int ewald32_sigma_long_points_and_shift (
 		complex double *target_sigmalr_y, // must be c->nelem_sc long
 		double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k, double unitcell_area,
 		size_t npoints, const point2d *Kpoints,
 		point2d beta,
@@ -214,7 +213,7 @@ int ewald32_sigma_long_points_and_shift (
 int ewald32_sigma_long_shiftedpoints_rordered(//NI
 		complex double *target_sigmalr_y, // must be c->nelem_sc long
 		double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k, double unitcell_area,
 		const points2d_rordered_t *Kpoints_plus_beta_rordered,
 		point2d particle_shift
@@ -223,7 +222,7 @@ int ewald32_sigma_long_shiftedpoints_rordered(//NI
 int ewald32_sigma_short_shiftedpoints(
 		complex double *target_sigmasr_y, // must be c->nelem_sc long
 		double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c, // N.B. not too useful here
+		const qpms_ewald3_constants_t *c, // N.B. not too useful here
 		double eta, double k,
 		size_t npoints, const point2d *Rpoints_plus_particle_shift,
 		point2d beta,
@@ -232,7 +231,7 @@ int ewald32_sigma_short_shiftedpoints(
 int ewald32_sigma_short_points_and_shift(
 		complex double *target_sigmasr_y, // must be c->nelem_sc long
 		double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c, // N.B. not too useful here
+		const qpms_ewald3_constants_t *c, // N.B. not too useful here
 		double eta, double k,
 		size_t npoints, const point2d *Rpoints, 
 		point2d beta,
@@ -241,7 +240,7 @@ int ewald32_sigma_short_points_and_shift(
 int ewald32_sigma_short_points_rordered(//NI
 		complex double *target_sigmasr_y, // must be c->nelem_sc long
 		double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c, // N.B. not too useful here
+		const qpms_ewald3_constants_t *c, // N.B. not too useful here
 		double eta, double k,
 		const points2d_rordered_t *Rpoints_plus_particle_shift_rordered,
 		point2d particle_shift    // used only in the very end to multiply it by the phase
@@ -252,7 +251,7 @@ int ewald32_sigma_short_points_rordered(//NI
 int ewald31z_sigma_long_points_and_shift (
 		complex double *target_sigmalr_y, // must be c->nelem_sc long
 		double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k, double unitcell_area,
 		size_t npoints, const double *Kpoints,
 		double beta,
@@ -261,16 +260,17 @@ int ewald31z_sigma_long_points_and_shift (
 int ewald31z_sigma_short_points_and_shift(
 		complex double *target_sigmasr_y, // must be c->nelem_sc long
 		double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
-		const qpms_ewald32_constants_t *c, // N.B. not too useful here
+		const qpms_ewald3_constants_t *c, // N.B. not too useful here
 		double eta, double k,
 		size_t npoints, const double *Rpoints, 
 		double beta,
 		double particle_shift
 		);
 int ewald31z_sigma0(complex double *result, double *err, 
-		const qpms_ewald32_constants_t *c,
+		const qpms_ewald3_constants_t *c,
 		double eta, double k
 		); // exactly the same as ewald32_sigma0
 
+#endif // EWALD_LEGACY
 
 #endif //EWALD_H
diff --git a/qpms/ewald_legacy.c b/qpms/ewald_legacy.c
new file mode 100644
index 0000000..c056456
--- /dev/null
+++ b/qpms/ewald_legacy.c
@@ -0,0 +1,434 @@
+#include "ewald.h"
+#include <stdlib.h>
+#include "indexing.h"
+#include "kahansum.h"
+#include <assert.h>
+#include <string.h>
+#include <complex.h>
+#include "tiny_inlines.h"
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_sf_legendre.h>
+#include <gsl/gsl_sf_expint.h>
+
+// parameters for the quadrature of integral in (4.6)
+#ifndef INTEGRATION_WORKSPACE_LIMIT
+#define INTEGRATION_WORKSPACE_LIMIT 30000
+#endif
+
+#ifndef INTEGRATION_EPSABS
+#define INTEGRATION_EPSABS 1e-13
+#endif
+
+#ifndef INTEGRATION_EPSREL
+#define INTEGRATION_EPSREL 1e-13
+#endif
+
+#ifndef M_SQRTPI
+#define M_SQRTPI 1.7724538509055160272981674833411452
+#endif
+
+// sloppy implementation of factorial
+static inline double factorial(const int n) {
+  assert(n >= 0);
+  if (n < 0)
+    return 0; // should not happen in the functions below. (Therefore the assert above)
+  else if (n <= 20) {
+    double fac = 1;
+    for (int i = 1; i <= n; ++i)
+      fac *= i;
+    return fac;
+  }
+  else 
+    return tgamma(n + 1); // hope it's precise and that overflow does not happen
+}
+
+static inline complex double csq(complex double x) { return x * x; }
+static inline double sq(double x) { return x * x; }
+
+
+typedef enum {
+  EWALD32_CONSTANTS_ORIG, // As in [1, (4,5)], NOT USED right now.
+  EWALD32_CONSTANTS_AGNOSTIC /* Not depending on the spherical harmonic sign/normalisation
+                              * convention – the $e^{im\alpha_pq}$ term in [1,(4.5)] being
+                              * replaced by the respective $Y_n^m(\pi/2,\alpha)$ 
+                              * spherical harmonic. See notes/ewald.lyx.
+                              */
+
+} ewald3_constants_option;
+
+static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC;
+
+
+
+int ewald32_sigma0(complex double *result, double *err,
+    const qpms_ewald3_constants_t *c,
+    const double eta, const double k)
+{
+  return ewald3_sigma0(result, err, c, eta, k);
+}
+
+int ewald32_sigma_long_shiftedpoints (
+    complex double *target, // must be c->nelem_sc long
+    double *err,
+    const qpms_ewald3_constants_t *c,
+    const double eta, const double k, const double unitcell_area,
+    const size_t npoints, const point2d *Kpoints_plus_beta,
+    const point2d beta, // not needed
+    const point2d particle_shift // target - src
+    ) 
+{
+  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_l_t lMax = c->lMax;
+  
+  // Manual init of the ewald summation targets
+  complex double *target_c = calloc(nelem_sc, sizeof(complex double));
+  memset(target, 0, nelem_sc * sizeof(complex double));
+  double *err_c = NULL;
+  if (err) {
+    err_c = calloc(nelem_sc, sizeof(double));
+    memset(err, 0, nelem_sc * sizeof(double));
+  }
+
+  const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC)
+  assert(commonfac > 0);
+
+  // space for Gamma_pq[j]'s
+  qpms_csf_result Gamma_pq[lMax/2+1];
+
+  // CHOOSE POINT BEGIN
+  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
+    const point2d beta_pq = Kpoints_plus_beta[i];
+    const point2d K_pq = {beta_pq.x - beta.x, beta_pq.y - beta.y};
+    const double rbeta_pq = cart2norm(beta_pq);
+  // CHOOSE POINT END
+
+    const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!!
+    const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT
+
+    // R-DEPENDENT BEGIN
+    const complex double gamma_pq = lilgamma(rbeta_pq/k);
+    const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
+    for(qpms_l_t j = 0; j <= lMax/2; ++j) {
+      int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
+      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
+      if(creal(z) < 0) 
+        Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
+      if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
+    }
+    // R-DEPENDENT END
+    
+    // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
+    // and just fetched for each n, m pair
+    for(qpms_l_t n = 0; n <= lMax; ++n)
+      for(qpms_m_t m = -n; m <= n; ++m) {
+        if((m+n) % 2 != 0) // odd coefficients are zero.
+          continue;
+        const qpms_y_t y = qpms_mn2y_sc(m, n);
+        const complex double e_imalpha_pq = cexp(I*m*arg_pq);
+        complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
+        double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
+        assert((n-abs(m))/2 == c->s1_jMaxes[y]);
+        for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
+          complex double summand = pow(rbeta_pq/k, n-2*j) 
+            * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
+            * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
+            * c->s1_constfacs[y][j];
+          if(err) {
+            // FIXME include also other errors than Gamma_pq's relative error
+             kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand));
+          }
+          summand *= Gamma_pq[j].val; // GGG
+          ckahanadd(&jsum, &jsum_c, summand);
+        }
+        jsum *= phasefac; // PFC
+        ckahanadd(target + y, target_c + y, jsum);
+        if(err) kahanadd(err + y, err_c + y, jsum_err);
+      }
+  } // END POINT LOOP
+ 
+  free(err_c);
+  free(target_c);
+  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
+    target[y] *= commonfac;
+  if(err)
+    for(qpms_y_t y = 0; y < nelem_sc; ++y)
+      err[y] *= commonfac;
+  return 0;
+}
+
+int ewald32_sigma_long_points_and_shift (
+    complex double *target, // must be c->nelem_sc long
+    double *err,
+    const qpms_ewald3_constants_t *c,
+    const double eta, const double k, const double unitcell_area,
+    const size_t npoints, const point2d *Kpoints,
+    const point2d beta, 
+    const point2d particle_shift // target - src
+    ) 
+{
+  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_l_t lMax = c->lMax;
+  
+  // Manual init of the ewald summation targets
+  complex double *target_c = calloc(nelem_sc, sizeof(complex double));
+  memset(target, 0, nelem_sc * sizeof(complex double));
+  double *err_c = NULL;
+  if (err) {
+    err_c = calloc(nelem_sc, sizeof(double));
+    memset(err, 0, nelem_sc * sizeof(double));
+  }
+
+  const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC)
+  assert(commonfac > 0);
+
+  // space for Gamma_pq[j]'s
+  qpms_csf_result Gamma_pq[lMax/2+1];
+
+  // CHOOSE POINT BEGIN
+  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
+    // Only these following two lines differ from ewald32_sigma_long_points_and_shift()!!!  WTFCOMMENT?!
+    const point2d K_pq = Kpoints[i]; 
+    const point2d beta_pq = {K_pq.x + beta.x, K_pq.y + beta.y};
+    const double rbeta_pq = cart2norm(beta_pq);
+  // CHOOSE POINT END
+
+    const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!!
+    const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT
+
+    // R-DEPENDENT BEGIN
+    const complex double gamma_pq = lilgamma(rbeta_pq/k);
+    const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné
+    for(qpms_l_t j = 0; j <= lMax/2; ++j) {
+      int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j);
+      // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above
+      if(creal(z) < 0) 
+        Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above
+      if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
+    }
+    // R-DEPENDENT END
+    
+    // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
+    // and just fetched for each n, m pair
+    for(qpms_l_t n = 0; n <= lMax; ++n)
+      for(qpms_m_t m = -n; m <= n; ++m) {
+        if((m+n) % 2 != 0) // odd coefficients are zero.
+          continue;
+        const qpms_y_t y = qpms_mn2y_sc(m, n);
+        const complex double e_imalpha_pq = cexp(I*m*arg_pq);
+        complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
+        double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
+        assert((n-abs(m))/2 == c->s1_jMaxes[y]);
+        for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
+          complex double summand = pow(rbeta_pq/k, n-2*j) 
+            * e_imalpha_pq  * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
+            * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
+            * c->s1_constfacs[y][j];
+          if(err) {
+            // FIXME include also other errors than Gamma_pq's relative error
+             kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand));
+          }
+          summand *= Gamma_pq[j].val; // GGG
+          ckahanadd(&jsum, &jsum_c, summand);
+        }
+        jsum *= phasefac; // PFC
+        ckahanadd(target + y, target_c + y, jsum);
+        if(err) kahanadd(err + y, err_c + y, jsum_err);
+      }
+  } // END POINT LOOP
+ 
+  free(err_c);
+  free(target_c);
+  for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
+    target[y] *= commonfac;
+  if(err)
+    for(qpms_y_t y = 0; y < nelem_sc; ++y)
+      err[y] *= commonfac;
+  return 0;
+}
+
+struct sigma2_integrand_params {
+  int n;
+  double k, R;
+};
+
+static double sigma2_integrand(double ksi, void *params) {
+  struct sigma2_integrand_params *p = (struct sigma2_integrand_params *) params;
+  return exp(-sq(p->R * ksi) + sq(p->k / ksi / 2)) * pow(ksi, 2*p->n);
+}
+
+static int ewald32_sr_integral(double r, double k, int n, double eta,
+    double *result, double *err, gsl_integration_workspace *workspace)
+{
+  struct sigma2_integrand_params p;
+
+
+  const double R0 = r; // Maybe could be chosen otherwise, but fuck it for now.
+  p.n = n;
+  eta *= R0;
+  p.k = k * R0;
+  p.R = r / R0;  // i.e. p.R = 1
+
+  gsl_function F;
+  F.function = sigma2_integrand;
+  F.params = &p;
+
+  int retval = gsl_integration_qagiu(&F, eta, INTEGRATION_EPSABS, 
+      INTEGRATION_EPSREL, INTEGRATION_WORKSPACE_LIMIT,
+      workspace, result, err);
+  double normfac = pow(R0, -2*p.n - 1);
+  *result *= normfac;
+  *err *= normfac;
+  return retval;
+}
+
+
+int ewald32_sigma_short_shiftedpoints(
+    complex double *target, // must be c->nelem_sc long
+    double *err,
+    const qpms_ewald3_constants_t *c, // N.B. not too useful here
+    const double eta, const double k,
+    const size_t npoints, const point2d *Rpoints_plus_particle_shift,
+    const point2d beta,
+    const point2d particle_shift           // used only in the very end to multiply it by the phase
+    )
+{
+  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_l_t lMax = c->lMax;
+  gsl_integration_workspace *workspace = 
+    gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
+  
+  // Manual init of the ewald summation targets
+  complex double * const target_c = calloc(nelem_sc, sizeof(complex double));
+  memset(target, 0, nelem_sc * sizeof(complex double));
+  double *err_c = NULL;
+  if (err) {
+    err_c = calloc(nelem_sc, sizeof(double));
+    memset(err, 0, nelem_sc * sizeof(double));
+  }
+  
+
+// CHOOSE POINT BEGIN
+  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
+    const point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
+    const double r_pq_shifted = cart2norm(Rpq_shifted);
+// CHOOSE POINT END
+
+    const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT
+    const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
+    
+    for(qpms_l_t n = 0; n <= lMax; ++n) {
+      const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
+      const double R_pq_pown = pow(r_pq_shifted, n);
+      // TODO the integral here
+      double intres, interr;
+      int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta,
+          &intres, &interr, workspace);
+      if (retval) abort();
+      for (qpms_m_t m = -n; m <= n; ++m){
+        if((m+n) % 2 != 0) // odd coefficients are zero.
+          continue; // nothing needed, already done by memset
+        const complex double e_imf = cexp(I*m*Rpq_shifted_arg);
+        const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
+        const qpms_y_t y = qpms_mn2y_sc(m,n);
+        if(err)
+          kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
+              * interr)); // TODO include also other errors
+        ckahanadd(target + y, target_c + y,
+            prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m));
+      }
+
+    }
+  }
+
+  gsl_integration_workspace_free(workspace);
+  if(err) free(err_c);
+  free(target_c);
+  return 0;
+}
+
+int ewald32_sigma_short_points_and_shift(
+    complex double *target, // must be c->nelem_sc long
+    double *err,
+    const qpms_ewald3_constants_t *c, // N.B. not too useful here
+    const double eta, const double k,
+    const size_t npoints, const point2d *Rpoints,
+    const point2d beta,
+    const point2d particle_shift           // used only in the very end to multiply it by the phase
+    )
+{
+  const qpms_y_t nelem_sc = c->nelem_sc;
+  const qpms_l_t lMax = c->lMax;
+  gsl_integration_workspace *workspace = 
+    gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
+  
+  // Manual init of the ewald summation targets
+  complex double * const target_c = calloc(nelem_sc, sizeof(complex double));
+  memset(target, 0, nelem_sc * sizeof(complex double));
+  double *err_c = NULL;
+  if (err) {
+    err_c = calloc(nelem_sc, sizeof(double));
+    memset(err, 0, nelem_sc * sizeof(double));
+  }
+  
+
+// CHOOSE POINT BEGIN
+  for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP
+    //const point2d Rpq_shifted = Rpoints_plus_particle_shift[i];
+    const point2d Rpq_shifted = cart2_add(Rpoints[i], cart2_scale(-1,particle_shift)); // CHECKSIGN!!!!
+    const double r_pq_shifted = cart2norm(Rpq_shifted);
+// CHOOSE POINT END
+
+    const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT
+    const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
+    
+    for(qpms_l_t n = 0; n <= lMax; ++n) {
+      const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
+      const double R_pq_pown = pow(r_pq_shifted, n);
+      // TODO the integral here
+      double intres, interr;
+      int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta,
+          &intres, &interr, workspace);
+      if (retval) abort();
+      for (qpms_m_t m = -n; m <= n; ++m){
+        if((m+n) % 2 != 0) // odd coefficients are zero.
+          continue; // nothing needed, already done by memset
+        const complex double e_imf = cexp(I*m*Rpq_shifted_arg);
+        const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))];
+        const qpms_y_t y = qpms_mn2y_sc(m,n);
+        if(err)
+          kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
+              * interr)); // TODO include also other errors
+        ckahanadd(target + y, target_c + y,
+            prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m));
+      }
+
+    }
+  }
+
+  gsl_integration_workspace_free(workspace);
+  if(err) free(err_c);
+  free(target_c);
+  return 0;
+}
+
+
+ 
+#if 0
+
+int ewald32_sigma_long_shiftedpoints_rordered(
+    complex double *target_sigmalr_y, // must be c->nelem_sc long
+    const qpms_ewald3_constants_t *c,
+    double eta, double k, double unitcell_area,
+    const points2d_rordered_t *Kpoints_plus_beta_rordered,
+    point2d particle_shift
+    );
+int ewald32_sigma_short_points_rordered(
+    complex double *target_sigmasr_y, // must be c->nelem_sc long
+    const qpms_ewald3_constants_t *c, // N.B. not too useful here
+    double eta, double k,
+    const points2d_rordered_t *Rpoints_plus_particle_shift_rordered,
+    point2d particle_shift    // used only in the very end to multiply it by the phase
+    );
+
+#endif
diff --git a/tests/ewalds.c b/tests/ewalds.c
index 5623403..32af265 100644
--- a/tests/ewalds.c
+++ b/tests/ewalds.c
@@ -185,7 +185,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 
-  qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase);
+  qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 
   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
       8*DBL_EPSILON, 8*DBL_EPSILON);
@@ -252,7 +252,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   }
 
   points2d_rordered_free(Kpoints_plus_beta);
-  qpms_ewald32_constants_free(c);
+  qpms_ewald3_constants_free(c);
   triangular_lattice_gen_free(Klg);
   triangular_lattice_gen_free(Rlg);
   ++ewaldtest_counter;
diff --git a/tests/ewaldshift.c b/tests/ewaldshift.c
index 0d00e87..3bb2984 100644
--- a/tests/ewaldshift.c
+++ b/tests/ewaldshift.c
@@ -266,7 +266,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 
-  qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase);
+  qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 
   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
       8*DBL_EPSILON, 8*DBL_EPSILON);
@@ -332,7 +332,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   }
 
   points2d_rordered_free(Kpoints_plus_beta);
-  qpms_ewald32_constants_free(c);
+  qpms_ewald3_constants_free(c);
   triangular_lattice_gen_free(Klg);
   triangular_lattice_gen_free(Rlg);
   ++ewaldtest_counter;
diff --git a/tests/ewaldshift2.c b/tests/ewaldshift2.c
index 77bd250..f9460b7 100644
--- a/tests/ewaldshift2.c
+++ b/tests/ewaldshift2.c
@@ -295,7 +295,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 
-  qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase);
+  qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 
   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
       8*DBL_EPSILON, 8*DBL_EPSILON);
@@ -361,7 +361,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   }
 
   points2d_rordered_free(Kpoints_plus_beta);
-  qpms_ewald32_constants_free(c);
+  qpms_ewald3_constants_free(c);
   triangular_lattice_gen_free(Klg);
   triangular_lattice_gen_free(Rlg);
   ++ewaldtest_counter;
diff --git a/tests/ewaldshift_3g.c b/tests/ewaldshift_3g.c
index 83b51e9..7b493b9 100644
--- a/tests/ewaldshift_3g.c
+++ b/tests/ewaldshift_3g.c
@@ -307,7 +307,7 @@ ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) {
   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 
-  qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase);
+  qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 
   //points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
   //    8*DBL_EPSILON, 8*DBL_EPSILON);
@@ -395,7 +395,7 @@ ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) {
 
 
   //points2d_rordered_free(Kpoints_plus_beta);
-  qpms_ewald32_constants_free(c);
+  qpms_ewald3_constants_free(c);
   //triangular_lattice_gen_free(Klg);
   //triangular_lattice_gen_free(Rlg);
   ++ewaldtest_counter;
diff --git a/tests/ewaldshift_ng.c b/tests/ewaldshift_ng.c
index e3d0991..c5172ad 100644
--- a/tests/ewaldshift_ng.c
+++ b/tests/ewaldshift_ng.c
@@ -295,7 +295,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 
-  qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase);
+  qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 
   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
       8*DBL_EPSILON, 8*DBL_EPSILON);
@@ -361,7 +361,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
   }
 
   points2d_rordered_free(Kpoints_plus_beta);
-  qpms_ewald32_constants_free(c);
+  qpms_ewald3_constants_free(c);
   triangular_lattice_gen_free(Klg);
   triangular_lattice_gen_free(Rlg);
   ++ewaldtest_counter;