Dense matrix uvswf repre of a general O(3) transformation.

Former-commit-id: 64cf762991956475bb1a3b208292262fa37f9365

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(
+// Used in the functions below to ensure memory allocation and checks for bspec validity
-    complex double *target,
+static inline complex double *ensure_alloc(complex double *target,
-    const qpms_vswf_set_spec_t *bspec) 
+    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;
+}
+
+ 

qpms/symmetries.h

@@ -22,6 +22,7 @@
 #ifndef SYMMETRIES_H
 #define SYMMETRIES_H
 #include "vswf.h"
+#include "qpms_types.h"
 #include <cblas.h>
 
 /// 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