From 2af9b83e355bc8b9703e1db40da3821fd6aa6876 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 18 Mar 2019 15:04:21 +0200 Subject: [PATCH] Move T-matrix -related C code to scatsystem.[ch] Former-commit-id: e9a162c14fe5d91281d79d9af08014787cd7ed13 --- qpms/qpms_cdefs.pxd | 50 ++--- qpms/qpms_types.h | 18 ++ qpms/scatsystem.c | 304 +---------------------------- qpms/scatsystem.h | 298 +---------------------------- qpms/tmatrices.c | 456 ++++++++++++++++++++++++++++++++++++++++++++ qpms/tmatrices.h | 286 ++++++++++++++++++++++++++- qpms/tmatrix_io.c | 145 -------------- setup.py | 2 +- 8 files changed, 783 insertions(+), 776 deletions(-) create mode 100644 qpms/tmatrices.c delete mode 100644 qpms/tmatrix_io.c diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 0ccfb82..be67b51 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -75,6 +75,10 @@ cdef extern from "qpms_types.h": ctypedef int qpms_gmi_t ctypedef int qpms_iri_t ctypedef const char * qpms_permutation_t + struct qpms_tmatrix_t: + qpms_vswf_set_spec_t *spec + cdouble *m + bint owns_m # FIXME in fact bool # maybe more if needed cdef extern from "indexing.h": @@ -221,17 +225,7 @@ cdef extern from "gsl/gsl_interp.h": const gsl_interp_type *gsl_interp_cspline # ^^^ These are probably the only relevant ones. -cdef extern from "scatsystem.h": - struct qpms_tmatrix_t: - qpms_vswf_set_spec_t *spec - cdouble *m - bint owns_m # FIXME in fact bool - struct qpms_particle_t: - cart3_t pos - const qpms_tmatrix_t *tmatrix - struct qpms_particle_tid_t: - cart3_t pos - qpms_ss_tmi_t tmatrix_id +cdef extern from "tmatrices.h": struct qpms_tmatrix_interpolator_t: const qpms_vswf_set_spec_t *bspec void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp) @@ -239,19 +233,6 @@ cdef extern from "scatsystem.h": qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, double *freqs, const qpms_tmatrix_t *tmatrices_array, const gsl_interp_type *iptype) void qpms_tmatrix_free(qpms_tmatrix_t *tmatrix) - struct qpms_scatsys_t: - qpms_tmatrix_t **tm - qpms_ss_tmi_t tm_count - qpms_particle_tid_t *p - qpms_ss_pi_t p_count - # We shouldn't need more to construct a symmetric scatsystem ^^^ - size_t fecv_size - size_t *saecv_sizes - const qpms_finite_group_t *sym - qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) - void qpms_scatsys_free(qpms_scatsys_t *s) - qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI - qpms_scatsys_t *qpms_scatsys_load(char *path) #NI qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B, const double rtol, const double atol) qpms_errno_t qpms_symmetrise_tmdata_irot3arr( @@ -277,6 +258,27 @@ cdef extern from "scatsystem.h": qpms_errno_t qpms_load_scuff_tmatrix(const char *path, const qpms_vswf_set_spec_t *bspec, size_t *n, double **freqs, double **freqs_su, qpms_tmatrix_t **tmatrices_array, cdouble **tmdata) + +cdef extern from "scatsystem.h": + struct qpms_particle_t: + cart3_t pos + const qpms_tmatrix_t *tmatrix + struct qpms_particle_tid_t: + cart3_t pos + qpms_ss_tmi_t tmatrix_id + struct qpms_scatsys_t: + qpms_tmatrix_t **tm + qpms_ss_tmi_t tm_count + qpms_particle_tid_t *p + qpms_ss_pi_t p_count + # We shouldn't need more to construct a symmetric scatsystem ^^^ + size_t fecv_size + size_t *saecv_sizes + const qpms_finite_group_t *sym + qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) + void qpms_scatsys_free(qpms_scatsys_t *s) + qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI + qpms_scatsys_t *qpms_scatsys_load(char *path) #NI cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed, const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 58a08b3..4c1585f 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -330,6 +330,24 @@ typedef int qpms_iri_t; */ typedef const char * qpms_permutation_t; +/// A T-matrix. +/** In the future, I might rather use a more abstract approach in which T-matrix + * is a mapping (function) of the field expansion coefficients. + * So the interface might change. + * For now, let me stick to the square dense matrix representation. + */ +typedef struct qpms_tmatrix_t { + /** \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. + bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free() +} qpms_tmatrix_t; #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l)) #endif // QPMS_TYPES diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index b76a54d..f3d421f 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -6,7 +6,6 @@ #include "vswf.h" #include "groups.h" #include "symmetries.h" -#include #include #include #include "vectors.h" @@ -14,314 +13,13 @@ #include #include "qpms_error.h" #include "translations.h" +#include "tmatrices.h" #include #define SQ(x) ((x)*(x)) #define QPMS_SCATSYS_LEN_RTOL 1e-13 #define QPMS_SCATSYS_TMATRIX_ATOL 1e-14 #define QPMS_SCATSYS_TMATRIX_RTOL 1e-12 -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(); - t->owns_m = true; - } - return t; -} - -qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) { - qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); - size_t n = T->spec->n; - for(size_t i = 0; i < n*n; ++i) - t->m = T->m; - return t; -} - -void qpms_tmatrix_free(qpms_tmatrix_t *t){ - if(t && t->owns_m) free(t->m); - free(t); -} - -qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( - 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_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_errno_t qpms_symmetrise_tmdata_irot3arr( - complex double *tmdata, const size_t tmcount, - const qpms_vswf_set_spec_t *bspec, - const size_t n_symops, const qpms_irot3_t *symops) { - const size_t n = bspec->n; - qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec); - complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops); - if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "malloc() failed."); - for (size_t i = 0; i < n_symops; ++i) - qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]); - complex double tmp[n][n]; - const complex double one = 1, zero = 0; - for (size_t tmi = 0; tmi < tmcount; ++tmi) { - // Move the data in tmcopy; we will then write the sum directly into tmdata. - memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double)); - memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double)); - for (size_t i = 0; i < n_symops; ++i) { - const complex double *const M = symop_matrices + i*n*n; - // tmp = M T - cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n); - // tmdata[...] += tmp M* = M T M* - cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, - n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n); - } - for (size_t ii = 0; ii < n*n; ++ii) - tmdata[n*n*tmi+ii] /= n_symops; - } - free(symop_matrices); - qpms_tmatrix_free(tmcopy); - return QPMS_SUCCESS; -} - -qpms_errno_t qpms_symmetrise_tmdata_finite_group( - complex double *tmdata, const size_t tmcount, - const qpms_vswf_set_spec_t *bspec, - const qpms_finite_group_t *pointgroup) { - if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "This function requires pointgroup->rep3d to be set correctly!"); - return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec, - pointgroup->order, pointgroup->rep3d); -} - -qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( - qpms_tmatrix_t *T, - size_t n_symops, - const qpms_irot3_t *symops - ) { - if(qpms_symmetrise_tmdata_irot3arr(T->m, 1, - T->spec, n_symops, symops) != QPMS_SUCCESS) - return NULL; - else return T; -} - -qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( - qpms_tmatrix_t *T, - const qpms_finite_group_t *pointgroup - ) { - if(qpms_symmetrise_tmdata_finite_group(T->m, 1, - T->spec, pointgroup) != QPMS_SUCCESS) - return NULL; - else return T; -} - -qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace( - qpms_tmatrix_t *T, - const complex double *M - ) -{ - qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M); - const size_t n = T->spec->n; - for(size_t i = 0; i < n*n; ++i) - T->m[i] = 0.5 * (t->m[i] + T->m[i]); - qpms_tmatrix_free(t); - 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_copy(T); - return qpms_tmatrix_symmetrise_C_inf_inplace(t); -} - -qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) { - 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) - ;// No-op // 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_copy(T); - return qpms_tmatrix_symmetrise_C_N_inplace(t, N); -} - -qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) { - 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; -} - -bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B, - const double rtol, const double atol) -{ - if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec)) - return false; - if (A->m == B->m) - return true; - const size_t n = A->spec->n; - for (size_t i = 0; i < n*n; ++i) { - const double tol = atol + rtol * (cabs(B->m[i])); - if ( cabs(B->m[i] - A->m[i]) > tol ) - return false; - } - return true; -} - - -qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount, - const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec - ) { - if (incount <= 0) return NULL; - qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t)); - /* - if (copy_bspec) { - ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec); - ip->owns_bspec = true; - } - else { - */ - ip->bspec = ta[0].spec; - // ip->owns_bspec = false; - //} - const size_t n = ip->bspec->n; - - // check if all matrices have the same bspec - for (size_t i = 0; i < incount; ++i) - if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec)) - abort(); - - if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort(); - if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort(); - for (size_t row = 0; row < n; ++row) - for (size_t col = 0; col < n; ++col) { - double y_real[incount], y_imag[incount]; - bool n0_real = false, n0_imag = false; - for (size_t i = 0; i < incount; ++i) { - complex double telem = ta[i].m[n * row + col]; - if ((y_real[i] = creal(telem))) n0_real = true; - if ((y_imag[i] = cimag(telem))) n0_imag = true; - } - if (n0_real) { - gsl_spline *s = - ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount); - if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort(); - } - else ip->splines_real[n * row + col] = NULL; - if (n0_imag) { - gsl_spline *s = - ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount); - if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort(); - } - else ip->splines_imag[n * row + col] = NULL; - } - return ip; -} - -void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) { - if (ip) { - const size_t n = ip->bspec->n; - for (size_t i = 0; i < n*n; ++i) { - if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]); - if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]); - } - //if (ip->owns_bspec) - // qpms_vswf_set_spec_free(ip->bspec); - free(ip); - } -} - -qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) { - qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec); - const size_t n = ip->bspec->n; - for (size_t i = 0; i < n*n; ++i){ - if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/); - if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/); - } - return t; -} - - // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() ------------- diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index f5bedb3..8609185 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -12,44 +12,6 @@ #define QPMS_SCATSYSTEM_H #include "qpms_types.h" #include -#include -#include // only because of qpms_read_scuff_tmatrix() - -/// A T-matrix. -/** In the future, I might rather use a more abstract approach in which T-matrix - * is a mapping (function) of the field expansion coefficients. - * So the interface might change. - * For now, let me stick to the square dense matrix representation. - */ -typedef struct qpms_tmatrix_t { - /** \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. - bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free() -} qpms_tmatrix_t; - -struct qpms_finite_group_t; -typedef struct qpms_finite_group_t qpms_finite_group_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); -} - -/// Initialises a zero T-matrix. -qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec); - -/// Copies a T-matrix, allocating new array for the T-matrix data. -qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T); - -/// 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 { @@ -58,227 +20,8 @@ typedef struct qpms_particle_t { const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t. } qpms_particle_t; -/// Check T-matrix equality/similarity. -/** - * This function actually checks for identical vswf specs. - * TODO define constants with "default" atol, rtol for this function. - */ -bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, - const double rtol, const double atol); - -/// 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 - ); - -/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data. -/** 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_inplace( - 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 - ); - -/// Symmetrizes a T-matrix with an involution symmetry operation, rewriting the original one. -/** 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_inplace( - 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, rewriting the original one. -/** - * \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f] - */ -qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace( - qpms_tmatrix_t *T ///< the original T-matrix - ); -/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix, rewriting the original one. -/** - * \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_inplace( - qpms_tmatrix_t *T, ///< the original T-matrix - int N ///< number of z-axis rotations in the group - ); - -/// Reads an open scuff-tmatrix generated file. -/** - * \a *freqs, \a *freqs_su, \a *tmatrices_array and \a *tmdata - * arrays are allocated by this function - * and have to be freed by the caller after use. - * \a freqs_su and \a tmatrices_array can be NULL, in that case - * the respective arrays are not filled nor allocated. - * - * The contents of tmatrices_array is NOT - * supposed to be freed element per element. - * - * TODO more checks and options regarding NANs etc. - * - */ -qpms_errno_t qpms_load_scuff_tmatrix( - const char *path, ///< Path to the TMatrix file - const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec - size_t *n, ///< Number of successfully loaded t-matrices - double **freqs, ///< Frequencies in SI units.. - double **freqs_su, ///< Frequencies in SCUFF units (optional). - /// The resulting T-matrices (optional). - qpms_tmatrix_t **tmatrices_array, - complex double **tmdata ///< The T-matrices raw contents - ); -/// Loads a scuff-tmatrix generated file. -/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a - * path instead of open FILE. - */ -qpms_errno_t qpms_read_scuff_tmatrix( - FILE *f, ///< An open stream with the T-matrix data. - const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec - size_t *n, ///< Number of successfully loaded t-matrices - double **freqs, ///< Frequencies in SI units. - double **freqs_su, ///< Frequencies in SCUFF units (optional). - /// The resulting T-matrices (optional). - qpms_tmatrix_t **tmatrices_array, - /// The T-matrices raw contents. - /** The coefficient of outgoing wave defined by - * \a bspec->ilist[desti] as a result of incoming wave - * \a bspec->ilist[srci] at frequency \a (*freqs)[fi] - * is accessed via - * (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci]. - */ - complex double ** tmdata - ); - -/// In-place application of point group elements on raw T-matrix data. -/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix(). - * The \a symops array should always contain all elements of a finite - * point (sub)group, including the identity operation. - * - * TODO more doc. - */ -qpms_errno_t qpms_symmetrise_tmdata_irot3arr( - complex double *tmdata, const size_t tmcount, - const qpms_vswf_set_spec_t *bspec, - size_t n_symops, - const qpms_irot3_t *symops - ); - -/// In-place application of a point group on raw T-matrix data. -/** This does the same as qpms_symmetrise_tmdata_irot3arr(), - * but takes a valid finite point group as an argument. - * - * TODO more doc. - */ -qpms_errno_t qpms_symmetrise_tmdata_finite_group( - complex double *tmdata, const size_t tmcount, - const qpms_vswf_set_spec_t *bspec, - const qpms_finite_group_t *pointgroup - ); - -/// In-place application of point group elements on a T-matrix. -/** The \a symops array should always contain all elements of a finite - * point (sub)group, including the identity operation. - * - * TODO more doc. - */ -qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( - qpms_tmatrix_t *T, - size_t n_symops, - const qpms_irot3_t *symops - ); - -/// In-place application of point group elements on a T-matrix. -/** This does the same as qpms_tmatrix_symmetrise_irot3arr(), - * but takes a valid finite point group as an argument. - * - * TODO more doc. - */ -qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( - qpms_tmatrix_t *T, - const qpms_finite_group_t *pointgroup - ); - -/* Fuck this, include the whole -typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct -typedef struct gsl_interp_type gsl_interp_type; -extern const gsl_interp_type * gsl_interp_linear; -extern const gsl_interp_type * gsl_interp_polynomial; -extern const gsl_interp_type * gsl_interp_cspline; -extern const gsl_interp_type * gsl_interp_cspline_periodic; -extern const gsl_interp_type * gsl_interp_akima; -extern const gsl_interp_type * gsl_interp_akima_periodic; -extern const gsl_interp_type * gsl_interp_steffen; -*/ - -// struct gsl_interp_accel; // use if lookup proves to be too slow -typedef struct qpms_tmatrix_interpolator_t { - const qpms_vswf_set_spec_t *bspec; - //bool owns_bspec; - gsl_spline **splines_real; ///< There will be a spline object for each nonzero element - gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element - // gsl_interp_accel **accel_real; - // gsl_interp_accel **accel_imag; -} qpms_tmatrix_interpolator_t; - -/// Free a T-matrix interpolator. -void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp); - -/// Evaluate a T-matrix interpolated value. -/** The result is to be freed using qpms_tmatrix_free().*/ -qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq); - -/// Create a T-matrix interpolator from frequency and T-matrix arrays. -qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided. - const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers! - const gsl_interp_type *iptype - //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix. - /*, ...? */); - +struct qpms_finite_group_t; +typedef struct qpms_finite_group_t qpms_finite_group_t; /// A particle, defined by its T-matrix INDEX and position, to be used in qpms_scatsys_t. typedef struct qpms_particle_tid_t { @@ -571,41 +314,4 @@ complex double *qpms_orbit_irrep_basis( /// The index of the irreducible representation of sym. const qpms_iri_t iri); - - -#if 0 -// Abstract types that describe T-matrix/particle/scatsystem symmetries -// To be implemented later. See also the thoughts in the beginning of groups.h. - -typedef *char qpms_tmatrix_id_t; ///< Maybe I want some usual integer type instead. - -///Abstract T-matrix type draft. -/** - * TODO. - */ -typedef struct qpms_abstract_tmatrix_t{ - qpms_tmatrix_id_t id; - /// Generators of the discrete point group under which T-matrix is invariant. - qpms_irot3_t *invar_gens; - /// Length of invar_gens. - qpms_gmi_t invar_gens_size; - -} qpms_abstract_tmatrix_t; - - -typedef struct qpms_abstract_particle_t{ -} qpms_abstract_particle_t; - -/// An abstract particle, defined by its position and abstract T-matrix. -typedef struct qpms_abstract_particle_t { - cart3_t pos; ///< Particle position in cartesian coordinates. - const qpms_abstract_tmatrix_t *tmatrix; ///< T-matrix; not owned by this. -} qpms_abstract_particle_t; - - -/** This is just an alias, as the same index can be used for - * abstract T-matrices as well. - */ -typedef qpms_particle_tid_t qpms_abstract_particle_tid_t; -#endif // 0 #endif //QPMS_SCATSYSTEM_H diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c new file mode 100644 index 0000000..4dfd6a1 --- /dev/null +++ b/qpms/tmatrices.c @@ -0,0 +1,456 @@ +#define _POSIX_C_SOURCE 200809L // for getline() +#include +#include +#include +#include +#include "scatsystem.h" +#include "indexing.h" +#include "vswf.h" +#include "groups.h" +#include "symmetries.h" +#include +#include +#include +#include "vectors.h" +#include "wigner.h" +#include +#include "qpms_error.h" +#include "tmatrices.h" + +#define HBAR (1.05457162825e-34) +#define ELECTRONVOLT (1.602176487e-19) +#define SCUFF_OMEGAUNIT (3e14) + +#define SQ(x) ((x)*(x)) +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(); + t->owns_m = true; + } + return t; +} + +qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) { + qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); + size_t n = T->spec->n; + for(size_t i = 0; i < n*n; ++i) + t->m = T->m; + return t; +} + +void qpms_tmatrix_free(qpms_tmatrix_t *t){ + if(t && t->owns_m) free(t->m); + free(t); +} + +qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( + 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_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_errno_t qpms_symmetrise_tmdata_irot3arr( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const size_t n_symops, const qpms_irot3_t *symops) { + const size_t n = bspec->n; + qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec); + complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops); + if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "malloc() failed."); + for (size_t i = 0; i < n_symops; ++i) + qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]); + complex double tmp[n][n]; + const complex double one = 1, zero = 0; + for (size_t tmi = 0; tmi < tmcount; ++tmi) { + // Move the data in tmcopy; we will then write the sum directly into tmdata. + memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double)); + memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double)); + for (size_t i = 0; i < n_symops; ++i) { + const complex double *const M = symop_matrices + i*n*n; + // tmp = M T + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n); + // tmdata[...] += tmp M* = M T M* + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, + n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n); + } + for (size_t ii = 0; ii < n*n; ++ii) + tmdata[n*n*tmi+ii] /= n_symops; + } + free(symop_matrices); + qpms_tmatrix_free(tmcopy); + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_symmetrise_tmdata_finite_group( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *pointgroup) { + if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "This function requires pointgroup->rep3d to be set correctly!"); + return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec, + pointgroup->order, pointgroup->rep3d); +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( + qpms_tmatrix_t *T, + size_t n_symops, + const qpms_irot3_t *symops + ) { + if(qpms_symmetrise_tmdata_irot3arr(T->m, 1, + T->spec, n_symops, symops) != QPMS_SUCCESS) + return NULL; + else return T; +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( + qpms_tmatrix_t *T, + const qpms_finite_group_t *pointgroup + ) { + if(qpms_symmetrise_tmdata_finite_group(T->m, 1, + T->spec, pointgroup) != QPMS_SUCCESS) + return NULL; + else return T; +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace( + qpms_tmatrix_t *T, + const complex double *M + ) +{ + qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M); + const size_t n = T->spec->n; + for(size_t i = 0; i < n*n; ++i) + T->m[i] = 0.5 * (t->m[i] + T->m[i]); + qpms_tmatrix_free(t); + 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_copy(T); + return qpms_tmatrix_symmetrise_C_inf_inplace(t); +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) { + 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) + ;// No-op // 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_copy(T); + return qpms_tmatrix_symmetrise_C_N_inplace(t, N); +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) { + 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; +} + +bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B, + const double rtol, const double atol) +{ + if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec)) + return false; + if (A->m == B->m) + return true; + const size_t n = A->spec->n; + for (size_t i = 0; i < n*n; ++i) { + const double tol = atol + rtol * (cabs(B->m[i])); + if ( cabs(B->m[i] - A->m[i]) > tol ) + return false; + } + return true; +} + +qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount, + const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec + ) { + if (incount <= 0) return NULL; + qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t)); + /* + if (copy_bspec) { + ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec); + ip->owns_bspec = true; + } + else { + */ + ip->bspec = ta[0].spec; + // ip->owns_bspec = false; + //} + const size_t n = ip->bspec->n; + + // check if all matrices have the same bspec + for (size_t i = 0; i < incount; ++i) + if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec)) + abort(); + + if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort(); + if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort(); + for (size_t row = 0; row < n; ++row) + for (size_t col = 0; col < n; ++col) { + double y_real[incount], y_imag[incount]; + bool n0_real = false, n0_imag = false; + for (size_t i = 0; i < incount; ++i) { + complex double telem = ta[i].m[n * row + col]; + if ((y_real[i] = creal(telem))) n0_real = true; + if ((y_imag[i] = cimag(telem))) n0_imag = true; + } + if (n0_real) { + gsl_spline *s = + ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount); + if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort(); + } + else ip->splines_real[n * row + col] = NULL; + if (n0_imag) { + gsl_spline *s = + ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount); + if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort(); + } + else ip->splines_imag[n * row + col] = NULL; + } + return ip; +} + +void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) { + if (ip) { + const size_t n = ip->bspec->n; + for (size_t i = 0; i < n*n; ++i) { + if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]); + if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]); + } + //if (ip->owns_bspec) + // qpms_vswf_set_spec_free(ip->bspec); + free(ip); + } +} + +qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) { + qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec); + const size_t n = ip->bspec->n; + for (size_t i = 0; i < n*n; ++i){ + if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/); + if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/); + } + return t; +} + + + +double qpms_SU2eV(double e_SU) { + return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR); +} + +double qpms_SU2SI(double e_SU) { + return e_SU * SCUFF_OMEGAUNIT; +} + +/// TODO doc and more checks +qpms_errno_t qpms_read_scuff_tmatrix( + FILE *f, ///< file handle + const qpms_vswf_set_spec_t * bs, ///< VSWF set spec + size_t *const n, ///< Number of successfully loaded t-matrices + double* *const freqs, ///< Frequencies in SI units + double* *const freqs_su, ///< Frequencies in SCUFF units (optional) + qpms_tmatrix_t* *const tmatrices_array, ///< The resulting T-matrices (optional). + complex double* *const tmdata + ) { + if (!(freqs && n && tmdata)) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "freqs, n, and tmdata are mandatory arguments and must not be NULL."); + if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" + " norm supported right now."); + int n_alloc = 128; // First chunk to allocate + *n = 0; + *freqs = malloc(n_alloc * sizeof(double)); + if (freqs_su) *freqs_su = malloc(n_alloc * sizeof(double)); + *tmdata = malloc(sizeof(complex double) * bs->n * bs->n * n_alloc); + if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "malloc() failed."); + size_t linebufsz = 256; + char *linebuf = malloc(linebufsz); + ssize_t readchars; + double lastfreq_su = NAN; + while((readchars = getline(&linebuf, &linebufsz, f)) != -1) { + if (linebuf[0] == '#') continue; + int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta; + double currentfreq_su, tr, ti; + if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf", + ¤tfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha, + &Beta, &LBeta, &MBeta, &PBeta, &tr, &ti)) + abort(); // Malformed T-matrix file + if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix + ++*n; + lastfreq_su = currentfreq_su; + if(*n > n_alloc) { + n_alloc *= 2; + *freqs = realloc(*freqs, n_alloc * sizeof(double)); + if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); + *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); + if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "realloc() failed."); + } + if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su; + (*freqs)[*n-1] = qpms_SU2SI(currentfreq_su); + + for(size_t i = 0; i < bs->n * bs->n; ++i) + (*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN; + } + qpms_vswf_type_t TAlpha, TBeta; + switch(PAlpha) { + case 0: TAlpha = QPMS_VSWF_MAGNETIC; break; + case 1: TAlpha = QPMS_VSWF_ELECTRIC; break; + default: assert(0); + } + switch(PBeta) { + case 0: TBeta = QPMS_VSWF_MAGNETIC; break; + case 1: TBeta = QPMS_VSWF_ELECTRIC; break; + default: assert(0); + } + qpms_uvswfi_t srcui = qpms_tmn2uvswfi(TAlpha, MAlpha, LAlpha), + destui = qpms_tmn2uvswfi(TBeta, MBeta, LBeta); + ssize_t srci = qpms_vswf_set_spec_find_uvswfi(bs, srcui), + desti = qpms_vswf_set_spec_find_uvswfi(bs, destui); + if (srci == -1 || desti == -1) + /* This element has not been requested in bs->ilist. */ + continue; + else + (*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti; + } + free(linebuf); + // free some more memory + n_alloc = *n; + *freqs = realloc(*freqs, n_alloc * sizeof(double)); + if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); + if (tmatrices_array) *tmatrices_array = realloc(*tmatrices_array, n_alloc * sizeof(qpms_tmatrix_t)); + *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); + if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "realloc() failed."); + if (tmatrices_array) { + *tmatrices_array = malloc(n_alloc * sizeof(qpms_tmatrix_t)); + if (!*tmatrices_array) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "malloc() failed."); + for (size_t i = 0; i < *n; ++i) { + (*tmatrices_array)[i].spec = bs; + (*tmatrices_array)[i].m = (*tmdata) + i * bs->n * bs->n; + (*tmatrices_array)[i].owns_m = false; + } + } + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_load_scuff_tmatrix( + const char *path, ///< file path + const qpms_vswf_set_spec_t * bs, ///< VSWF set spec + size_t *const n, ///< Number of successfully loaded t-matrices + double **const freqs, ///< Frequencies in SI units + double ** const freqs_su, ///< Frequencies in SCUFF units (optional) + qpms_tmatrix_t ** const tmatrices_array, ///< The resulting T-matrices (optional). + complex double ** const tmdata + ) { + FILE *f = fopen(path, "r"); + if (!f) + qpms_pr_error_at_line(__FILE__, __LINE__, __func__, + "Could not open T-matrix file %s", path); + qpms_errno_t retval = + qpms_read_scuff_tmatrix(f, bs, n, freqs, freqs_su, tmatrices_array, tmdata); + if(fclose(f)) qpms_pr_error_at_line(__FILE__, __LINE__, __func__, + "Could not close the T-matrix file %s (well, that's weird, " + "since it's read only).", path); + + return retval; +} + diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 305c3e1..d8bdcec 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -1,12 +1,284 @@ +/* \file tmatrices.h + * \brief T-matrices for scattering systems. + */ #ifndef TMATRICES_H #define TMATRICES_H -/* TODO - * This file will contain declarations of functions providing - * a) Mie T-matrix for spherical particle - * i) using Drude model - * ii) using interpolated material data - * b) T-matrix from scuff-tmatrix output, using interpolation - */ +#include "qpms_types.h" +#include +struct qpms_finite_group_t; +typedef struct qpms_finite_group_t qpms_finite_group_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); +} + +/// Initialises a zero T-matrix. +qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec); + +/// Copies a T-matrix, allocating new array for the T-matrix data. +qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T); + +/// Destroys a T-matrix. +void qpms_tmatrix_free(qpms_tmatrix_t *t); + +/// Check T-matrix equality/similarity. +/** + * This function actually checks for identical vswf specs. + * TODO define constants with "default" atol, rtol for this function. + */ +bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, + const double rtol, const double atol); + +/// 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 + ); + +/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data. +/** 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_inplace( + 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 + ); + +/// Symmetrizes a T-matrix with an involution symmetry operation, rewriting the original one. +/** 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_inplace( + 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, rewriting the original one. +/** + * \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f] + */ +qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace( + qpms_tmatrix_t *T ///< the original T-matrix + ); +/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix, rewriting the original one. +/** + * \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_inplace( + qpms_tmatrix_t *T, ///< the original T-matrix + int N ///< number of z-axis rotations in the group + ); + +/// Reads an open scuff-tmatrix generated file. +/** + * \a *freqs, \a *freqs_su, \a *tmatrices_array and \a *tmdata + * arrays are allocated by this function + * and have to be freed by the caller after use. + * \a freqs_su and \a tmatrices_array can be NULL, in that case + * the respective arrays are not filled nor allocated. + * + * The contents of tmatrices_array is NOT + * supposed to be freed element per element. + * + * TODO more checks and options regarding NANs etc. + * + */ +qpms_errno_t qpms_load_scuff_tmatrix( + const char *path, ///< Path to the TMatrix file + const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec + size_t *n, ///< Number of successfully loaded t-matrices + double **freqs, ///< Frequencies in SI units.. + double **freqs_su, ///< Frequencies in SCUFF units (optional). + /// The resulting T-matrices (optional). + qpms_tmatrix_t **tmatrices_array, + complex double **tmdata ///< The T-matrices raw contents + ); +/// Loads a scuff-tmatrix generated file. +/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a + * path instead of open FILE. + */ +qpms_errno_t qpms_read_scuff_tmatrix( + FILE *f, ///< An open stream with the T-matrix data. + const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec + size_t *n, ///< Number of successfully loaded t-matrices + double **freqs, ///< Frequencies in SI units. + double **freqs_su, ///< Frequencies in SCUFF units (optional). + /// The resulting T-matrices (optional). + qpms_tmatrix_t **tmatrices_array, + /// The T-matrices raw contents. + /** The coefficient of outgoing wave defined by + * \a bspec->ilist[desti] as a result of incoming wave + * \a bspec->ilist[srci] at frequency \a (*freqs)[fi] + * is accessed via + * (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci]. + */ + complex double ** tmdata + ); + +/// In-place application of point group elements on raw T-matrix data. +/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix(). + * The \a symops array should always contain all elements of a finite + * point (sub)group, including the identity operation. + * + * TODO more doc. + */ +qpms_errno_t qpms_symmetrise_tmdata_irot3arr( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + size_t n_symops, + const qpms_irot3_t *symops + ); + +/// In-place application of a point group on raw T-matrix data. +/** This does the same as qpms_symmetrise_tmdata_irot3arr(), + * but takes a valid finite point group as an argument. + * + * TODO more doc. + */ +qpms_errno_t qpms_symmetrise_tmdata_finite_group( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *pointgroup + ); + +/// In-place application of point group elements on a T-matrix. +/** The \a symops array should always contain all elements of a finite + * point (sub)group, including the identity operation. + * + * TODO more doc. + */ +qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( + qpms_tmatrix_t *T, + size_t n_symops, + const qpms_irot3_t *symops + ); + +/// In-place application of point group elements on a T-matrix. +/** This does the same as qpms_tmatrix_symmetrise_irot3arr(), + * but takes a valid finite point group as an argument. + * + * TODO more doc. + */ +qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( + qpms_tmatrix_t *T, + const qpms_finite_group_t *pointgroup + ); + +/* Fuck this, include the whole +typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct +typedef struct gsl_interp_type gsl_interp_type; +extern const gsl_interp_type * gsl_interp_linear; +extern const gsl_interp_type * gsl_interp_polynomial; +extern const gsl_interp_type * gsl_interp_cspline; +extern const gsl_interp_type * gsl_interp_cspline_periodic; +extern const gsl_interp_type * gsl_interp_akima; +extern const gsl_interp_type * gsl_interp_akima_periodic; +extern const gsl_interp_type * gsl_interp_steffen; +*/ + +// struct gsl_interp_accel; // use if lookup proves to be too slow +typedef struct qpms_tmatrix_interpolator_t { + const qpms_vswf_set_spec_t *bspec; + //bool owns_bspec; + gsl_spline **splines_real; ///< There will be a spline object for each nonzero element + gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element + // gsl_interp_accel **accel_real; + // gsl_interp_accel **accel_imag; +} qpms_tmatrix_interpolator_t; + +/// Free a T-matrix interpolator. +void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp); + +/// Evaluate a T-matrix interpolated value. +/** The result is to be freed using qpms_tmatrix_free().*/ +qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq); + +/// Create a T-matrix interpolator from frequency and T-matrix arrays. +qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided. + const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers! + const gsl_interp_type *iptype + //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix. + /*, ...? */); + + +#if 0 +// Abstract types that describe T-matrix/particle/scatsystem symmetries +// To be implemented later. See also the thoughts in the beginning of groups.h. + +typedef *char qpms_tmatrix_id_t; ///< Maybe I want some usual integer type instead. + +///Abstract T-matrix type draft. +/** + * TODO. + */ +typedef struct qpms_abstract_tmatrix_t{ + qpms_tmatrix_id_t id; + /// Generators of the discrete point group under which T-matrix is invariant. + qpms_irot3_t *invar_gens; + /// Length of invar_gens. + qpms_gmi_t invar_gens_size; + +} qpms_abstract_tmatrix_t; + + +typedef struct qpms_abstract_particle_t{ +} qpms_abstract_particle_t; + +/// An abstract particle, defined by its position and abstract T-matrix. +typedef struct qpms_abstract_particle_t { + cart3_t pos; ///< Particle position in cartesian coordinates. + const qpms_abstract_tmatrix_t *tmatrix; ///< T-matrix; not owned by this. +} qpms_abstract_particle_t; + + +/** This is just an alias, as the same index can be used for + * abstract T-matrices as well. + */ +typedef qpms_particle_tid_t qpms_abstract_particle_tid_t; +#endif // 0 #endif //TMATRICES_H diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c deleted file mode 100644 index ac77e26..0000000 --- a/qpms/tmatrix_io.c +++ /dev/null @@ -1,145 +0,0 @@ -#define _POSIX_C_SOURCE 200809L // for getline() -#include "scatsystem.h" -#include -#include -#include -#include "indexing.h" -#include "vswf.h" // qpms_vswf_set_spec_find_uvswfi() -#include "qpms_error.h" - - -#define HBAR (1.05457162825e-34) -#define ELECTRONVOLT (1.602176487e-19) -#define SCUFF_OMEGAUNIT (3e14) - -double qpms_SU2eV(double e_SU) { - return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR); -} - -double qpms_SU2SI(double e_SU) { - return e_SU * SCUFF_OMEGAUNIT; -} - -/// TODO doc and more checks -qpms_errno_t qpms_read_scuff_tmatrix( - FILE *f, ///< file handle - const qpms_vswf_set_spec_t * bs, ///< VSWF set spec - size_t *const n, ///< Number of successfully loaded t-matrices - double* *const freqs, ///< Frequencies in SI units - double* *const freqs_su, ///< Frequencies in SCUFF units (optional) - qpms_tmatrix_t* *const tmatrices_array, ///< The resulting T-matrices (optional). - complex double* *const tmdata - ) { - if (!(freqs && n && tmdata)) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "freqs, n, and tmdata are mandatory arguments and must not be NULL."); - if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" - " norm supported right now."); - int n_alloc = 128; // First chunk to allocate - *n = 0; - *freqs = malloc(n_alloc * sizeof(double)); - if (freqs_su) *freqs_su = malloc(n_alloc * sizeof(double)); - *tmdata = malloc(sizeof(complex double) * bs->n * bs->n * n_alloc); - if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "malloc() failed."); - size_t linebufsz = 256; - char *linebuf = malloc(linebufsz); - ssize_t readchars; - double lastfreq_su = NAN; - while((readchars = getline(&linebuf, &linebufsz, f)) != -1) { - if (linebuf[0] == '#') continue; - int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta; - double currentfreq_su, tr, ti; - if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf", - ¤tfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha, - &Beta, &LBeta, &MBeta, &PBeta, &tr, &ti)) - abort(); // Malformed T-matrix file - if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix - ++*n; - lastfreq_su = currentfreq_su; - if(*n > n_alloc) { - n_alloc *= 2; - *freqs = realloc(*freqs, n_alloc * sizeof(double)); - if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); - *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); - if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "realloc() failed."); - } - if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su; - (*freqs)[*n-1] = qpms_SU2SI(currentfreq_su); - - for(size_t i = 0; i < bs->n * bs->n; ++i) - (*tmdata)[(*n-1)*bs->n*bs->n + i] = NAN + I*NAN; - } - qpms_vswf_type_t TAlpha, TBeta; - switch(PAlpha) { - case 0: TAlpha = QPMS_VSWF_MAGNETIC; break; - case 1: TAlpha = QPMS_VSWF_ELECTRIC; break; - default: assert(0); - } - switch(PBeta) { - case 0: TBeta = QPMS_VSWF_MAGNETIC; break; - case 1: TBeta = QPMS_VSWF_ELECTRIC; break; - default: assert(0); - } - qpms_uvswfi_t srcui = qpms_tmn2uvswfi(TAlpha, MAlpha, LAlpha), - destui = qpms_tmn2uvswfi(TBeta, MBeta, LBeta); - ssize_t srci = qpms_vswf_set_spec_find_uvswfi(bs, srcui), - desti = qpms_vswf_set_spec_find_uvswfi(bs, destui); - if (srci == -1 || desti == -1) - /* This element has not been requested in bs->ilist. */ - continue; - else - (*tmdata)[(*n-1)*bs->n*bs->n + desti*bs->n + srci] = tr + I*ti; - } - free(linebuf); - // free some more memory - n_alloc = *n; - *freqs = realloc(*freqs, n_alloc * sizeof(double)); - if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); - if (tmatrices_array) *tmatrices_array = realloc(*tmatrices_array, n_alloc * sizeof(qpms_tmatrix_t)); - *tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); - if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "realloc() failed."); - if (tmatrices_array) { - *tmatrices_array = malloc(n_alloc * sizeof(qpms_tmatrix_t)); - if (!*tmatrices_array) - qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "malloc() failed."); - for (size_t i = 0; i < *n; ++i) { - (*tmatrices_array)[i].spec = bs; - (*tmatrices_array)[i].m = (*tmdata) + i * bs->n * bs->n; - (*tmatrices_array)[i].owns_m = false; - } - } - return QPMS_SUCCESS; -} - -qpms_errno_t qpms_load_scuff_tmatrix( - const char *path, ///< file path - const qpms_vswf_set_spec_t * bs, ///< VSWF set spec - size_t *const n, ///< Number of successfully loaded t-matrices - double **const freqs, ///< Frequencies in SI units - double ** const freqs_su, ///< Frequencies in SCUFF units (optional) - qpms_tmatrix_t ** const tmatrices_array, ///< The resulting T-matrices (optional). - complex double ** const tmdata - ) { - FILE *f = fopen(path, "r"); - if (!f) - qpms_pr_error_at_line(__FILE__, __LINE__, __func__, - "Could not open T-matrix file %s", path); - qpms_errno_t retval = - qpms_read_scuff_tmatrix(f, bs, n, freqs, freqs_su, tmatrices_array, tmdata); - if(fclose(f)) qpms_pr_error_at_line(__FILE__, __LINE__, __func__, - "Could not close the T-matrix file %s (well, that's weird, " - "since it's read only).", path); - - return retval; -} - - diff --git a/setup.py b/setup.py index 52775cf..cf067dd 100755 --- a/setup.py +++ b/setup.py @@ -27,7 +27,7 @@ qpms_c = Extension('qpms_c', 'qpms/scatsystem.c', 'qpms/vswf.c', # FIXME many things from vswf.c are not required by this module, but they have many dependencies (following in this list); maybe I want to move all the "basespec stuff" 'qpms/legendre.c', - 'qpms/tmatrix_io.c', + 'qpms/tmatrices.c', 'qpms/error.c'], extra_compile_args=['-std=c99','-ggdb', '-O3', '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required