oldtests/test_qpms_p.py

   1 """
   2 Unit tests for qpms_p
   3 =====================
   4
   5 Covered functions
   6 -----------------
   7 plane_pq_y vs. vswf_yr1
   8
   9 Not covered
  10 -----------
  11 Everything else
  12
  13 """
  14
  15 import unittest
  16 import qpms
  17 import numpy as np
  18 from numpy import newaxis as ň
  19 import warnings
  20
  21 # Some constants go here.
  22 lengthOrdersOfMagnitude = [2.**i for i in range(-15,10,2)]
  23
  24 class PlaneWaveDecompositionTests(unittest.TestCase):
  25     """
  26     covers plane_pq_y and vswf_yr1
  27     """
  28     def testRandomInc(self):
  29         # The "maximum" argument of the Bessel's functions, i.e. maximum wave number times the distance,
  30         # for the "locally strongly varying fields"
  31         maxx = 10
  32         rfailtol = 0.01 # how much of the randomized test will be tolerated
  33         lMax = 80 # To which order we decompose the waves
  34         rtol = 1e-5 # relative required precision
  35         atol = 1. # absolute tolerance, does not really play a role
  36         nsamples = 4 # (frequency, direction, polarisation) samples per order of magnitude and test
  37         npoints = 15 # points to evaluate per sample
  38
  39
  40         failcounter = 0
  41         passcounter = 0
  42         for oom in lengthOrdersOfMagnitude:
  43             k = np.random.randn(nsamples, 3) / oom
  44             ksiz = np.linalg.norm(k, axis=-1)
  45             kdir = k / ksiz[...,ň]
  46             E_0 = np.cross(np.random.randn(nsamples, 3), k) * oom # ensure orthogonality
  47             for s in range(nsamples):
  48                 testpoints = oom * maxx * np.random.randn(npoints, 3)
  49                 p, q = qpms.plane_pq_y(lMax, k[s], E_0[s])
  50                 planewave_1 =  np.exp(1j*np.dot(testpoints,k[s]))[:,ň] * E_0[s,:]
  51                 for i in range(npoints):
  52                     sph = qpms.cart2sph(ksiz[s]*testpoints[i])
  53                     M̃_y, Ñ_y = qpms.vswf_yr1(sph, lMax, 1)
  54                     planewave_2_p = -1j*qpms.sph_loccart2cart(np.dot(p,Ñ_y)+np.dot(q,M̃_y),sph)
  55                     #self.assertTrue(np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol))
  56                     if not np.allclose(planewave_2_p, planewave_1[i], rtol=rtol, atol=atol):
  57                         False and warnings.warn('Planewave expansion test not passed; r = '
  58                                 +str(testpoints[i])+', k = '+str(k[s])
  59                                 +', E_0 = '+str(E_0[s])+', (original) E = '
  60                                 +str(planewave_1[i])+', (reexpanded) E = '
  61                                 +str(planewave_2_p)
  62                                 +', x = '+str(np.dot(testpoints[i],k[s]))
  63                                 +'; distance = '
  64                                 +str(np.linalg.norm(planewave_1[i]-planewave_2_p))
  65                                 +', required relative precision = '
  66                                 +str(rtol)+'.')
  67                         failcounter += 1
  68                     else:
  69                         passcounter += 1
  70         self.assertLess(failcounter / (failcounter + passcounter), rfailtol,
  71                 '%d / %d (%.2e) randomized numerical tests failed (tolerance %.2e)'
  72                 % (failcounter, failcounter + passcounter,
  73                     failcounter / (failcounter + passcounter), rfailtol))
  74         return
  75
  76     def testCornerCases(self):
  77         pass
  78
  79
  80 class SphericalWaveTranslationTests(unittest.TestCase):
  81     def testRandom1to1(self):
  82         # The "maximum" argument of the Bessel's functions, i.e. maximum wave number times the distance,
  83         # for the "locally strongly varying fields"
  84         maxx = 10
  85         rfailtol = 0.01 # how much of the randomized test fail proportion will be tolerated
  86         lMax = 50 # To which order we decompose the waves
  87         lMax_outgoing = 4 # To which order we try the outgoing waves
  88         rtol = 1e-5 # relative required precision
  89         atol = 1. # absolute tolerance, does not really play a role
  90         nsamples = 4 # frequency samples per order of magnitude and test
  91         npoints = 15 # points to evaluate per frequency and center
  92
  93         ncentres = 3 # number of spherical coordinate centres between which the translations are to be made
  94         maxxd = 2000 # the center position standard deviation
  95
  96         failcounter = 0
  97         passcounter = 0
  98         my, ny = qpms.get_mn_y(lMax)
  99         nelem_full = len(my)
 100         nelem_out = lMax_outgoing * (lMax_outgoing + 2)
 101         for oom in lengthOrdersOfMagnitude:
 102             centres = oom * maxxd * np.random.randn(ncentres, 3)
 103             ksizs = np.random.randn(nsamples)
 104             for ksiz in ksizs:
 105                 for i in range(ncentres): # "source"
 106                     Rs = centres[i]
 107                     testr = oom * maxx * np.random.randn(npoints, 3)
 108                     for j in range(ncentres): # "destination"
 109                         if j == i:
 110                             continue
 111                         Rd = centres[j]
 112                         shift = Rd - Rs
 113                         shift_sph = qpms.cart2sph(shift)
 114                         shift_kr = ksiz * shift_sph[0]
 115                         shift_theta = shift_sph[1]
 116                         shift_phi = shift_sph[2]
 117
 118                         A_yd_ys = np.empty((nelem_full,nelem_out), dtype = np.complex_)
 119                         B_yd_ys = np.empty((nelem_full,nelem_out), dtype = np.complex_)
 120                         for yd in range(nelem_full):
 121                             for ys in range(nelem_out):
 122                                 A_yd_ys[yd, ys] = qpms.Ã(my[yd],ny[yd],my[ys],ny[ys],shift_kr, shift_theta, shift_phi, True, 1)
 123                                 B_yd_ys[yd, ys] = qpms.B̃(my[yd],ny[yd],my[ys],ny[ys],shift_kr, shift_phi, shift_phi, True, 1)
 124                         for r in testr:
 125                             sph_ssys = qpms.cart2sph(r+Rd-Rs)
 126                             M_ssys, N_ssys = qpms.vswf_yr1(np.array([ksiz * sph_ssys[0], sph_ssys[1], sph_ssys[2]]), lMax_outgoing, J=1)
 127                             sph_dsys = qpms.cart2sph(r)
 128                             M_dsys, N_dsys = qpms.vswf_yr1(np.array([ksiz * sph_dsys[0], sph_dsys[1], sph_dsys[2]]), lMax, J=1)
 129                             for ys in range(nelem_out):
 130                                 # Electrical waves
 131                                 E_1 = -1j*qpms.sph_loccart2cart(N_ssys[ys], sph_ssys)
 132                                 E_2 = -1j*qpms.sph_loccart2cart(np.dot(A_yd_ys[:,ys],N_dsys)+np.dot(B_yd_ys[:,ys],M_dsys),sph_dsys)
 133                                 if not np.allclose(E_1, E_2, rtol=rtol, atol=atol):
 134                                     failcounter += 1
 135                                 else:
 136                                     passcounter += 1
 137                                 # Magnetic waves
 138                                 E_1 = -1j*qpms.sph_loccart2cart(M_ssys[ys], sph_ssys)
 139                                 E_2 = -1j*qpms.sph_loccart2cart(np.dot(A_yd_ys[:,ys],M_dsys)+np.dot(B_yd_ys[:,ys],N_dsys),sph_dsys)
 140                                 if not np.allclose(E_1, E_2, rtol=rtol, atol=atol):
 141                                     failcounter += 1
 142                                 else:
 143                                     passcounter += 1
 144         self.assertLess(failcounter / (failcounter + passcounter), rfailtol,
 145                 '%d / %d (%.2e) randomized numerical tests failed (tolerance %.2e)'
 146                 % (failcounter, failcounter + passcounter,
 147                     failcounter / (failcounter + passcounter), rfailtol))
 148         return
 149
 150     def testRandom3to1(self):
 151         pass
 152
 153 def main():
 154     unittest.main()
 155
 156 if __name__ == '__main__':
 157     main()
 158