diff --git a/qpms/ewald.c b/qpms/ewald.c
index c992c50..7cf52ca 100644
--- a/qpms/ewald.c
+++ b/qpms/ewald.c
@@ -257,7 +257,8 @@ int ewald3_21_xy_sigma_long (
 #endif
   // recycleable values if rbeta_pq stays the same:
   complex double gamma_pq;
-  complex double z;
+  complex double z; // I?*κ*γ*particle_shift.z
+  complex double x; // κ*γ/(2*η)
   complex double factor1d = 1; // the "additional" factor for the 1D case (then it is not 1)
   // space for Gamma_pq[j]'s
   qpms_csf_result Gamma_pq[lMax/2+1];
@@ -291,17 +292,28 @@ int ewald3_21_xy_sigma_long (
 
 
     // R-DEPENDENT BEGIN
+    //void ewald3_2_sigma_long_Delta(complex double *target, int maxn, complex double x, complex double z) {
     if (new_rbeta_pq) {
       gamma_pq = clilgamma(rbeta_pq/k);
-      z = csq(gamma_pq*k/(2*eta)); 
-      for(qpms_l_t j = 0; j <= lMax/2; ++j) {
-        int retval = complex_gamma_inc_e(0.5-j, z,
-          // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
-          QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j);
-        QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW);
+      complex double x = gamma_pq*k/(2*eta);
+      complex double x2 = x*x;
+      if(particle_shift.z == 0) {
+        for(qpms_l_t j = 0; j <= lMax/2; ++j) {
+          int retval = complex_gamma_inc_e(0.5-j, x2,
+            // we take the branch which is principal for the Re x < 0, Im x < 0 quadrant, cf. [Linton, p. 642 in the middle]
+            QPMS_LIKELY(creal(x2) < 0) && !signbit(cimag(x2)) ? -1 : 0, Gamma_pq+j);
+          QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW);
+        } 
+        if (latdim & LAT1D)
+          factor1d =  M_SQRT1_2 * .5 * k * gamma_pq;
+      } else if (latdim & LAT2D) {
+        QPMS_UNTESTED;
+        // TODO check the branches/phases!
+        complex double z = I * k * gamma_pq * particle_shift.z;
+        ewald3_2_sigma_long_Delta(Gamma_pq, lMax/2, x, z);
+      } else {
+        QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented");
       }
-      if (latdim & LAT1D)
-        factor1d =  M_SQRT1_2 * .5 * k * gamma_pq;
     }
     // R-DEPENDENT END
     
diff --git a/qpms/ewald.h b/qpms/ewald.h
index adde9d7..722e838 100644
--- a/qpms/ewald.h
+++ b/qpms/ewald.h
@@ -210,6 +210,12 @@ int hyperg_2F2_series(double a, double b, double c, double d, double x,
 int ewald32_sr_integral(double r, double k, double n, double eta, double *result, double *err, gsl_integration_workspace *workspace);
 #endif
 
+/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
+/** \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f]
+ *
+ * NOTE: The error part is not yet implemented, NANs are returned instead!
+ */
+void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z);
 
 // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space
 
diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c
index ddf0bc1..69bce2e 100644
--- a/qpms/ewaldsf.c
+++ b/qpms/ewaldsf.c
@@ -12,6 +12,7 @@
 #include <gsl/gsl_errno.h>
 #include <float.h>
 #include <stdbool.h>
+#include <Faddeeva.h>
 
 #ifndef COMPLEXPART_REL_ZERO_LIMIT
 #define COMPLEXPART_REL_ZERO_LIMIT 1e-14
@@ -289,5 +290,28 @@ int hyperg_2F2_series(const double a, const double b, const double c, const doub
   return GSL_SUCCESS;
 }
 
+// The Delta_n factor from [Kambe II], Appendix 3
+// \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f]
+void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z) {
+  complex double expfac = cexp(-x + 0.25 * z*z / x);
+  complex double sqrtx = csqrt(x); // TODO check carefully, which branch is needed
+  // These are used in the first two recurrences
+  complex double w_plus  = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0);
+  complex double w_minus = Faddeeva_w(-z/(2*sqrtx) + I*sqrtx, 0);
+  QPMS_ASSERT(maxn >= 0); 
+  if (maxn >= 0) {
+    target[0].val = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus);
+    target[0].err = NAN; //TODO
+  }
+  if (maxn >= 1) {
+    target[1].val = -I / z * M_SQRTPI * expfac * (w_minus - w_plus);
+    target[1].err = NAN; //TODO
+  }
+  for(int n = 1; n < maxn; ++n) { // The rest via recurrence
+    // TODO The cpow(x, 0.5 - n) might perhaps better be replaced with a recurrently computed variant
+    target[n+1].val = (4 / (z*z)) * ((0.5 - n) * target[n].val + target[n-1].val - cpow(x, 0.5 - n) * expfac);
+    target[n+1].err = NAN; //TODO
+  }
+}