From 147bed193f6691bbd5f888855d705d3ef56a9e14 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 25 Feb 2019 18:39:16 +0200
Subject: [PATCH] Dense matrix uvswf repre of a general O(3) transformation.

Former-commit-id: 64cf762991956475bb1a3b208292262fa37f9365
---
 qpms/symmetries.c | 65 +++++++++++++++++++++++++++++++++++------------
 qpms/symmetries.h |  7 +++++
 2 files changed, 56 insertions(+), 16 deletions(-)

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 <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