252 lines
7.2 KiB
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;
|
|
}
|
|
}
|
|
}
|
|
|