some tests and calculate the translation matrix separately; still not working

Maybe the problem is with the x-y flips?


Former-commit-id: 6a1ffe4c9697c7ece90a359211ee4a82b2dbb538

qpms/qpms_c.pyx

@@ -1463,6 +1463,15 @@ cdef class ScatteringSystem:
         qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k)
         return target
 
+    def translation_matrix_full(self, double k):
+        cdef size_t flen = self.s[0].fecv_size
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (flen,flen),dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        qpms_scatsys_build_translation_matrix_full(&target_view[0][0], self.s, k)
+        return target
+
+
 
 def tlm2uvswfi(t, l, m):
     ''' TODO doc

qpms/qpms_cdefs.pxd

@@ -287,6 +287,8 @@ cdef extern from "scatsystem.h":
             const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
     cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target,
             const qpms_scatsys_t *ss, double k)
+    cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target,
+            const qpms_scatsys_t *ss, double k)
     cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target,
             const qpms_scatsys_t *ss, qpms_iri_t iri, double k)
 

qpms/scatsystem.c

@@ -1032,6 +1032,43 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
   return target_full; 
 }
 
+complex double *qpms_scatsys_build_translation_matrix_full(
+    /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
+    complex double *target,
+    const qpms_scatsys_t *ss,
+    double k ///< Wave number to use in the translation matrix.
+    )
+{
+  const size_t full_len = ss->fecv_size;
+  if(!target)
+    QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
+  memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary?
+  { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC)
+    size_t fullvec_offsetR = 0;
+    for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) {
+      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+      const cart3_t posR = ss->p[piR].pos;
+      size_t fullvec_offsetC = 0;
+      for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) {
+        if(piC == piR) continue; // The diagonal will be dealt with later.
+        const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+        const cart3_t posC = ss->p[piC].pos;
+        QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
+              target + fullvec_offsetR*full_len + fullvec_offsetC,
+              bspecR, full_len, bspecC, 1,
+              k, posR, posC));      
+        fullvec_offsetC += bspecC->n;
+      }
+      assert(fullvec_offsetC = full_len);
+      fullvec_offsetR += bspecR->n;
+    }
+    assert(fullvec_offsetR == full_len);
+  }
+  
+  return target;
+}
+
+
 complex double *qpms_scatsys_build_modeproblem_matrix_full(
     /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
     complex double *target,

qpms/scatsystem.h

@@ -450,6 +450,12 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
 		const complex double *orig_packed, const qpms_scatsys_t *ss,
 		qpms_iri_t iri, bool add);
 
+complex double *qpms_scatsys_build_translation_matrix_full(
+		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
+		complex double *target,
+		const qpms_scatsys_t *ss,
+		double k ///< Wave number to use in the translation matrix.
+		);
 complex double *qpms_scatsys_build_modeproblem_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,

tests/bspectransl.c

@@ -0,0 +1,42 @@
+//  c99 -g -I.. bspectransl.c ../qpms/translations.c ../qpms/legendre.c ../qpms/vswf.c ../qpms/gaunt.c ../qpms/error.c -lm -lgsl -lblas
+#include <qpms/translations.h>
+#include <qpms/vswf.h>
+#include <complex.h>
+#include <stdio.h>
+
+int main() {
+  cart3_t pos1={0,1,2}, pos2={3,5,2};
+  qpms_vswf_set_spec_t 
+    *bspec1 = qpms_vswf_set_spec_from_lMax(1, QPMS_NORMALISATION_POWER),
+    *bspec2 = qpms_vswf_set_spec_from_lMax(2, QPMS_NORMALISATION_POWER);
+
+  qpms_trans_calculator *c = qpms_trans_calculator_init(3, QPMS_NORMALISATION_POWER);
+
+  complex double s_2_1[bspec2->n][bspec1->n];
+  complex double s_1_2[bspec1->n][bspec2->n];
+  const double k = 1.8;
+
+  qpms_trans_calculator_get_trans_array_lc3p(c, s_2_1[0], bspec2, bspec1->n,
+      bspec1, 1, k, pos2, pos1);
+
+  qpms_trans_calculator_get_trans_array_lc3p(c, s_1_2[0], bspec1, bspec2->n,
+      bspec2, 1, k, pos1, pos2);
+
+  for(size_t R = 0; R < bspec2->n; ++R) {
+    for(size_t C = 0; C < bspec1->n; ++C)
+      printf("%.3lg+%.3lgj\t", creal(s_2_1[R][C]), cimag(s_2_1[C][R]));
+    putchar('\n');
+  }
+  putchar('\n');
+
+  for(size_t R = 0; R < bspec1->n; ++R) {
+    for(size_t C = 0; C < bspec2->n; ++C)
+      printf("%.3lg+%.3lgj\t", creal(s_1_2[R][C]), cimag(s_1_2[R][C]));
+    putchar('\n');
+  }
+
+  qpms_trans_calculator_free(c);
+  qpms_vswf_set_spec_free(bspec1);
+  qpms_vswf_set_spec_free(bspec2);
+  return 0;
+}