diff --git a/qpms/ewald.c b/qpms/ewald.c
index 02b8ffd..603fe44 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -124,6 +124,7 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
                                      : isq(sz_n_lMax);
   }
   c->s1_constfacs_1Dz_base = malloc(s1_constfacs_1Dz_sz * sizeof(complex double));
+  c->s1_constfacs_1Dz = malloc((lMax+1)*sizeof(complex double *));
 
   size_t s1_constfacs_1Dz_sz_cumsum = 0;
   for (qpms_l_t n = 0; n <= lMax; ++n) {
@@ -131,7 +132,7 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
     for (qpms_l_t j = 0; j <= n/2; ++j) {
       switch(type) {
         case EWALD32_CONSTANTS_AGNOSTIC:
-          c->s1_constfacs[n][j] = -ipow(n+1) * min1pow(j) * factorial(n)
+          c->s1_constfacs_1Dz[n][j] = -ipow(n+1) * min1pow(j) * factorial(n)
             / (factorial(j) * pow(2, 2*j) * factorial(n - 2*j));
           break;
         default:
@@ -145,6 +146,8 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
 
   c->legendre_csphase = csphase;
   c->legendre0 = malloc(gsl_sf_legendre_array_n(lMax) * sizeof(double));
+  c->legendre_plus1 = malloc(gsl_sf_legendre_array_n(lMax) * sizeof(double));
+  c->legendre_minus1 = malloc(gsl_sf_legendre_array_n(lMax) * sizeof(double));
   // N.B. here I use the GSL_SF_LEGENRE_NONE, in order to be consistent with translations.c
   c->legendre_normconv = GSL_SF_LEGENDRE_NONE;
   // Moreover, using this approach (i.e. gsl) takes about 64kB extra memory
@@ -159,9 +162,12 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
 
 void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) {
   free(c->legendre0);
+  free(c->legendre_plus1);
+  free(c->legendre_minus1);
   free(c->s1_constfacs);
   free(c->s1_constfacs_base);
   free(c->s1_constfacs_1Dz_base);
+  free(c->s1_constfacs_1Dz);
   free(c->s1_jMaxes);
   free(c);
 }
@@ -378,7 +384,7 @@ int ewald3_21_xy_sigma_long (
     const double eta, const double k,
     const double unitcell_volume /* with the corresponding lattice dimensionality */,
     const LatticeDimensionality latdim,
-    PGenSph *pgen_K, const bool pgen_generates_shifted_points
+    PGen *pgen_K, const bool pgen_generates_shifted_points
     /* If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift,
      * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
      * and adds beta to the generated points before calculations.
@@ -417,7 +423,8 @@ int ewald3_21_xy_sigma_long (
   qpms_csf_result Gamma_pq[lMax/2+1];
 
   // CHOOSE POINT BEGIN
-  while ((pgen_retdata = PGenSph_next(pgen_K)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
+  // TODO mayby PGen_next_sph is not the best coordinate system choice here
+  while ((pgen_retdata = PGen_next_sph(pgen_K)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
     cart3_t K_pq_cart;
     sph_t beta_pq_sph;
     if (pgen_generates_shifted_points) {
@@ -511,7 +518,7 @@ int ewald3_1_z_sigma_long (
     const double eta, const double k,
     const double unitcell_volume /* length (periodicity) in this case */,
     const LatticeDimensionality latdim,
-    PGenSph *pgen_K, const bool pgen_generates_shifted_points
+    PGen *pgen_K, const bool pgen_generates_shifted_points
     /* If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift,
      * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
      * and adds beta to the generated points before calculations.
@@ -546,7 +553,8 @@ int ewald3_1_z_sigma_long (
 
   PGenSphReturnData pgen_retdata;
   // CHOOSE POINT BEGIN
-  while ((pgen_retdata = PGenSph_next(pgen_K)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
+  // TODO FIXME USE PGen_next_z
+  while ((pgen_retdata = PGen_next_sph(pgen_K)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
     assert(pgen_retdata.flags & PGEN_AT_Z);
     double K_z, beta_mu_z;
     if (pgen_generates_shifted_points) {
@@ -626,7 +634,7 @@ int ewald3_sigma_long (
     const double eta, const double k,
     const double unitcell_volume /* with the corresponding lattice dimensionality */,
     const LatticeDimensionality latdim,
-    PGenSph *pgen_K, const bool pgen_generates_shifted_points
+    PGen *pgen_K, const bool pgen_generates_shifted_points
     /* If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift,
      * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
      * and adds beta to the generated points before calculations.
@@ -817,7 +825,7 @@ int ewald3_sigma_short(
                 const qpms_ewald32_constants_t *c,
                 const double eta, const double k,
                 const LatticeDimensionality latdim, // apart from asserts and possible optimisations ignored, as the SR formula stays the same
-                PGenSph *pgen_R, const bool pgen_generates_shifted_points
+                PGen *pgen_R, const bool pgen_generates_shifted_points
                 /* If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift,
                  * so the function assumes that the generated points correspond to the unshifted Bravais lattice,
                  * and adds particle_shift to the generated points before calculations.
@@ -853,7 +861,8 @@ int ewald3_sigma_short(
   double intres[lMax+1], interr[lMax+1];
 
 // CHOOSE POINT BEGIN
-  while ((pgen_retdata = PGenSph_next(pgen_R)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
+// TODO check whether _next_sph is the optimal coordinate system choice here
+  while ((pgen_retdata = PGen_next_sph(pgen_R)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
 // CHOOSE POINT END
     cart3_t Rpq_shifted_cart; // I will need both sph and cart representations anyway...
     sph_t Rpq_shifted_sph;
diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c
index 539bf02..2c308f9 100644
--- a/qpms/ewaldsf.c
+++ b/qpms/ewaldsf.c
@@ -113,8 +113,8 @@ int complex_expint_n_e(int n, complex double x, qpms_csf_result *result) {
   } else {
     int retval = complex_gamma_inc_e(-n+1, x, result);
     complex double f = cpow(x, 2*n-2);
-    retval.val *= f;
-    retval.err *= cabs(f);
+    result->val *= f;
+    result->err *= cabs(f);
     return retval;
   }
 }
diff --git a/qpms/lattices.h b/qpms/lattices.h
index 30aae93..578655f 100644
--- a/qpms/lattices.h
+++ b/qpms/lattices.h
@@ -159,6 +159,12 @@ static inline void PGen_destroy(PGen *g) {
 	assert(g->stateData == NULL); // this should be done by the destructor
 }
 
+static inline PGenZReturnData PGen_next_z(PGen *g) {
+	if (g->c->next_z)
+		return g->c->next_z(g);
+	else abort();
+}
+
 static inline PGenSphReturnData PGen_next_sph(PGen *g) {
 	// TODO maybe some asserts around here
 	if (g->c->next_sph) 
diff --git a/qpms/vectors.h b/qpms/vectors.h
index 99f63ea..3557c4c 100644
--- a/qpms/vectors.h
+++ b/qpms/vectors.h
@@ -8,6 +8,9 @@
 
 //static inline double vectors_h_sq(double x) {return x*x;}
 
+static const cart2_t CART2_ZERO = {0, 0};
+static const cart3_t CART3_ZERO = {0, 0, 0};
+
 static inline cart2_t cart2_add(const cart2_t a, const cart2_t b) {
 	cart2_t res = {a.x+b.x, a.y+b.y};
 	return res;
@@ -66,6 +69,11 @@ static inline cart3_t cart22cart3xy(const cart2_t a) {
 	return c;
 }
 
+static inline cart2_t cart3xy2cart2(const cart3_t a) {
+	cart2_t c = {a.x, a.y};
+	return c;
+}
+
 static inline double cart3norm(const cart3_t v) {
 	return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
 }