diff --git a/qpms/symmetries.c b/qpms/symmetries.c index f2dd050..24c473d 100644 --- a/qpms/symmetries.c +++ b/qpms/symmetries.c @@ -1,6 +1,7 @@ #include "symmetries.h" #include "tiny_inlines.h" #include "indexing.h" +#include "wigner.h" // TODO at some point, maybe support also other norms. // (perhaps just use qpms_normalisation_t_factor() at the right places) @@ -16,15 +17,25 @@ static inline void check_norm_compat(const qpms_vswf_set_spec_t *s) } } -complex double *qpms_zflip_uvswi_dense( - complex double *target, - const qpms_vswf_set_spec_t *bspec) -{ +// Used in the functions below to ensure memory allocation and checks for bspec validity +static inline complex double *ensure_alloc(complex double *target, + const qpms_vswf_set_spec_t *bspec) { check_norm_compat(bspec); const size_t n = bspec->n; if (target == NULL) target = malloc(n * n * sizeof(complex double)); if (target == NULL) abort(); + return target; +} + + + +complex double *qpms_zflip_uvswi_dense( + complex double *target, + const qpms_vswf_set_spec_t *bspec) +{ + target = ensure_alloc(target, bspec); + const size_t n = bspec->n; for (size_t row = 0; row < n; row++) { qpms_vswf_type_t rt; @@ -58,11 +69,8 @@ complex double *qpms_yflip_uvswi_dense( complex double *target, const qpms_vswf_set_spec_t *bspec) { - check_norm_compat(bspec); + target = ensure_alloc(target, bspec); const size_t n = bspec->n; - if (target == NULL) - target = malloc(n * n * sizeof(complex double)); - if (target == NULL) abort(); for (size_t row = 0; row < n; row++) { qpms_vswf_type_t rt; @@ -96,11 +104,8 @@ complex double *qpms_xflip_uvswi_dense( complex double *target, const qpms_vswf_set_spec_t *bspec) { - check_norm_compat(bspec); + target = ensure_alloc(target, bspec); const size_t n = bspec->n; - if (target == NULL) - target = malloc(n * n * sizeof(complex double)); - if (target == NULL) abort(); for (size_t row = 0; row < n; row++) { qpms_vswf_type_t rt; @@ -137,11 +142,8 @@ complex double *qpms_zrot_uvswi_dense( double phi ///< Rotation angle ) { - check_norm_compat(bspec); + target = ensure_alloc(target, bspec); const size_t n = bspec->n; - if (target == NULL) - target = malloc(n * n * sizeof(complex double)); - if (target == NULL) abort(); for (size_t row = 0; row < n; row++) { qpms_vswf_type_t rt; @@ -175,3 +177,34 @@ complex double *qpms_zrot_rational_uvswi_dense( return qpms_zrot_uvswi_dense(target, bspec, phi); } +complex double *qpms_irot3_uvswfi_dense( + complex double *target, + const qpms_vswf_set_spec_t *bspec, + const qpms_irot3_t t) +{ + target = ensure_alloc(target, bspec); + const size_t n = bspec->n; + + for (size_t row = 0; row < n; row++) { + qpms_vswf_type_t rt; + qpms_l_t rl; + qpms_m_t rm; + qpms_uvswfi2tmn(bspec->ilist[row], &rt, &rm, &rl); + for (size_t col = 0; col < n; col++) { + qpms_vswf_type_t ct; + qpms_l_t cl; + qpms_m_t cm; + if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); + if (rl == cl && rm == cm && rt == ct) + // TODO qpms_vswf_irot_elem_from_irot3 might be slow and not too accurate for large l + target[n*row + col] = // Checkme rm and cm order + qpms_vswf_irot_elem_from_irot3(t, + rl, rm /* CHECKME here */, cm /* and here */, + rt == QPMS_VSWF_MAGNETIC); + else target[n*row + col] = 0; + } + } + return target; +} + + diff --git a/qpms/symmetries.h b/qpms/symmetries.h index 8b64713..4548ff0 100644 --- a/qpms/symmetries.h +++ b/qpms/symmetries.h @@ -22,6 +22,7 @@ #ifndef SYMMETRIES_H #define SYMMETRIES_H #include "vswf.h" +#include "qpms_types.h" #include /// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring). @@ -52,4 +53,10 @@ complex double *qpms_zrot_rational_uvswi_dense( int w ); +/// Dense matrix (uvswi-indexed) representation of any O(3) transformation. +complex double *qpms_irot3_uvswfi_dense( + complex double *target, ///< If NULL, a new array is allocated. + const qpms_vswf_set_spec_t *bspec, + const qpms_irot3_t transf); + #endif // SYMMETRIES_H