Fix the integrand for axially symmetric T-matrices.

Forgotten jacobian, FFS.


Former-commit-id: f9a409bd1f5af78bb460d4bd3a9aadf78cebc2a7

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

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

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);
 }
 

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