ss build translation matrix: enable other Bessel funcs.

Former-commit-id: 66d4993430322b7de00a408406b89972b7dd7c48

qpms/qpms_c.pyx

@@ -1571,12 +1571,12 @@ cdef class ScatteringSystem:
             qpms_scatsys_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.s, iri, k)
         return target
 
-    def translation_matrix_full(self, double k):
+    def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS):
         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)
+        qpms_scatsys_build_translation_matrix_full_e(&target_view[0][0], self.s, k, J)
         return target
     
     def fullvec_psizes(self):

qpms/qpms_cdefs.pxd

@@ -367,6 +367,8 @@ cdef extern from "scatsystem.h":
             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_translation_matrix_full_e(cdouble *target,
+            const qpms_scatsys_t *ss, double k, qpms_bessel_t J)
     cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target,
             const qpms_scatsys_t *ss, qpms_iri_t iri, double k)
     cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(

qpms/scatsystem.c

@@ -964,6 +964,18 @@ complex double *qpms_scatsys_build_translation_matrix_full(
     const qpms_scatsys_t *ss,
     double k ///< Wave number to use in the translation matrix.
     )
+{
+  return qpms_scatsys_build_translation_matrix_full_e(
+      target, ss, k, QPMS_HANKEL_PLUS);
+}
+
+complex double *qpms_scatsys_build_translation_matrix_full_e(
+    /// 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.
+    qpms_bessel_t J ///< Bessel function type.
+    )
 {
   const size_t full_len = ss->fecv_size;
   if(!target)
@@ -982,7 +994,7 @@ complex double *qpms_scatsys_build_translation_matrix_full(
           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));      
+                k, posR, posC, J));      
         }
         fullvec_offsetC += bspecC->n;
       }
@@ -1025,7 +1037,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
           QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
                 tmp, // tmp is S(piR<-piC)
                 bspecR, bspecC->n, bspecC, 1,
-                k, posR, posC));      
+                k, posR, posC, QPMS_HANKEL_PLUS));      
           cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
               bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 
               &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
@@ -1117,7 +1129,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
             QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
                   Sblock, // Sblock is S(piR->piC)
                   bspecR, bspecC->n, bspecC, 1,
-                  k, posR, posC));
+                  k, posR, posC, QPMS_HANKEL_PLUS));
 
 						cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
 								bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
@@ -1236,7 +1248,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
               QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
                     Sblock, // Sblock is S(piR->piC)
                     bspecR, bspecC->n, bspecC, 1,
-                    k, posR, posC));
+                    k, posR, posC, QPMS_HANKEL_PLUS));
 
               cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                   bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
@@ -1368,7 +1380,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
               QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
                     Sblock, // Sblock is S(piR->piC)
                     bspecR, bspecC->n, bspecC, 1,
-                    a->k, posR, posC));
+                    a->k, posR, posC, QPMS_HANKEL_PLUS));
 
               SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                   bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,

qpms/scatsystem.h

@@ -234,6 +234,18 @@ complex double *qpms_scatsys_build_translation_matrix_full(
 		const qpms_scatsys_t *ss,
 		double k ///< Wave number to use in the translation matrix.
 		);
+
+/// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
+/** Might be useful for evaluation of cross sections and testing.
+ */
+complex double *qpms_scatsys_build_translation_matrix_full_e(
+		/// 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.
+		qpms_bessel_t J
+		);
+
 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/translations.c

@@ -1119,8 +1119,6 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
   return retval;
 }
 
-/// Version with \a k and cartesian particle positions
-/// and with automatic \a kdlj = false and \a J = QPMS_HANKEL_PLUS.
 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( 
     const qpms_trans_calculator *c,
     complex double *target,
@@ -1128,7 +1126,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
     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,
-    double k, cart3_t destpos, cart3_t srcpos
+    double k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J
     /// Workspace has to be at least 2 * c->neleme**2 long
     )
 {
@@ -1136,7 +1134,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
   kdlj.r *= k;
   return qpms_trans_calculator_get_trans_array(c, target,
       destspec, deststride, srcspec, srcstride, kdlj, 
-      false, QPMS_HANKEL_PLUS);
+      false, J);
 }
 
 #ifdef LATTICESUMS31

qpms/translations.h

@@ -155,7 +155,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
 		sph_t kdlj, bool r_ge_d, qpms_bessel_t J);
 
 /// Version with \a k and cartesian particle positions
-/// and with automatic \a kdlj = false and \a J = QPMS_HANKEL_PLUS.
+/// and with automatic \a r_ge_d = `false`.
 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
 		const qpms_trans_calculator *c,
 		complex double *target, 
@@ -163,7 +163,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
 		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,
-		double k, cart3_t destpos, cart3_t srcpos
+		double k, cart3_t destpos, cart3_t srcpos,
+		qpms_bessel_t J
 		/// Workspace has to be at least 2 * c->neleme**2 long
 		);