Plane wave test, docs.

Former-commit-id: 86bdabea5c1a424401d3b93ea98af16bd2e963be

README.rst

@@ -14,3 +14,15 @@ After all dependencies are installed, install qpms to your local python library
 
 python3 setup.py install --user
 
+
+Easiest installation ever 
+=========================
+(Just skip those you have already installed.)
+
+pip3 install --user numpy
+pip3 install --user scipy
+pip3 install --user cython
+pip3 install --user git+https://github.com/moble/quaternion.git
+pip3 install --user git+https://github.com/moble/spherical_functions.git
+pip3 install --user git+https://github.com/texnokrates/py_gmm.git@standalone_mie
+python3 setup.py install --user

qpms/qpms_p.py

@@ -10,8 +10,10 @@ import math
 import cmath
 import quaternion, spherical_functions as sf # because of the Wigner matrices
 import sys, time
+from numba import jit
 
 # Coordinate transforms for arrays of "arbitrary" shape
+#@jit
 def cart2sph(cart,axis=-1):
     if (cart.shape[axis] != 3):
         raise ValueError("The converted array has to have dimension 3"
@@ -23,6 +25,7 @@ def cart2sph(cart,axis=-1):
     φ = np.arctan2(y,x) # arctan2 handles zeroes correctly itself
     return np.concatenate((r,θ,φ),axis=axis)
 
+#@jit
 def sph2cart(sph, axis=-1):
     if (sph.shape[axis] != 3):
         raise ValueError("The converted array has to have dimension 3"
@@ -34,6 +37,7 @@ def sph2cart(sph, axis=-1):
     z = r * np.cos(θ)
     return np.concatenate((x,y,z),axis=axis)
 
+#@jit
 def sph_loccart2cart(loccart, sph, axis=-1):
     """
     Transformation of vector specified in local orthogonal coordinates 
@@ -83,6 +87,7 @@ def sph_loccart2cart(loccart, sph, axis=-1):
     out=inr̂*r̂+inθ̂*θ̂+inφ̂*φ̂
     return out
 
+#@jit
 def sph_loccart_basis(sph, sphaxis=-1, cartaxis=None):
     """
     Returns the local cartesian basis in terms of global cartesian basis.
@@ -118,6 +123,7 @@ def sph_loccart_basis(sph, sphaxis=-1, cartaxis=None):
     out = np.concatenate((x,y,z),axis=cartaxis)
     return out
 
+@jit
 def lpy(nmax, z):
     """
     Associated legendre function and its derivatative at z in the 'y-indexing'.
@@ -194,14 +200,18 @@ def vswf_yr(pos_sph,nmax,J=1):
     pass
 
 from scipy.special import sph_jn, sph_yn
+@jit
 def _sph_zn_1(n,z):
     return sph_jn(n,z)
+@jit
 def _sph_zn_2(n,z):
     return sph_yn(n,z)
+@jit
 def _sph_zn_3(n,z):
     besj=sph_jn(n,z)
     besy=sph_yn(n,z)
     return (besj[0] + 1j*besy[0],besj[1] + 1j*besy[1])
+@jit
 def _sph_zn_4(n,z):
     besj=sph_jn(n,z)
     besy=sph_yn(n,z)
@@ -294,6 +304,7 @@ def get_π̃τ̃_y1(θ,nmax):
     τ̃_y = prenorm * dPy * (- math.sin(θ))  # TADY BACHA!!!!!!!!!! * (- math.sin(pos_sph[1])) ???
     return (π̃_y,τ̃_y)
     
+#@jit
 def vswf_yr1(pos_sph,nmax,J=1):
     """
     As vswf_yr, but evaluated only at single position (i.e. pos_sph has
@@ -376,6 +387,27 @@ def plane_pq_y(nmax, kdir_cart, E_cart):
     the phase and polarisation state). If E_cart and kdir_cart are
     not orthogonal, the result should correspond to the k-normal part
     of E_cart.
+
+    The Taylor's convention on the expansion is used; therefore,
+    the electric field is expressed as
+    E(r) = E_cart * exp(1j * k ∙ r)
+         = -1j * ∑_y ( p_y * Ñ_y(|k| * r) + q_y * M̃_y(|k| * r))
+    (note the -1j factor).
+
+    Parameters
+    ----------
+    nmax: int
+        The order of the expansion.
+    kdir_cart: (3,)-shaped real array
+        The wave vector (its magnitude does not play a role).
+    E_cart: (3,)-shaped complex array
+        Amplitude of the plane wave
+
+    Returns
+    -------
+    p_y, q_y:
+        The expansion coefficients for the electric (Ñ) and magnetic
+        (M̃) waves, respectively.
     """
     if np.iscomplexobj(kdir_cart):
         warnings.warn("The direction vector for the plane wave coefficients should be real. I am discarding the imaginary part now.")
@@ -419,6 +451,7 @@ def a_q(m,n,μ,ν,q = None):
 
 # All arguments are single numbers (for now)
 # ZDE VYCHÁZEJÍ DIVNÁ ZNAMÉNKA
+#@jit
 def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
     exponent=(math.lgamma(2*n+1)-math.lgamma(n+2)+math.lgamma(2*ν+3)-math.lgamma(ν+2) 
                 +math.lgamma(n+ν+m-μ+1)-math.lgamma(n-m+1)-math.lgamma(ν+μ+1)
@@ -450,6 +483,7 @@ def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
     return presum * np.sum(summandq)
     
 # ZDE OPĚT JINAK ZNAMÉNKA než v Xu (J. comp. phys 127, 285)
+#@jit
 def B̃(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
     exponent=(math.lgamma(2*n+3)-math.lgamma(n+2)+math.lgamma(2*ν+3)-math.lgamma(ν+2) 
                 +math.lgamma(n+ν+m-μ+2)-math.lgamma(n-m+1)-math.lgamma(ν+μ+1)
@@ -579,6 +613,7 @@ def mie_coefficients(a, nmax,  #ω, ε_i, ε_e=1, J_ext=1, J_scat=3
     TH = -(( η_inv_e * že * zs - η_inv_e * ze * žs)/(-η_inv_i * ži * zs + η_inv_e * zi * žs)) 
     return (RH, RV, TH, TV)
 
+#@jit
 def G_Mie_scat_precalc_cart_new(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3):
     """
     Implementation according to Kristensson, page 50

tests/test_qpms_p.py

@@ -1,16 +1,70 @@
+"""
+Unit tests for qpms_p
+=====================
+
+Covered functions
+-----------------
+plane_pq_y vs. vswf_yr1
+
+Not covered
+-----------
+Everything else
+
+"""
+
 import unittest
 import qpms
 import numpy as np
+from numpy import newaxis as ň
+import warnings
 
 # Some constants go here.
-maxx = 20 # The maximum argument of the Bessel's functions, i.e. maximum wave number times the distance
-lMax = 25 # To which order we decompose the waves
 
+# 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):
-        pass
+        for oom in lengthOrdersOfMagnitude:
+            k = np.random.randn(nsamples, 3) / oom
+            ksiz = np.linalg.norm(k, axis=-1)
+            kdir = k / ksiz[...,ň]
+            E_0 = np.cross(np.random.randn(nsamples, 3), k) * oom # ensure orthogonality
+            for s in range(nsamples):
+                testpoints = oom * maxx * np.random.randn(npoints, 3)
+                p, q = qpms.plane_pq_y(lMax, k[s], E_0[s])
+                planewave_1 =  np.exp(1j*np.dot(testpoints,k[s]))[:,ň] * E_0[s,:]
+                for i in range(npoints):
+                    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
 
     def testCornerCases(self):
         pass