From 67f93e461cb352318f327d3d89567f7b908b2a9a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 23 Sep 2019 08:59:18 +0300
Subject: [PATCH] Ewald translation calculator function using base specs.

Former-commit-id: d49880caa280cc089fa688d0cdeccb95db7cb64b
---
 qpms/translations.c | 55 +++++++++++++++++++++++++++++++++++++++++++++
 qpms/translations.h | 14 ++++++++++++
 2 files changed, 69 insertions(+)

diff --git a/qpms/translations.c b/qpms/translations.c
index 2832cdc..20a58d8 100644
--- a/qpms/translations.c
+++ b/qpms/translations.c
@@ -705,6 +705,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
   assert(c->normalisation == destspec->norm && c->normalisation == srcspec->norm);
   assert(c->lMax >= destspec->lMax && c->lMax >= srcspec->lMax);
   assert(destspec->lMax_L < 0 && srcspec->lMax_L < 0);
+  // TODO don't use c->lMax etc. if both destspec->lMax and srcspec->lMax are smaller
   complex double A[c->nelem][c->nelem];
   complex double B[c->nelem][c->nelem];
   qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays(c,
@@ -727,6 +728,60 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
   return retval;
 }
 
+qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
+    complex double *target, double *err,
+    /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
+    const qpms_vswf_set_spec_t *destspec, size_t deststride,
+    /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
+    const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
+    const double eta, const complex double k,
+    cart2_t b1, cart2_t b2,
+    const cart2_t beta,
+    const cart2_t particle_shift,
+    double maxR, double maxK
+    )
+{
+  TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c->normalisation);
+  QPMS_ENSURE(c->normalisation == destspec->norm && c->normalisation == srcspec->norm,
+      "The normalisation conventions must be the same");
+  assert(c->lMax >= destspec->lMax && c->lMax >= srcspec->lMax);
+  assert(destspec->lMax_L < 0 && srcspec->lMax_L < 0);
+  // TODO don't use c->lMax etc. if both destspec->lMax and srcspec->lMax are smaller
+  const ptrdiff_t ldAB = c->nelem;
+  complex double *A, *B;
+  double *Aerr = NULL, *Berr = NULL;
+  QPMS_CRASHING_MALLOC(A, c->nelem*c->nelem*sizeof(complex double));
+  QPMS_CRASHING_MALLOC(B, c->nelem*c->nelem*sizeof(complex double));
+  if(err) {
+    QPMS_CRASHING_MALLOC(Aerr, c->nelem*c->nelem*sizeof(double));
+    QPMS_CRASHING_MALLOC(Berr, c->nelem*c->nelem*sizeof(double));
+  }
+  qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays_e32(c,
+      A, Aerr, B, Berr, ldAB, 1, 
+      eta, k, b1, b2, beta, particle_shift, maxR, maxK);
+  for (size_t desti = 0; desti < destspec->n; ++desti) {
+    qpms_y_t desty; qpms_vswf_type_t destt;
+    if(QPMS_SUCCESS != qpms_uvswfi2ty(destspec->ilist[desti], &destt, &desty))
+        qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,
+          "Invalid u. vswf index %llx.", destspec->ilist[desti]);
+    for (size_t srci = 0; srci < srcspec->n; ++srci){
+      qpms_y_t srcy; qpms_vswf_type_t srct;
+      if(QPMS_SUCCESS != qpms_uvswfi2ty(srcspec->ilist[srci], &srct, &srcy))
+          qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,
+            "Invalid u. vswf index %llx.", srcspec->ilist[srci]);
+      target[srci * srcstride + desti * deststride]
+        = (srct == destt) ? A[ldAB*desty + srcy] : B[ldAB*desty + srcy];
+      if(err) err[srci * srcstride + desti * deststride]
+        = (srct == destt) ? Aerr[ldAB*desty + srcy] : Berr[ldAB*desty + srcy];
+    }
+  }
+  free(A); free(B);
+  if (err) { free(Aerr); free(Berr); }
+  return retval;
+}
+
+
+
 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( 
     const qpms_trans_calculator *c,
     complex double *target,
diff --git a/qpms/translations.h b/qpms/translations.h
index ef8cb48..0921dad 100644
--- a/qpms/translations.h
+++ b/qpms/translations.h
@@ -166,6 +166,20 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
 		double maxR, double maxK
 		);
 
+// Convenience functions using VSWF base specs
+qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
+		complex double *target, double *err,
+		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
+		const qpms_vswf_set_spec_t *destspec, size_t deststride,
+		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
+		const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
+		const double eta, const complex double k,
+		cart2_t b1, cart2_t b2,
+		const cart2_t beta,
+		const cart2_t particle_shift,
+		double maxR, double maxK
+		);
+
 #endif //LATTICESUMS32
 
 #ifdef LATTICESUMS31