From 86b87c955cf51b3af418a01167c9d0b35ace1b6b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 17 Mar 2019 23:32:51 +0000
Subject: [PATCH] Fix the alternative order modematrix implementation

Former-commit-id: 7daaf636ca803565e2a72b03e1ce184ece24e4c0
---
 qpms/scatsystem.c                    |  4 +--
 tests/modeproblem_implementations.py | 44 ++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 2 deletions(-)
 create mode 100644 tests/modeproblem_implementations.py

diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index 4878506..69acca6 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -1472,7 +1472,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
   const complex double one = 1, zero = 0;
 
   for(qpms_ss_pi_t opistartR = 0; opistartR < ss->p_count; 
-      opistartR += ss->orbit_types[ss->p_orbitinfo[opistartR].t].size //orbit_p_countR; might write a while() instead
+      opistartR += ss->orbit_types[ss->p_orbitinfo[ss->p_by_orbit[opistartR]].t].size //orbit_p_countR; might write a while() instead
      ) {
     const qpms_ss_pi_t orbitstartpiR = ss->p_by_orbit[opistartR];
     const qpms_ss_oti_t otiR = ss->p_orbitinfo[orbitstartpiR].t;
@@ -1495,7 +1495,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
       for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) {
         qpms_ss_pi_t piR = ss->p_by_orbit[opistartR + opiR];
         assert(opiR == ss->p_orbitinfo[piR].p);
-        const qpms_ss_oti_t otiR = ss->p_orbitinfo[piR].t;
+        assert(otiR == ss->p_orbitinfo[piR].t);
         assert(ss->p_orbitinfo[piR].osn == osnR);
         const cart3_t posR = ss->p[piR].pos;
         // dest particle T-matrix
diff --git a/tests/modeproblem_implementations.py b/tests/modeproblem_implementations.py
new file mode 100644
index 0000000..bc2de15
--- /dev/null
+++ b/tests/modeproblem_implementations.py
@@ -0,0 +1,44 @@
+from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c
+from qpms.symmetries import point_group_info
+import numpy as np
+import os
+import sys
+nm = 1e-9
+
+sym = FinitePointGroup(point_group_info['D2h'])
+bspec = BaseSpec(lMax = 2)
+#tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
+tmfile = '/home/mmn/repo/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
+#outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_new'
+outputdatadir = '/tmp/u/46/necadam1/unix/project/AaroBECfinite_new'
+os.makedirs(outputdatadir, exist_ok = True)
+interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
+# There is only one t-matrix in the system for each frequency. We initialize the matrix with the lowest frequency data.
+# Later, we can replace it using the tmatrix[...] = interp(freq) and s.update_tmatrices NOT YET; TODO
+
+omega = 1.475 * eV/hbar
+sv_threshold = 0.5
+
+# Now place the particles and set background index.
+px = 571*nm; py = 621*nm
+n = 1.51
+Nx = 5
+Ny = 7
+
+orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
+orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
+
+orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
+
+tmatrix = interp(omega)
+particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
+
+
+ss = ScatteringSystem(particles, sym)
+k = n * omega / c
+
+for iri in range(ss.nirreps):
+    mm_iri_orig = ss.modeproblem_matrix_packed(k, iri)
+    mm_iri_alt = ss.modeproblem_matrix_packed(k, iri, version='R')
+    print(np.amax(abs(mm_iri_orig-mm_iri_alt)))
+