scatsystem: T-matrix manipulation.

Former-commit-id: 978a4d3b7de4f3ba3b606e8b6e8c6c7c871ba50a

qpms/indexing.h

@@ -65,11 +65,19 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
 		qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
 	*t = u & 3;
 	qpms_y_sc_t y_sc = u / 4;
+	qpms_y2mn_sc_p(y_sc, m, n);
 	// Test validity
 	if (*t == 3) return QPMS_ERROR; // VSWF type code invalid, TODO WARN
 	if (*t && !y_sc) return QPMS_ERROR; // l == 0 for transversal wave, TODO WARN
-	qpms_y2mn_sc_p(y_sc, m, n);
 	return QPMS_SUCCESS;
 }
 
+/// Extract degree \a m from an universal VSWF index \a u.
+static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
+	qpms_vswf_type_t t; qpms_m_t m; qpms_l_t n;
+	qpms_uvswfi2tmn(u, &t,&m,&n);
+	return m;
+}
+
+
 #endif //QPMS_INDEXING_H

qpms/scatsystem.c

@@ -0,0 +1,100 @@
+#include "scatsystem.h"
+
+#include <cblas.h>
+
+qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
+  qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
+  if (!t) abort();
+  else {
+    t->spec = bspec;
+    size_t n = bspec->n;
+    t->m = calloc(n*n, sizeof(complex double));
+    if (!t->m) abort();
+  }
+  return t;
+}
+
+void qpms_tmatrix_free(qpms_tmatrix_t *t){
+  if(t) free(t->m);
+  free(t);
+}
+
+qpms_tmatrix_t *qpms_tmatrix_apply_symop(
+                const qpms_tmatrix_t *T, 
+                const complex double *M 
+                )
+{
+  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  const size_t n = T->spec->n;
+  complex double tmp[n][n];
+  // tmp = M T
+  const complex double one = 1, zero = 0;
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasNoTrans,
+      n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
+  // t->m = tmp M* = M T M*
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasConjTrans,
+      n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
+  return t;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
+                const qpms_tmatrix_t *T, 
+                const complex double *M 
+                )
+{
+  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  const size_t n = T->spec->n;
+  complex double tmp[n][n];
+  // tmp = M T
+  const complex double one = 1, zero = 0;
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasNoTrans,
+      n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
+  // t->m = tmp M* = M T M*
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasConjTrans,
+      n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
+  for(size_t i = 0; i < n*n; ++i)
+    t->m[i] = 0.5 * (t->m[i] + T->m[i]);
+  return t;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(const qpms_tmatrix_t *T) {
+  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  const size_t n = T->spec->n;
+  for (size_t row = 0; row < n; row++) {
+    qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
+    for (size_t col = 0; col < n; col++) {
+      qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
+      if (rm == cm)
+        t->m[n*row + col] = T->m[n*row + col];
+      else
+        t->m[n*row + col] = 0;
+    }
+  }
+  return t;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(const qpms_tmatrix_t *T, int N) {
+  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  const size_t n = T->spec->n;
+  for (size_t row = 0; row < n; row++) {
+    qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
+    for (size_t col = 0; col < n; col++) {
+      qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
+      if (((rm - cm) % N) == 0)
+        t->m[n*row + col] = T->m[n*row + col];
+      else
+        t->m[n*row + col] = 0;
+    }
+  }
+  return t;
+}
+
+

qpms/scatsystem.h

@@ -3,6 +3,7 @@
  */
 #ifndef QPMS_SCATSYSTEM_H
 #define QPMS_SCATSYSTEM_H
+#include "vswf.h"
 
 /// A T-matrix.
 /** In the future, I might rather use a more abstract approach in which T-matrix
@@ -11,18 +12,27 @@
  *  For now, let me stick to the square dense matrix representation.
  */
 typedef struct qpms_tmatrix_t {
-	const qpms_vswf_set_spec_t *spec; ///< NOT owned by qpms_tmatrix_t by default.
-	complex double *m; ///< Matrix elements 
+	/** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default.
+	 *
+	 *  Usually not checked for meaningfulness by the functions (methods),
+	 *  so the caller should take care that \a spec->ilist does not 
+	 *  contain any duplicities and that for each wave with order \a m
+	 *  there is also one with order \a −m.
+	 */
+	const qpms_vswf_set_spec_t *spec; 
+	complex double *m; ///< Matrix elements in row-major order.
 } qpms_tmatrix_t;
+
 /// Returns a pointer to the beginning of the T-matrix row \a rowno.
 static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
 	return t->m + (t->spec->n * rowno);
 }
-/// NOT IMPLEMENTED Initialises a zero T-matrix.
-qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
-/// NOT IMPLEMENTED Destroys a T-matrix.
-void qpms_tmatrix_free(qpms_tmatrix_t *t);
 
+/// Initialises a zero T-matrix.
+qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
+
+/// Destroys a T-matrix.
+void qpms_tmatrix_free(qpms_tmatrix_t *t);
 
 /// A particle, defined by its T-matrix and position.
 typedef struct qpms_particle_t {
@@ -31,4 +41,47 @@ typedef struct qpms_particle_t {
 	const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t.
 } qpms_particle_t;
 
+
+
+/// Creates a T-matrix from another matrix and a symmetry operation.
+/** The symmetry operation is expected to be a unitary (square) 
+ *  matrix \a M with the same
+ *  VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
+ *  correspond to CHECKME \f[ T' = MTM^\dagger \f] 
+ */
+qpms_tmatrix_t *qpms_tmatrix_apply_symop(
+		const qpms_tmatrix_t *T, ///< the original T-matrix
+		const complex double *M ///< the symmetry op matrix in row-major format
+		);
+
+/// Symmetrizes a T-matrix with an involution symmetry operation.
+/** The symmetry operation is expected to be a unitary (square) 
+ *  matrix \a M with the same
+ *  VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
+ *  correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f] 
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
+		const qpms_tmatrix_t *T, ///< the original T-matrix
+		const complex double *M ///< the symmetry op matrix in row-major format
+		);
+
+/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
+/**
+ *  \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(
+		const qpms_tmatrix_t *T ///< the original T-matrix
+		);
+/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix.
+/**
+ *  \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
+ *  T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
+ *  0 & (m-\mu)\mod N\ne0
+ *  \end{cases} . \f]
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
+		const qpms_tmatrix_t *T, ///< the original T-matrix
+		int N ///< number of z-axis rotations in the group
+		);
+
 #endif //QPMS_SCATSYSTEM_H