Correct flips for magnetic part. Some tests.

Former-commit-id: 998f4b8fa5b94447ae75220fbb43c0143d4083a5

qpms/qpms_p.py

@@ -915,7 +915,19 @@ def xflip_yy(lmax):
         elems[b_in:e_in,b_in:e_in] = np.eye(2*l+1)[::-1,:]
         b_in = e_in
     return elems
-    
+
+def xflip_tyy(lmax):
+    fl_yy = xflip_yy(lmax)
+    return np.array([fl_yy,-fl_yy])
+
+def xflip_tyty(lmax):
+    fl_yy = xflip_yy(lmax)
+    nelem = fl_yy.shape[0]
+    fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
+    fl_tyty[0,:,0,:] = fl_yy
+    fl_tyty[1,:,1,:] = -fl_yy
+    return fl_tyty
+
 def yflip_yy(lmax):
     """
     TODO doc
@@ -927,6 +939,18 @@ def yflip_yy(lmax):
     elems[(my % 2)==1] = elems[(my % 2)==1] * -1 # Obvious sign of tiredness (this is correct but ugly; FIXME)
     return elems
 
+def yflip_tyy(lmax):
+    fl_yy = yflip_yy(lmax)
+    return np.array([fl_yy,-fl_yy])
+
+def yflip_tyty(lmax):
+    fl_yy = yflip_yy(lmax)
+    nelem = fl_yy.shape[0]
+    fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
+    fl_tyty[0,:,0,:] = fl_yy
+    fl_tyty[1,:,1,:] = -fl_yy
+    return fl_tyty
+
 def zflip_yy(lmax):
     """
     TODO doc
@@ -942,6 +966,18 @@ def zflip_yy(lmax):
         b_in = e_in
     return elems
 
+def zflip_tyy(lmax):
+    fl_yy = zflip_yy(lmax)
+    return np.array([fl_yy,-fl_yy])
+
+def zflip_tyty(lmax):
+    fl_yy = zflip_yy(lmax)
+    nelem = fl_yy.shape[0]
+    fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
+    fl_tyty[0,:,0,:] = fl_yy
+    fl_tyty[1,:,1,:] = -fl_yy
+    return fl_tyty
+
 def parity_yy(lmax):
     """
     Parity operator (flip in x,y,z)
@@ -1294,6 +1330,7 @@ def scatter_plane_wave_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, k_
         return ab
     returnlist = [ab]
     if (return_pq_0):
+        pq_0_arr.shape = ab.shape
         returnlist.append(pq_0_arr)
     if (return_pq):
         warnings.warn("return_pq not implemented, ignoring")

tests/test_qpms_p.py

@@ -19,26 +19,26 @@ from numpy import newaxis as ň
 import warnings
 
 # Some constants go here.
-
+lengthOrdersOfMagnitude = [2.**i for i in range(-15,10,2)]
-# The "maximum" argument of the Bessel's functions, i.e. maximum wave number times the distance,
-# for the "locally strongly varying fields"
-maxx = 3
-# The "maximum" argument of the Bessel's function for reexpansion of the waves into regular waves
-# in another center
-maxxd = 2000
-lMax = 50 # To which order we decompose the waves
-
-lengthOrdersOfMagnitude = [2.**i for i in range(-15,10)]
-nsamples = 4 # (frequency, direction, polarisation) samples per order of magnitude and test
-npoints = 40 # points to evaluate per sample
-rtol = 1e-7 # relative required precision
-atol = 1. # absolute tolerance, does not really play a role
 
 class PlaneWaveDecompositionTests(unittest.TestCase):
     """
     covers plane_pq_y and vswf_yr1
     """
     def testRandomInc(self):
+        # The "maximum" argument of the Bessel's functions, i.e. maximum wave number times the distance,
+        # for the "locally strongly varying fields"
+        maxx = 10
+        rfailtol = 0.01 # how much of the randomized test will be tolerated
+        lMax = 80 # To which order we decompose the waves
+        rtol = 1e-5 # relative required precision
+        atol = 1. # absolute tolerance, does not really play a role
+        nsamples = 4 # (frequency, direction, polarisation) samples per order of magnitude and test
+        npoints = 15 # points to evaluate per sample
+
+    
+        failcounter = 0
+        passcounter = 0
         for oom in lengthOrdersOfMagnitude:
             k = np.random.randn(nsamples, 3) / oom
             ksiz = np.linalg.norm(k, axis=-1)
@@ -52,27 +52,102 @@ class PlaneWaveDecompositionTests(unittest.TestCase):
                     sph = qpms.cart2sph(ksiz[s]*testpoints[i])
                     M̃_y, Ñ_y = qpms.vswf_yr1(sph, lMax, 1)
                     planewave_2_p = -1j*qpms.sph_loccart2cart(np.dot(p,Ñ_y)+np.dot(q,M̃_y),sph)
-                    self.assertTrue(np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol))
+                    #self.assertTrue(np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol))
-#                    if not np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol):
+                    if not np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol):
-#                        warnings.warn('Planewave expansion test not passed; r = '
+                        False and warnings.warn('Planewave expansion test not passed; r = '
-#                                +str(testpoints[i])+', k = '+str(k[s])
+                                +str(testpoints[i])+', k = '+str(k[s])
-#                                +', E_0 = '+str(E_0[s])+', (original) E = '
+                                +', E_0 = '+str(E_0[s])+', (original) E = '
-#                                +str(planewave_1[i])+', (reexpanded) E = '
+                                +str(planewave_1[i])+', (reexpanded) E = '
-#                                +str(planewave_2_p)
+                                +str(planewave_2_p)
-#                                +', x = '+str(np.dot(testpoints[i],k[s]))
+                                +', x = '+str(np.dot(testpoints[i],k[s]))
-#                                +'; distance = '
+                                +'; distance = '
-#                                +str(np.linalg.norm(planewave_1[i]-planewave_2_p))
+                                +str(np.linalg.norm(planewave_1[i]-planewave_2_p))
-#                                +', required relative precision = '
+                                +', required relative precision = '
-#                                +str(relprecision)+'.')
+                                +str(rtol)+'.')
-        return
+                        failcounter += 1
+                    else:
+                        passcounter += 1
+        self.assertLess(failcounter / (failcounter + passcounter), rfailtol, 
+                '%d / %d (%.2e) randomized numerical tests failed (tolerance %.2e)' 
+                % (failcounter, failcounter + passcounter, 
+                    failcounter / (failcounter + passcounter), rfailtol))
+        return 
 
     def testCornerCases(self):
         pass
 
 
 class SphericalWaveTranslationTests(unittest.TestCase):
+    def testRandom1to1(self):
+        # The "maximum" argument of the Bessel's functions, i.e. maximum wave number times the distance,
+        # for the "locally strongly varying fields"
+        maxx = 10
+        rfailtol = 0.01 # how much of the randomized test fail proportion will be tolerated
+        lMax = 50 # To which order we decompose the waves
+        lMax_outgoing = 4 # To which order we try the outgoing waves
+        rtol = 1e-5 # relative required precision
+        atol = 1. # absolute tolerance, does not really play a role
+        nsamples = 4 # frequency samples per order of magnitude and test
+        npoints = 15 # points to evaluate per frequency and center
+     
+        ncentres = 3 # number of spherical coordinate centres between which the translations are to be made
+        maxxd = 2000 # the center position standard deviation
 
-    def sometest(self):
+        failcounter = 0
+        passcounter = 0
+        my, ny = qpms.get_mn_y(lMax)
+        nelem_full = len(my)
+        nelem_out = lMax_outgoing * (lMax_outgoing + 2)
+        for oom in lengthOrdersOfMagnitude:
+            centres = oom * maxxd * np.random.randn(ncentres, 3)
+            ksizs = np.random.randn(nsamples)
+            for ksiz in ksizs:
+                for i in range(ncentres): # "source"
+                    Rs = centres[i]
+                    testr = oom * maxx * np.random.randn(npoints, 3)
+                    for j in range(ncentres): # "destination"
+                        if j == i:
+                            continue
+                        Rd = centres[j]
+                        shift = Rd - Rs
+                        shift_sph = qpms.cart2sph(shift)
+                        shift_kr = ksiz * shift_sph[0]
+                        shift_theta = shift_sph[1]
+                        shift_phi = shift_sph[2]
+
+                        A_yd_ys = np.empty((nelem_full,nelem_out), dtype = np.complex_)
+                        B_yd_ys = np.empty((nelem_full,nelem_out), dtype = np.complex_)
+                        for yd in range(nelem_full):
+                            for ys in range(nelem_out):
+                                A_yd_ys[yd, ys] = qpms.Ã(my[yd],ny[yd],my[ys],ny[ys],shift_kr, shift_theta, shift_phi, True, 1) 
+                                B_yd_ys[yd, ys] = qpms.B̃(my[yd],ny[yd],my[ys],ny[ys],shift_kr, shift_phi, shift_phi, True, 1)
+                        for r in testr:
+                            sph_ssys = qpms.cart2sph(r+Rd-Rs)
+                            M_ssys, N_ssys = qpms.vswf_yr1(np.array([ksiz * sph_ssys[0], sph_ssys[1], sph_ssys[2]]), lMax_outgoing, J=1)
+                            sph_dsys = qpms.cart2sph(r)
+                            M_dsys, N_dsys = qpms.vswf_yr1(np.array([ksiz * sph_dsys[0], sph_dsys[1], sph_dsys[2]]), lMax, J=1)
+                            for ys in range(nelem_out): 
+                                # Electrical waves
+                                E_1 = -1j*qpms.sph_loccart2cart(N_ssys[ys], sph_ssys)
+                                E_2 = -1j*qpms.sph_loccart2cart(np.dot(A_yd_ys[:,ys],N_dsys)+np.dot(B_yd_ys[:,ys],M_dsys),sph_dsys)
+                                if not np.allclose(E_1, E_2, rtol=rtol, atol=atol):
+                                    failcounter += 1
+                                else:
+                                    passcounter += 1
+                                # Magnetic waves
+                                E_1 = -1j*qpms.sph_loccart2cart(M_ssys[ys], sph_ssys)
+                                E_2 = -1j*qpms.sph_loccart2cart(np.dot(A_yd_ys[:,ys],M_dsys)+np.dot(B_yd_ys[:,ys],N_dsys),sph_dsys)
+                                if not np.allclose(E_1, E_2, rtol=rtol, atol=atol):
+                                    failcounter += 1
+                                else:
+                                    passcounter += 1
+        self.assertLess(failcounter / (failcounter + passcounter), rfailtol, 
+                '%d / %d (%.2e) randomized numerical tests failed (tolerance %.2e)' 
+                % (failcounter, failcounter + passcounter, 
+                    failcounter / (failcounter + passcounter), rfailtol))
+        return 
+
+    def testRandom3to1(self):
         pass
 
 def main():