qpms/qpms/symmetries.c

252 lines
7.2 KiB
C

#include "symmetries.h"
#include "tiny_inlines.h"
#include "indexing.h"
#include "quaternions.h"
#include "qpms_error.h"
// TODO at some point, maybe support also other norms.
// (perhaps just use qpms_normalisation_t_factor() at the right places)
static inline void check_norm_compat(const qpms_vswf_set_spec_t *s)
{
switch (s->norm & QPMS_NORMALISATION_NORM_BITS) {
case QPMS_NORMALISATION_NORM_POWER:
break;
case QPMS_NORMALISATION_NORM_SPHARM:
break;
default:
QPMS_NOT_IMPLEMENTED("At the moment, only spherical harmonics of spherical harmonics or power normalisations implemented.");
}
}
static inline void ONLY_EIMF_IMPLEMENTED(const qpms_normalisation_t norm)
{
if (norm & QPMS_NORMALISATION_SPHARM_REAL)
QPMS_NOT_IMPLEMENTED("Support for real spherical harmonics not implemented yet.");
}
// 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)
QPMS_CRASHING_MALLOC(target, n * n * sizeof(complex double));
return target;
}
complex double *qpms_zflip_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;
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)
switch(rt) {
case QPMS_VSWF_ELECTRIC:
case QPMS_VSWF_LONGITUDINAL:
target[n*row + col] = min1pow(cm + cl);
break;
case QPMS_VSWF_MAGNETIC:
target[n*row + col] = -min1pow(cm + cl);
break;
default:
abort();
}
else target[n*row + col] = 0;
}
}
return target;
}
complex double *qpms_yflip_uvswi_dense(
complex double *target,
const qpms_vswf_set_spec_t *bspec)
{
check_norm_compat(bspec);
ONLY_EIMF_IMPLEMENTED(bspec->norm);
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)
switch(rt) {
case QPMS_VSWF_ELECTRIC:
case QPMS_VSWF_LONGITUDINAL:
target[n*row + col] = min1pow(rm);
break;
case QPMS_VSWF_MAGNETIC:
target[n*row + col] = -min1pow(rm);
break;
default:
abort();
}
else target[n*row + col] = 0;
}
}
return target;
}
complex double *qpms_xflip_uvswi_dense(
complex double *target,
const qpms_vswf_set_spec_t *bspec)
{
check_norm_compat(bspec);
ONLY_EIMF_IMPLEMENTED(bspec->norm);
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)
switch(rt) {
case QPMS_VSWF_ELECTRIC:
case QPMS_VSWF_LONGITUDINAL:
target[n*row + col] = 1;
break;
case QPMS_VSWF_MAGNETIC:
target[n*row + col] = -1;
break;
default:
abort();
}
else target[n*row + col] = 0;
}
}
return target;
}
// Dense matrix representation of a rotation around the z-axis
complex double *qpms_zrot_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
double phi ///< Rotation angle
)
{
QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well.
check_norm_compat(bspec);
ONLY_EIMF_IMPLEMENTED(bspec->norm);
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 COMPARE WITH PYTHON
target[n*row + col] = cexp(/* - ?*/I * rm * phi);
else target[n*row + col] = 0;
}
}
return target;
}
// Dense matrix representation of a "rational" rotation around the z-axis
/* Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
*/
complex double *qpms_zrot_rational_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
int N,
int w
)
{
double phi = 2 * M_PI * w / N;
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)
{
QPMS_UNTESTED; // not sure about the C.-S. phase. Don't forget documenting it as well.
check_norm_compat(bspec);
ONLY_EIMF_IMPLEMENTED(bspec->norm);
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 && 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;
}
size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol) {
size_t changed = 0;
for(size_t i = 0; i < nmemb; ++i)
if(fabs(arr[i]) <= atol) {
arr[i] = 0;
++changed;
}
return changed;
}
size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol) {
size_t changed = 0;
for(size_t i = 0; i < nmemb; ++i) {
if(fabs(creal(arr[i])) <= atol) {
arr[i] = I*cimag(arr[i]);
++changed;
}
if(fabs(cimag(arr[i])) <= atol) {
arr[i] = creal(arr[i]);
++changed;
}
}
}