From 324303478c7b46cba47fe0ece3b3bfa931a697ae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 17 Dec 2018 10:58:52 +0200
Subject: [PATCH] W file processing for variable particle number, some
 convenience funcs.

Former-commit-id: f75687d81961da077bce21fdfe4f1721173e00b7
---
 misc/processWfiles.py      |  7 ++++---
 misc/processWfiles2part.py | 10 ++++++++++
 qpms/qpms_p.py             | 12 ++++++++----
 qpms/tmatrices.py          | 13 +++++++++++++
 4 files changed, 35 insertions(+), 7 deletions(-)
 create mode 100755 misc/processWfiles2part.py

diff --git a/misc/processWfiles.py b/misc/processWfiles.py
index 9704ddc..eca9291 100755
--- a/misc/processWfiles.py
+++ b/misc/processWfiles.py
@@ -3,8 +3,9 @@
 import sys
 from qpms import processWfiles_sameKs
 
-dest = sys.argv[1]
-srcs = sys.argv[2:]
+npart = int(sys.argv[1])
+dest = sys.argv[2]
+srcs = sys.argv[3:]
 
-processWfiles_sameKs(srcs, dest, f='d')
+processWfiles_sameKs(srcs, dest, f='d', nparticles=npart)
 
diff --git a/misc/processWfiles2part.py b/misc/processWfiles2part.py
new file mode 100755
index 0000000..9704ddc
--- /dev/null
+++ b/misc/processWfiles2part.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+
+import sys
+from qpms import processWfiles_sameKs
+
+dest = sys.argv[1]
+srcs = sys.argv[2:]
+
+processWfiles_sameKs(srcs, dest, f='d')
+
diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py
index 6abb789..4fab6b9 100644
--- a/qpms/qpms_p.py
+++ b/qpms/qpms_p.py
@@ -742,13 +742,13 @@ def G0_sum_1_slow(source_cart, dest_cart, k, nmax):
     return G_Mie_scat_precalc_cart(source_cart, dest_cart, RH, RV, a=0.001, nmax=nmax, k_i=1, k_e=k, μ_i=1, μ_e=1, J_ext=1, J_scat=3)
 
 
-def loadWfile(fileName, lMax = None, fatForm = True):
+def loadWfile(fileName, lMax = None, fatForm = True, nparticles = 2):
     '''
     Load the translation operator lattice sums generated by hexlattice_ewald.c
     fatForm has different indices (so that it can be easily multiplied with T-Matrix
     and twice the size.
     '''
-    nparticles = 2 # TODO generalize
+    #nparticles = 2 # TODO generalize
     um = 1e-6 # micrometer in SI units
     data = np.loadtxt(fileName)
     nsamples = data.shape[0]
@@ -777,6 +777,7 @@ def loadWfile(fileName, lMax = None, fatForm = True):
     return {
             'lMax' : lMax,
             'nelem' : nelem,
+            'npart' : nparticles,
             'freqs' : freqs,
             'ks' : ks,
             'k0_effs' : k0s,
@@ -784,7 +785,7 @@ def loadWfile(fileName, lMax = None, fatForm = True):
             'freqs_weirdunits' : freqs_weirdunits,
             }
 
-def processWfiles_sameKs(freqfilenames, destfilename, lMax = None, f='npz'):
+def processWfiles_sameKs(freqfilenames, destfilename, nparticles = 2, lMax = None, f='npz'):
     '''
     Processes translation operator data in different files; each file is supposed to contain one frequency.
     The Ks in the different files are expected to be exactly the same and in the same order; 
@@ -794,7 +795,10 @@ def processWfiles_sameKs(freqfilenames, destfilename, lMax = None, f='npz'):
     'npz' creates npz archive, 'd' creates a directory with individual npy files, where the W-matrix data
     are writen using memmap.
     '''
-    nparticles = 2 #TODO generalize
+
+    if (nparticles is None):
+        # TODO here read nparticles from file;
+        raise # NOT IMPLEMENTED
     
     um = 1e-6 # micrometre in SI units
     
diff --git a/qpms/tmatrices.py b/qpms/tmatrices.py
index 03550d9..d58cf75 100644
--- a/qpms/tmatrices.py
+++ b/qpms/tmatrices.py
@@ -126,6 +126,19 @@ def zflip_tyty(lmax):
     fl_tyty[1,:,1,:] = -fl_yy
     return fl_tyty
 
+def zrotN_yy(N, lMax):
+    return WignerD_yy_fromvector(lMax, np.array([0,0,pi * (2/N)]))
+
+def op_yy2tyty(yyop):
+    '''
+    Broadcasts an yy operator to tyty operator without considering mirroring effects.
+    Good (maybe only) for rotations.
+    '''
+    return np.moveaxis(np.eye(2)[:,:,ň,ň] * yyop, 2,1)
+
+def zrotN_tyty(N, lMax):
+    return op_yy2tyty(zrotN_yy(N, lMax))
+
 def parity_yy(lmax):
     """
     Parity operator (flip in x,y,z)