QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

dudom 

Former-commit-id: 77250aa7b2ea2931e3cf1e3f401919cd19ba829d
This commit is contained in:
Marek Nečada 2016-12-15 16:15:01 +02:00
parent 82259b77f5
commit c18e1209f4
4 changed files with 120 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
qpms/__init__.py
View File

@ -3,3 +3,4 @@ __version__ = get_distribution('qpms').version
from qpms_c import * from qpms_c import *
from .qpms_p import * from .qpms_p import *
from .lattices import *

74
qpms/lattices.py
View File

@ -5,8 +5,8 @@ import numpy as np
import time import time
import scipy import scipy
import sys import sys
from qpms_c import * # TODO be explicit about what is imported from qpms_c import get_mn_y # TODO be explicit about what is imported
from .qpms_p import nelem2lMax # TODO be explicit about what is imported from .qpms_p import cart2sph, nelem2lMax, Ã, B̃ # TODO be explicit about what is imported
class Scattering(object): class Scattering(object):
''' '''
@ -39,36 +39,82 @@ class Scattering(object):
k_0 (float): Wave number for the space between scatterers. k_0 (float): Wave number for the space between scatterers.
lMax (int): Absolute maximum l for all scatterers. Depending on implementation, lMax (int): Absolute maximum l for all scatterers. Depending on implementation,
lMax can be smaller for some individual scatterers in certain subclasses. lMax can be smaller for some individual scatterers in certain subclasses.
FIXME: here it is still implemented as constant lMax for all sites, see #!
prepared (bool): Keeps information whether the interaction matrix has prepared (bool): Keeps information whether the interaction matrix has
already been built and factorized. already been built and factorized.
''' '''
def __init__(self, positions, TMatrices, k_0, lMax = None, verbose=False):
def __init__(self, positions, TMatrices, k_0, lMax = None, verbose=False, J_scat=3):
self.J_scat = J_scat
self.positions = positions self.positions = positions
self.TMatrices = TMatrices self.interaction_matrix = None
self.N = positions.shape[0]
self.k_0 = k_0 self.k_0 = k_0
self.lMax = lMax ? lMax : nelem2lMax(TMatrices.shape[-1]) self.lMax = lMax if lMax else nelem2lMax(TMatrices.shape[-1])
nelem = lMax * (lMax + 2) #!
self.nelem = nelem #!
self.prepared = False self.prepared = False
self.TMatrices = np.broadcast_to(TMatrices, (self.N,2,nelem,2,nelem))
def prepare(self, keep_interaction_matrix = False, verbose=False): def prepare(self, keep_interaction_matrix = False, verbose=False):
if not self.prepared: if not self.prepared:
if not self.interaction_matrix: if not self.interaction_matrix:
self.build_interaction_matrix(verbose=verbose) self.build_interaction_matrix(verbose=verbose)
self.lupiv = scipy.linalg_lu_factor(interaction_matrix) self.lupiv = scipy.linalg.lu_factor(self.interaction_matrix,overwrite_a = not keep_interaction_matrix)
if not keep_interaction_matrix: if not keep_interaction_matrix:
self.interaction_matrix = None self.interaction_matrix = None
self.prepared = True self.prepared = True
def build_interaction_matrix(verbose = False): def build_interaction_matrix(self,verbose = False):
pass N = self.N
my, ny = get_mn_y(self.lMax)
nelem = len(my)
leftmatrix = np.zeros((N,2,nelem,N,2,nelem), dtype=complex)
for i in range(N):
for j in range(N):
for yi in range(nelem):
for yj in range(nelem):
if(i != j):
d_i2j = cart2sph(self.positions[j]-self.positions[i])
a = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat)
b = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat)
leftmatrix[j,0,yj,i,0,yi] = a
leftmatrix[j,1,yj,i,1,yi] = a
leftmatrix[j,0,yj,i,1,yi] = b
leftmatrix[j,1,yj,i,0,yi] = b
# at this point, leftmatrix is the translation matrix
n2id = np.identity(2*nelem)
n2id.shape = (2,nelem,2,nelem)
for j in range(N):
leftmatrix[j] = - np.tensordot(self.TMatrices[j],leftmatrix[j],axes=([-2,-1],[0,1]))
# at this point, jth row of leftmatrix is that of -MT
leftmatrix[j,:,:,j,:,:] += n2id
# now we are done, 1-MT
leftmatrix.shape=(N*2*nelem,N*2*nelem)
self.interaction_matrix = leftmatrix
def scatter(pq_0_c, verbose = False): def scatter_constmultipole(self, pq_0_c, verbose = False):
N = self.N
self.prepare(verbose=verbose) self.prepare(verbose=verbose)
nelem = self.nelem
if(pq_0_c ==1):
pq_0_c = np.full((2,nelem),1)
ab = np.empty((2,nelem,N*2*nelem), dtype=complex)
for N_or_M in range(2):
for yy in range(nelem):
pq_0 = np.zeros((2,nelem),dtype=np.complex_)
pq_0[N_or_M,yy] = pq_0_c[N_or_M,yy]
pq_0 = np.broadcast_to(pq_0, (N,2,nelem))
MP_0 = np.empty((N,2,nelem),dtype=np.complex_)
for j in range(N):
MP_0[j] = np.tensordot(self.TMatrices[j], pq_0[j],axes=([-2,-1],[-2,-1]))
MP_0.shape = (N*2*nelem,)
a[N_or_M,yy] = scipy.linalg.lu_solve(lupiv,MP_0)
ab.shape = (2,nelem,N,2,nelem)
return ab
class Scattering_lattice(Scattering):
def __init__(self):
pass pass

116
qpms/qpms_p.py
View File

@ -29,7 +29,7 @@ except ImportError:
Accordingly, we define our own jit decorator that handles Accordingly, we define our own jit decorator that handles
different versions of numba or does nothing if numba is not different versions of numba or does nothing if numba is not
present. Note that functions that include unicode identifiers present. Note that functions that include unicode identifiers
must be decorated with @ujit must be decorated with #@ujit
''' '''
#def dummywrap(f): #def dummywrap(f):
# return f # return f
@ -46,7 +46,7 @@ def ujit(f):
# Coordinate transforms for arrays of "arbitrary" shape # Coordinate transforms for arrays of "arbitrary" shape
@ujit #@ujit
def cart2sph(cart,axis=-1): def cart2sph(cart,axis=-1):
if (cart.shape[axis] != 3): if (cart.shape[axis] != 3):
raise ValueError("The converted array has to have dimension 3" raise ValueError("The converted array has to have dimension 3"
@ -58,7 +58,7 @@ def cart2sph(cart,axis=-1):
φ = np.arctan2(y,x) # arctan2 handles zeroes correctly itself φ = np.arctan2(y,x) # arctan2 handles zeroes correctly itself
return np.concatenate((r,θ,φ),axis=axis) return np.concatenate((r,θ,φ),axis=axis)
@ujit #@ujit
def sph2cart(sph, axis=-1): def sph2cart(sph, axis=-1):
if (sph.shape[axis] != 3): if (sph.shape[axis] != 3):
raise ValueError("The converted array has to have dimension 3" raise ValueError("The converted array has to have dimension 3"
@ -70,7 +70,7 @@ def sph2cart(sph, axis=-1):
z = r * np.cos(θ) z = r * np.cos(θ)
return np.concatenate((x,y,z),axis=axis) return np.concatenate((x,y,z),axis=axis)
@ujit #@ujit
def sph_loccart2cart(loccart, sph, axis=-1): def sph_loccart2cart(loccart, sph, axis=-1):
""" """
Transformation of vector specified in local orthogonal coordinates Transformation of vector specified in local orthogonal coordinates
@ -120,7 +120,7 @@ def sph_loccart2cart(loccart, sph, axis=-1):
out=inr̂*r̂+inθ̂*θ̂+inφ̂*φ̂ out=inr̂*r̂+inθ̂*θ̂+inφ̂*φ̂
return out return out
@ujit #@ujit
def sph_loccart_basis(sph, sphaxis=-1, cartaxis=None): def sph_loccart_basis(sph, sphaxis=-1, cartaxis=None):
""" """
Returns the local cartesian basis in terms of global cartesian basis. Returns the local cartesian basis in terms of global cartesian basis.
@ -156,7 +156,7 @@ def sph_loccart_basis(sph, sphaxis=-1, cartaxis=None):
out = np.concatenate((x,y,z),axis=cartaxis) out = np.concatenate((x,y,z),axis=cartaxis)
return out return out
@jit #@jit
def lpy(nmax, z): def lpy(nmax, z):
""" """
Associated legendre function and its derivatative at z in the 'y-indexing'. Associated legendre function and its derivatative at z in the 'y-indexing'.
@ -233,18 +233,18 @@ def vswf_yr(pos_sph,nmax,J=1):
pass pass
from scipy.special import sph_jn, sph_yn from scipy.special import sph_jn, sph_yn
@jit #@jit
def _sph_zn_1(n,z): def _sph_zn_1(n,z):
return sph_jn(n,z) return sph_jn(n,z)
@jit #@jit
def _sph_zn_2(n,z): def _sph_zn_2(n,z):
return sph_yn(n,z) return sph_yn(n,z)
@jit #@jit
def _sph_zn_3(n,z): def _sph_zn_3(n,z):
besj=sph_jn(n,z) besj=sph_jn(n,z)
besy=sph_yn(n,z) besy=sph_yn(n,z)
return (besj[0] + 1j*besy[0],besj[1] + 1j*besy[1]) return (besj[0] + 1j*besy[0],besj[1] + 1j*besy[1])
@jit #@jit
def _sph_zn_4(n,z): def _sph_zn_4(n,z):
besj=sph_jn(n,z) besj=sph_jn(n,z)
besy=sph_yn(n,z) besy=sph_yn(n,z)
@ -253,7 +253,7 @@ _sph_zn = [_sph_zn_1,_sph_zn_2,_sph_zn_3,_sph_zn_4]
# computes bessel/hankel functions for orders from 0 up to n; drops # computes bessel/hankel functions for orders from 0 up to n; drops
# the derivatives which are also included in scipy.special.sph_jn/yn # the derivatives which are also included in scipy.special.sph_jn/yn
@jit #@jit
def zJn(n, z, J=1): def zJn(n, z, J=1):
return _sph_zn[J-1](n=n,z=z) return _sph_zn[J-1](n=n,z=z)
@ -263,7 +263,7 @@ def zJn(n, z, J=1):
# FIXME: this can be expressed simply as: # FIXME: this can be expressed simply as:
# $$ -\frac{1}{2}\sqrt{\frac{2n+1}{4\pi}n\left(n+1\right)}(\delta_{m,1}+\delta_{m,-1}) $$ # $$ -\frac{1}{2}\sqrt{\frac{2n+1}{4\pi}n\left(n+1\right)}(\delta_{m,1}+\delta_{m,-1}) $$
@ujit #@ujit
def π̃_zerolim(nmax): # seems OK def π̃_zerolim(nmax): # seems OK
""" """
lim_{θ→ 0-} π̃(cos θ) lim_{θ→ 0-} π̃(cos θ)
@ -281,7 +281,7 @@ def π̃_zerolim(nmax): # seems OK
π̃_y = prenorm * π̃_y π̃_y = prenorm * π̃_y
return π̃_y return π̃_y
@ujit #@ujit
def π̃_pilim(nmax): # Taky OK, jen to možná není kompatibilní se vzorečky z mathematiky def π̃_pilim(nmax): # Taky OK, jen to možná není kompatibilní se vzorečky z mathematiky
""" """
lim_{θ→ π+} π̃(cos θ) lim_{θ→ π+} π̃(cos θ)
@ -301,7 +301,7 @@ def π̃_pilim(nmax): # Taky OK, jen to možná není kompatibilní se vzorečky
# FIXME: this can be expressed simply as # FIXME: this can be expressed simply as
# $$ -\frac{1}{2}\sqrt{\frac{2n+1}{4\pi}n\left(n+1\right)}(\delta_{m,1}-\delta_{m,-1}) $$ # $$ -\frac{1}{2}\sqrt{\frac{2n+1}{4\pi}n\left(n+1\right)}(\delta_{m,1}-\delta_{m,-1}) $$
@ujit #@ujit
def τ̃_zerolim(nmax): def τ̃_zerolim(nmax):
""" """
lim_{θ→ 0-} τ̃(cos θ) lim_{θ→ 0-} τ̃(cos θ)
@ -312,7 +312,7 @@ def τ̃_zerolim(nmax):
p0[minus1mmask] = -p0[minus1mmask] p0[minus1mmask] = -p0[minus1mmask]
return p0 return p0
@ujit #@ujit
def τ̃_pilim(nmax): def τ̃_pilim(nmax):
""" """
lim_{θ→ π+} τ̃(cos θ) lim_{θ→ π+} τ̃(cos θ)
@ -323,7 +323,7 @@ def τ̃_pilim(nmax):
t[plus1mmask] = -t[plus1mmask] t[plus1mmask] = -t[plus1mmask]
return t return t
@ujit #@ujit
def get_π̃τ̃_y1(θ,nmax): def get_π̃τ̃_y1(θ,nmax):
# TODO replace with the limit functions (below) when θ approaches # TODO replace with the limit functions (below) when θ approaches
# the extreme values at about 1e-6 distance # the extreme values at about 1e-6 distance
@ -343,7 +343,7 @@ def get_π̃τ̃_y1(θ,nmax):
τ̃_y = prenorm * dPy * (- math.sin(θ)) # TADY BACHA!!!!!!!!!! * (- math.sin(pos_sph[1])) ??? τ̃_y = prenorm * dPy * (- math.sin(θ)) # TADY BACHA!!!!!!!!!! * (- math.sin(pos_sph[1])) ???
return (π̃_y,τ̃_y) return (π̃_y,τ̃_y)
@ujit #@ujit
def vswf_yr1(pos_sph,nmax,J=1): def vswf_yr1(pos_sph,nmax,J=1):
""" """
As vswf_yr, but evaluated only at single position (i.e. pos_sph has As vswf_yr, but evaluated only at single position (i.e. pos_sph has
@ -400,7 +400,7 @@ def vswf_yr1(pos_sph,nmax,J=1):
# return 1j**ny * np.sqrt((2*ny+1)*factorial(ny-my) / # return 1j**ny * np.sqrt((2*ny+1)*factorial(ny-my) /
# (ny*(ny+1)*factorial(ny+my)) # (ny*(ny+1)*factorial(ny+my))
# ) # )
@ujit #@ujit
def zplane_pq_y(nmax, betap = 0): def zplane_pq_y(nmax, betap = 0):
""" """
The z-propagating plane wave expansion coefficients as in [1, (1.12)] The z-propagating plane wave expansion coefficients as in [1, (1.12)]
@ -419,7 +419,7 @@ def zplane_pq_y(nmax, betap = 0):
#import warnings #import warnings
@ujit #@ujit
def plane_pq_y(nmax, kdir_cart, E_cart): def plane_pq_y(nmax, kdir_cart, E_cart):
""" """
The plane wave expansion coefficients for any direction kdir_cart The plane wave expansion coefficients for any direction kdir_cart
@ -476,13 +476,13 @@ def plane_pq_y(nmax, kdir_cart, E_cart):
# Functions copied from scattering_xu, additionaly normalized # Functions copied from scattering_xu, additionaly normalized
from py_gmm.gmm import vec_trans as vc from py_gmm.gmm import vec_trans as vc
@ujit #@ujit
def q_max(m,n,μ,ν): def q_max(m,n,μ,ν):
return min(n,ν,(n+ν-abs(m+μ))/2) return min(n,ν,(n+ν-abs(m+μ))/2)
# returns array with indices corresponding to q # returns array with indices corresponding to q
# argument q does nothing for now # argument q does nothing for now
@ujit #@ujit
def a_q(m,n,μ,ν,q = None): def a_q(m,n,μ,ν,q = None):
qm=q_max(m,n,μ,ν) qm=q_max(m,n,μ,ν)
res, err= vc.gaunt_xu(m,n,μ,ν,qm) res, err= vc.gaunt_xu(m,n,μ,ν,qm)
@ -493,7 +493,7 @@ def a_q(m,n,μ,ν,q = None):
# All arguments are single numbers (for now) # All arguments are single numbers (for now)
# ZDE VYCHÁZEJÍ DIVNÁ ZNAMÉNKA # ZDE VYCHÁZEJÍ DIVNÁ ZNAMÉNKA
@ujit #@ujit
def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J): def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
""" """
The à translation coefficient for spherical vector waves. The à translation coefficient for spherical vector waves.
@ -552,7 +552,7 @@ def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
return presum * np.sum(summandq) return presum * np.sum(summandq)
# ZDE OPĚT JINAK ZNAMÉNKA než v Xu (J. comp. phys 127, 285) # ZDE OPĚT JINAK ZNAMÉNKA než v Xu (J. comp. phys 127, 285)
@ujit #@ujit
def B̃(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J): def B̃(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
""" """
The B̃ translation coefficient for spherical vector waves. The B̃ translation coefficient for spherical vector waves.
@ -617,7 +617,7 @@ def B̃(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J):
# In[7]: # In[7]:
# Material parameters # Material parameters
@ujit #@ujit
def ε_drude(ε_inf, ω_p, γ_p, ω): # RELATIVE permittivity, of course def ε_drude(ε_inf, ω_p, γ_p, ω): # RELATIVE permittivity, of course
return ε_inf - ω_p*ω_p/(ω*(ω+1j*γ_p)) return ε_inf - ω_p*ω_p/(ω*(ω+1j*γ_p))
@ -625,7 +625,7 @@ def ε_drude(ε_inf, ω_p, γ_p, ω): # RELATIVE permittivity, of course
# In[8]: # In[8]:
# Mie scattering # Mie scattering
@ujit #@ujit
def mie_coefficients(a, nmax, #ω, ε_i, ε_e=1, J_ext=1, J_scat=3 def mie_coefficients(a, nmax, #ω, ε_i, ε_e=1, J_ext=1, J_scat=3
k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3): k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3):
""" """
@ -705,7 +705,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)) TH = -(( η_inv_e * že * zs - η_inv_e * ze * žs)/(-η_inv_i * ži * zs + η_inv_e * zi * žs))
return (RH, RV, TH, TV) return (RH, RV, TH, TV)
@ujit #@ujit
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): 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 Implementation according to Kristensson, page 50
@ -742,7 +742,7 @@ def G_Mie_scat_precalc_cart_new(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_
RV[ny][:,ň,ň] * Ñlo_cart_y[:,:,ň].conj() * Ñhi_cart_y[:,ň,:]) / (ny * (ny+1))[:,ň,ň] RV[ny][:,ň,ň] * Ñlo_cart_y[:,:,ň].conj() * Ñhi_cart_y[:,ň,:]) / (ny * (ny+1))[:,ň,ň]
return 1j* k_e*np.sum(G_y,axis=0) return 1j* k_e*np.sum(G_y,axis=0)
@ujit #@ujit
def G_Mie_scat_precalc_cart(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3): def G_Mie_scat_precalc_cart(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3):
""" """
r1_cart (destination), r2_cart (source) and the result are in cartesian coordinates r1_cart (destination), r2_cart (source) and the result are in cartesian coordinates
@ -797,7 +797,7 @@ def G_Mie_scat_precalc_cart(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e,
G_source_dest = sph_loccart2cart(G_source_dest, sph=orig2dest_sph, axis=-1) G_source_dest = sph_loccart2cart(G_source_dest, sph=orig2dest_sph, axis=-1)
return G_source_dest return G_source_dest
@ujit #@ujit
def G_Mie_scat_cart(source_cart, dest_cart, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3): def G_Mie_scat_cart(source_cart, dest_cart, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3):
""" """
TODO TODO
@ -817,7 +817,7 @@ def cross_section_Mie(a, nmax, k_i, k_e, μ_i, μ_e,):
# In[9]: # In[9]:
# From PRL 112, 253601 (1) # From PRL 112, 253601 (1)
@ujit #@ujit
def Grr_Delga(nmax, a, r, k, ε_m, ε_b): def Grr_Delga(nmax, a, r, k, ε_m, ε_b):
om = k * c om = k * c
z = (r-a)/a z = (r-a)/a
@ -839,7 +839,7 @@ def Grr_Delga(nmax, a, r, k, ε_m, ε_b):
# Test if the decomposition of plane wave works also for absorbing environment (complex k). # Test if the decomposition of plane wave works also for absorbing environment (complex k).
# From PRL 112, 253601 (1) # From PRL 112, 253601 (1)
@ujit #@ujit
def Grr_Delga(nmax, a, r, k, ε_m, ε_b): def Grr_Delga(nmax, a, r, k, ε_m, ε_b):
om = k * c om = k * c
z = (r-a)/a z = (r-a)/a
@ -848,7 +848,7 @@ def Grr_Delga(nmax, a, r, k, ε_m, ε_b):
s = np.sum( (n+1)**2 * (ε_m-ε_b) / ((1+z)**(2*n+4) * (ε_m + ((n+1)/n)*ε_b))) s = np.sum( (n+1)**2 * (ε_m-ε_b) / ((1+z)**(2*n+4) * (ε_m + ((n+1)/n)*ε_b)))
return (g0 + s * c**2/(4*π*om**2*ε_b*a**3)) return (g0 + s * c**2/(4*π*om**2*ε_b*a**3))
@ujit #@ujit
def G0_dip_1(r_cart,k): def G0_dip_1(r_cart,k):
""" """
Free-space dyadic Green's function in terms of the spherical vector waves. Free-space dyadic Green's function in terms of the spherical vector waves.
@ -865,15 +865,15 @@ def G0_dip_1(r_cart,k):
# Free-space dyadic Green's functions from RMP 70, 2, 447 =: [1] # Free-space dyadic Green's functions from RMP 70, 2, 447 =: [1]
# (The numerical value is correct only at the regular part, i.e. r != 0) # (The numerical value is correct only at the regular part, i.e. r != 0)
@ujit #@ujit
def _P(z): def _P(z):
return (1-1/z+1/(z*z)) return (1-1/z+1/(z*z))
@ujit #@ujit
def _Q(z): def _Q(z):
return (-1+3/z-3/(z*z)) return (-1+3/z-3/(z*z))
# [1, (9)] FIXME The sign here is most likely wrong!!! # [1, (9)] FIXME The sign here is most likely wrong!!!
@ujit #@ujit
def G0_analytical(r #cartesian! def G0_analytical(r #cartesian!
, k): , k):
I=np.identity(3) I=np.identity(3)
@ -887,7 +887,7 @@ def G0_analytical(r #cartesian!
)) ))
# [1, (11)] # [1, (11)]
@ujit #@ujit
def G0L_analytical(r, k): def G0L_analytical(r, k):
I=np.identity(3) I=np.identity(3)
rn = sph_loccart2cart(np.array([1.,0.,0.]), cart2sph(r), axis=-1) rn = sph_loccart2cart(np.array([1.,0.,0.]), cart2sph(r), axis=-1)
@ -896,11 +896,11 @@ def G0L_analytical(r, k):
return (I-3*rnxrn)/(4*π*k*k*r**3)[...,ň,ň] return (I-3*rnxrn)/(4*π*k*k*r**3)[...,ň,ň]
# [1,(10)] # [1,(10)]
@jit #@jit
def G0T_analytical(r, k): def G0T_analytical(r, k):
return G0_analytical(r,k) - G0L_analytical(r,k) return G0_analytical(r,k) - G0L_analytical(r,k)
@ujit #@ujit
def G0_sum_1_slow(source_cart, dest_cart, k, nmax): def G0_sum_1_slow(source_cart, dest_cart, k, nmax):
my, ny = get_mn_y(nmax) my, ny = get_mn_y(nmax)
nelem = len(my) nelem = len(my)
@ -911,7 +911,7 @@ def G0_sum_1_slow(source_cart, dest_cart, k, nmax):
# Transformations of spherical bases # Transformations of spherical bases
@jit #@jit
def WignerD_mm(l, quat): def WignerD_mm(l, quat):
""" """
Calculates Wigner D matrix (as an numpy (2*l+1,2*l+1)-shaped array) Calculates Wigner D matrix (as an numpy (2*l+1,2*l+1)-shaped array)
@ -925,7 +925,7 @@ def WignerD_mm(l, quat):
Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1) Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1)
return Delems return Delems
@jit #@jit
def WignerD_mm_fromvector(l, vect): def WignerD_mm_fromvector(l, vect):
""" """
TODO doc TODO doc
@ -933,7 +933,7 @@ def WignerD_mm_fromvector(l, vect):
return WignerD_mm(l, quaternion.from_rotation_vector(vect)) return WignerD_mm(l, quaternion.from_rotation_vector(vect))
@jit #@jit
def WignerD_yy(lmax, quat): def WignerD_yy(lmax, quat):
""" """
TODO doc TODO doc
@ -948,7 +948,7 @@ def WignerD_yy(lmax, quat):
b_in = e_in b_in = e_in
return Delems return Delems
@jit #@jit
def WignerD_yy_fromvector(lmax, vect): def WignerD_yy_fromvector(lmax, vect):
""" """
TODO doc TODO doc
@ -956,7 +956,7 @@ def WignerD_yy_fromvector(lmax, vect):
return WignerD_yy(lmax, quaternion.from_rotation_vector(vect)) return WignerD_yy(lmax, quaternion.from_rotation_vector(vect))
@jit #@jit
def xflip_yy(lmax): def xflip_yy(lmax):
""" """
TODO doc TODO doc
@ -973,12 +973,12 @@ def xflip_yy(lmax):
b_in = e_in b_in = e_in
return elems return elems
@jit #@jit
def xflip_tyy(lmax): def xflip_tyy(lmax):
fl_yy = xflip_yy(lmax) fl_yy = xflip_yy(lmax)
return np.array([fl_yy,-fl_yy]) return np.array([fl_yy,-fl_yy])
@jit #@jit
def xflip_tyty(lmax): def xflip_tyty(lmax):
fl_yy = xflip_yy(lmax) fl_yy = xflip_yy(lmax)
nelem = fl_yy.shape[0] nelem = fl_yy.shape[0]
@ -987,7 +987,7 @@ def xflip_tyty(lmax):
fl_tyty[1,:,1,:] = -fl_yy fl_tyty[1,:,1,:] = -fl_yy
return fl_tyty return fl_tyty
@jit #@jit
def yflip_yy(lmax): def yflip_yy(lmax):
""" """
TODO doc TODO doc
@ -999,12 +999,12 @@ def yflip_yy(lmax):
elems[(my % 2)==1] = elems[(my % 2)==1] * -1 # Obvious sign of tiredness (this is correct but ugly; FIXME) elems[(my % 2)==1] = elems[(my % 2)==1] * -1 # Obvious sign of tiredness (this is correct but ugly; FIXME)
return elems return elems
@jit #@jit
def yflip_tyy(lmax): def yflip_tyy(lmax):
fl_yy = yflip_yy(lmax) fl_yy = yflip_yy(lmax)
return np.array([fl_yy,-fl_yy]) return np.array([fl_yy,-fl_yy])
@jit #@jit
def yflip_tyty(lmax): def yflip_tyty(lmax):
fl_yy = yflip_yy(lmax) fl_yy = yflip_yy(lmax)
nelem = fl_yy.shape[0] nelem = fl_yy.shape[0]
@ -1013,7 +1013,7 @@ def yflip_tyty(lmax):
fl_tyty[1,:,1,:] = -fl_yy fl_tyty[1,:,1,:] = -fl_yy
return fl_tyty return fl_tyty
@jit #@jit
def zflip_yy(lmax): def zflip_yy(lmax):
""" """
TODO doc TODO doc
@ -1029,12 +1029,12 @@ def zflip_yy(lmax):
b_in = e_in b_in = e_in
return elems return elems
@jit #@jit
def zflip_tyy(lmax): def zflip_tyy(lmax):
fl_yy = zflip_yy(lmax) fl_yy = zflip_yy(lmax)
return np.array([fl_yy,-fl_yy]) return np.array([fl_yy,-fl_yy])
@jit #@jit
def zflip_tyty(lmax): def zflip_tyty(lmax):
fl_yy = zflip_yy(lmax) fl_yy = zflip_yy(lmax)
nelem = fl_yy.shape[0] nelem = fl_yy.shape[0]
@ -1043,7 +1043,7 @@ def zflip_tyty(lmax):
fl_tyty[1,:,1,:] = -fl_yy fl_tyty[1,:,1,:] = -fl_yy
return fl_tyty return fl_tyty
@jit #@jit
def parity_yy(lmax): def parity_yy(lmax):
""" """
Parity operator (flip in x,y,z) Parity operator (flip in x,y,z)
@ -1061,7 +1061,7 @@ def parity_yy(lmax):
#----------------------------------------------------# #----------------------------------------------------#
# We don't really need this particular function anymore, but... # We don't really need this particular function anymore, but...
@jit #@jit
def _scuffTMatrixConvert_EM_01(EM): def _scuffTMatrixConvert_EM_01(EM):
#print(EM) #print(EM)
if (EM == b'E'): if (EM == b'E'):
@ -1071,7 +1071,7 @@ def _scuffTMatrixConvert_EM_01(EM):
else: else:
return None return None
@ujit #@ujit
def loadScuffTMatrices(fileName): def loadScuffTMatrices(fileName):
""" """
TODO doc TODO doc
@ -1107,7 +1107,7 @@ def loadScuffTMatrices(fileName):
# misc tensor maniputalion # misc tensor maniputalion
@jit #@jit
def apply_matrix_left(matrix, tensor, axis): def apply_matrix_left(matrix, tensor, axis):
""" """
TODO doc TODO doc
@ -1118,7 +1118,7 @@ def apply_matrix_left(matrix, tensor, axis):
tmp = np.tensordot(matrix, tensor, axes=(-1,axis)) tmp = np.tensordot(matrix, tensor, axes=(-1,axis))
return np.moveaxis(tmp, 0, axis) return np.moveaxis(tmp, 0, axis)
@jit #@jit
def apply_ndmatrix_left(matrix,tensor,axes): def apply_ndmatrix_left(matrix,tensor,axes):
""" """
Generalized apply_matrix_left, the matrix can have more (2N) abstract dimensions, Generalized apply_matrix_left, the matrix can have more (2N) abstract dimensions,
@ -1134,7 +1134,7 @@ def apply_ndmatrix_left(matrix,tensor,axes):
# Array simulations # Array simulations
#################### ####################
@jit #@jit
def nelem2lMax(nelem): def nelem2lMax(nelem):
""" """
Auxiliary inverse function to nelem(lMax) = (lMax + 2) * lMax. Returns 0 if Auxiliary inverse function to nelem(lMax) = (lMax + 2) * lMax. Returns 0 if
@ -1190,7 +1190,7 @@ def scatter_plane_wave(omega, epsilon_b, positions, Tmatrices, k_dirs, E_0s, #sa
pass pass
import warnings import warnings
@ujit #@ujit
def scatter_plane_wave_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, k_dirs, E_0s, def scatter_plane_wave_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, k_dirs, E_0s,
return_pq_0 = False, return_pq= False, return_xy = False, watch_time = False): return_pq_0 = False, return_pq= False, return_xy = False, watch_time = False):
""" """
@ -1422,7 +1422,7 @@ def scatter_plane_wave_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, k_
import warnings import warnings
@ujit #@ujit
def scatter_constmultipole_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, pq_0_c = 1, def scatter_constmultipole_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, pq_0_c = 1,
return_pq= False, return_xy = False, watch_time = False): return_pq= False, return_xy = False, watch_time = False):
""" """

2
setup.py
View File

@ -12,7 +12,7 @@ qpms_c = Extension('qpms_c',
sources = ['qpms/qpms_c.pyx']) sources = ['qpms/qpms_c.pyx'])
setup(name='qpms', setup(name='qpms',
version = "0.1.5", version = "0.1.6",
packages=['qpms'], packages=['qpms'],
# setup_requires=['setuptools_cython'], # setup_requires=['setuptools_cython'],
install_requires=['cython>=0.21','quaternion','spherical_functions','py_gmm'], install_requires=['cython>=0.21','quaternion','spherical_functions','py_gmm'],