From 03188d43f74b8c6c473ef49493684d5e740919d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Fri, 16 Aug 2019 10:05:23 +0300
Subject: [PATCH] Fix the integrand for axially symmetric T-matrices.

Forgotten jacobian, FFS.


Former-commit-id: f9a409bd1f5af78bb460d4bd3a9aadf78cebc2a7
---
 qpms/cymaterials.pyx             |  2 +-
 qpms/qpms_c.pyx                  |  1 -
 qpms/tmatrices.c                 |  2 +-
 tests/cytmatrices_consistence.py | 42 ++++++++++++++++++++++++++++++++
 4 files changed, 44 insertions(+), 3 deletions(-)
 create mode 100644 tests/cytmatrices_consistence.py

diff --git a/qpms/cymaterials.pyx b/qpms/cymaterials.pyx
index 3b1b7f3..55068b7 100644
--- a/qpms/cymaterials.pyx
+++ b/qpms/cymaterials.pyx
@@ -105,7 +105,7 @@ cdef class EpsMuGenerator:
             self.g.function = qpms_permittivity_interpolator_epsmu_g
             self.g.params = (<MaterialInterpolator?>self.holder).rawpointer()
         elif isinstance(what, EpsMuGenerator): # Copy constructor
-            self.holder = what.holder
+            self.holder = (<EpsMuGenerator?>what).holder
             self.g = (<EpsMuGenerator?>what).g
         elif callable(what):
             warnings.warn("Custom python (eps,mu) generators are an experimental feature")
diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index 9660c1e..2292f08 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -1,5 +1,4 @@
 """@package qpms_c
-            self.s.norm = <qpms_normalisation_t>(QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE)
 Cythonized parts of QPMS; mostly wrappers over the C data structures
 to make them available in Python.
 """
diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c
index 95d6298..b4522cf 100644
--- a/qpms/tmatrices.c
+++ b/qpms/tmatrices.c
@@ -687,7 +687,7 @@ static double tmatrix_axialsym_integrand(double theta, void *param) {
   complex double tp2 = nrc * (y2.thetac * v_in2.phic - y2.phic * v_in2.thetac)
                  + nthetac * (y2.phic   * v_in2.rc   - y2.rc   * v_in2.phic);
   double jac = SQ(rb.r) * sin(theta) / nrc; // Jacobian
-  complex double res = p->z/p->z_in * tp1 + tp2;
+  complex double res = jac * (p->z/p->z_in * tp1 + tp2);
   return p->realpart ? creal(res) : cimag(res);
 }
 
diff --git a/tests/cytmatrices_consistence.py b/tests/cytmatrices_consistence.py
new file mode 100644
index 0000000..4d0a689
--- /dev/null
+++ b/tests/cytmatrices_consistence.py
@@ -0,0 +1,42 @@
+from qpms import *
+import numpy as np
+R = 40e-9
+lMax = 2
+ω = 1.5*eV/ℏ
+
+spec = BaseSpec(lMax = lMax)
+
+inside = EpsMuGenerator(lorentz_drude['Au'])
+outside = EpsMuGenerator(EpsMu(2.3104,1))
+
+gensphere_arc = TMatrixGenerator.sphere_asarc(outside=outside, inside=inside, r=R, lMax_extend=lMax)
+
+np.set_printoptions(precision=3, suppress=True,linewidth=1000)
+
+QT = gensphere_arc.Q_transposed(ω, spec.norm)
+RT = gensphere_arc.R_transposed(ω, spec.norm)
+T = gensphere_arc(spec,ω)
+
+
+QT_corrected = np.array(QT)
+RT_corrected = np.array(RT)
+QT_corrected[8:,:8] = 0
+QT_corrected[:8,8:] = 0
+RT_corrected[8:,:8] = 0
+RT_corrected[:8,8:] = 0
+T_corrected = np.dot(np.linalg.inv(QT_corrected), RT_corrected)
+
+print("QT:")
+print(QT)
+
+print("RT:")
+print(RT)
+
+print("T:")
+print(T[:])
+
+print("T_corrected:")
+print(T_corrected)
+
+
+