diff --git a/qpms/ewald.c b/qpms/ewald.c
index bc2fd4a..3fec709 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -28,6 +28,7 @@
 #endif
 
 
+
 // sloppy implementation of factorial
 static inline double factorial(const int n) {
   assert(n >= 0);
@@ -133,6 +134,22 @@ void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) {
 
 
 
+int ewald32_sigma0(complex double *result, double *err,
+    const qpms_ewald32_constants_t *c,
+    const double eta, const double k)
+{
+  qpms_csf_result gam;
+  int retval = complex_gamma_inc_e(-0.5, -k/sq(2*eta), &gam);
+  if (0 != retval)
+    abort();
+  *result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI;
+  if(err) 
+    *err = gam.err * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI;
+  return 0;
+}
+
+
+
 int ewald32_sigma_long_shiftedpoints (
     complex double *target, // must be c->nelem long
     double *err,
diff --git a/qpms/ewald.h b/qpms/ewald.h
index 8497e92..2fb42af 100644
--- a/qpms/ewald.h
+++ b/qpms/ewald.h
@@ -78,6 +78,12 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
 
 #include "lattices.h"
 
+
+int ewald32_sigma0(complex double *result, double *err,
+		const qpms_ewald32_constants_t *c,
+		double eta, double k
+);
+
 // TODO make "compressed versions" where the (m+n)-odd terms (which are zero)
 // are not included.
 
@@ -114,6 +120,7 @@ int ewald32_sigma_short_shiftedpoints(
 		const qpms_ewald32_constants_t *c, // N.B. not too useful here
 		double eta, double k,
 		size_t npoints, const point2d *Rpoints_plus_particle_shift,
+		point2d beta,
 		point2d particle_shift           // used only in the very end to multiply it by the phase
 		);
 int ewald32_sigma_short_points_and_shift(//NI