From 8fb938961f5fd94c36fe177f8b5a33f3dc29b61c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 9 Aug 2016 14:39:45 +0300 Subject: [PATCH] Correct flips for magnetic part. Some tests. Former-commit-id: 998f4b8fa5b94447ae75220fbb43c0143d4083a5 --- qpms/qpms_p.py | 39 ++++++++++++- tests/test_qpms_p.py | 131 ++++++++++++++++++++++++++++++++++--------- 2 files changed, 141 insertions(+), 29 deletions(-) diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py index ea999b2..2205d9c 100644 --- a/qpms/qpms_p.py +++ b/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") diff --git a/tests/test_qpms_p.py b/tests/test_qpms_p.py index 63a422d..427594e 100644 --- a/tests/test_qpms_p.py +++ b/tests/test_qpms_p.py @@ -19,26 +19,26 @@ from numpy import newaxis as ň import warnings # Some constants go here. - -# 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 +lengthOrdersOfMagnitude = [2.**i for i in range(-15,10,2)] 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)) -# if not np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol): -# warnings.warn('Planewave expansion test not passed; r = ' -# +str(testpoints[i])+', k = '+str(k[s]) -# +', E_0 = '+str(E_0[s])+', (original) E = ' -# +str(planewave_1[i])+', (reexpanded) E = ' -# +str(planewave_2_p) -# +', x = '+str(np.dot(testpoints[i],k[s])) -# +'; distance = ' -# +str(np.linalg.norm(planewave_1[i]-planewave_2_p)) -# +', required relative precision = ' -# +str(relprecision)+'.') - return + #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): + False and warnings.warn('Planewave expansion test not passed; r = ' + +str(testpoints[i])+', k = '+str(k[s]) + +', E_0 = '+str(E_0[s])+', (original) E = ' + +str(planewave_1[i])+', (reexpanded) E = ' + +str(planewave_2_p) + +', x = '+str(np.dot(testpoints[i],k[s])) + +'; distance = ' + +str(np.linalg.norm(planewave_1[i]-planewave_2_p)) + +', required relative precision = ' + +str(rtol)+'.') + 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():