From d31b6c4cedd317375f16606e428b691d0b1bcead Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 21 Feb 2019 03:48:35 +0000 Subject: [PATCH] scatsystem: T-matrix manipulation. Former-commit-id: 978a4d3b7de4f3ba3b606e8b6e8c6c7c871ba50a --- qpms/indexing.h | 10 ++++- qpms/scatsystem.c | 100 ++++++++++++++++++++++++++++++++++++++++++++++ qpms/scatsystem.h | 65 +++++++++++++++++++++++++++--- 3 files changed, 168 insertions(+), 7 deletions(-) create mode 100644 qpms/scatsystem.c diff --git a/qpms/indexing.h b/qpms/indexing.h index 267eddb..89398b2 100644 --- a/qpms/indexing.h +++ b/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 diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c new file mode 100644 index 0000000..f6a3881 --- /dev/null +++ b/qpms/scatsystem.c @@ -0,0 +1,100 @@ +#include "scatsystem.h" + +#include + +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; +} + + diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 4ea278e..6680b77 100644 --- a/qpms/scatsystem.h +++ b/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