From 30013ac7fa1d66d54132df121492055809f80095 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 26 Jul 2016 00:42:16 +0300 Subject: [PATCH] Plane wave test, docs. Former-commit-id: 86bdabea5c1a424401d3b93ea98af16bd2e963be --- README.rst | 12 +++++++++ qpms/qpms_p.py | 35 +++++++++++++++++++++++++ tests/test_qpms_p.py | 62 +++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 105 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index cee538a..5b74a25 100644 --- a/README.rst +++ b/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 diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py index 1c64a41..8c677cc 100644 --- a/qpms/qpms_p.py +++ b/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 diff --git a/tests/test_qpms_p.py b/tests/test_qpms_p.py index 2b6b354..63a422d 100644 --- a/tests/test_qpms_p.py +++ b/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