1531 lines
63 KiB
Cython
1531 lines
63 KiB
Cython
"""@package qpms_c
|
|
Cythonized parts of QPMS; mostly wrappers over the C data structures
|
|
to make them available in Python.
|
|
"""
|
|
|
|
# Cythonized parts of QPMS here
|
|
# -----------------------------
|
|
|
|
import numpy as np
|
|
from .qpms_cdefs cimport *
|
|
from .cyquaternions cimport IRot3, CQuat
|
|
from .cybspec cimport BaseSpec
|
|
from .cycommon cimport make_c_string
|
|
from .cycommon import string_c2py, PointGroupClass, BesselType
|
|
from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
|
|
from .cymaterials cimport EpsMuGenerator, EpsMu
|
|
from libc.stdlib cimport malloc, free, calloc
|
|
import warnings
|
|
|
|
def qpms_library_version():
|
|
return "git_" + qpms_version_git_sha1.decode('utf-8')
|
|
|
|
from cython.parallel cimport prange, parallel
|
|
from cython cimport boundscheck, wraparound
|
|
|
|
# Set custom GSL error handler. N.B. this is obviously not thread-safe.
|
|
cdef char *pgsl_err_reason
|
|
cdef char *pgsl_err_file
|
|
cdef int pgsl_err_line
|
|
cdef int pgsl_errno = 0
|
|
cdef int *pgsl_errno_ignorelist = NULL # list of ignored error codes, terminated by zero
|
|
|
|
# This error handler only sets the variables above
|
|
cdef void pgsl_error_handler(const char *reason, const char *_file, const int line, const int gsl_errno):
|
|
global pgsl_err_reason, pgsl_err_file, pgsl_err_line, pgsl_errno, pgsl_errno_ignorelist
|
|
cdef size_t i
|
|
if(pgsl_errno_ignorelist):
|
|
i = 0
|
|
while pgsl_errno_ignorelist[i] != 0:
|
|
if gsl_errno == pgsl_errno_ignorelist[i]:
|
|
return
|
|
i += 1
|
|
pgsl_err_file = _file
|
|
pgsl_err_reason = reason
|
|
pgsl_errno = gsl_errno
|
|
pgsl_err_line = line
|
|
return
|
|
|
|
cdef const int* pgsl_set_ignorelist(const int *new_ignorelist):
|
|
global pgsl_errno_ignorelist
|
|
cdef const int *oldlist = pgsl_errno_ignorelist
|
|
pgsl_errno_ignorelist = new_ignorelist
|
|
return oldlist
|
|
|
|
cdef class pgsl_ignore_error():
|
|
'''Context manager for setting a temporary list of errnos ignored by pgsl_error_handler.
|
|
Always sets pgsl_error_handler.
|
|
|
|
Performs pgsl_check_err() on exit unless
|
|
'''
|
|
cdef const int *ignorelist_old
|
|
cdef gsl_error_handler_t *old_handler
|
|
cdef bint final_check
|
|
cdef object ignorelist_python
|
|
|
|
cdef int *ignorelist
|
|
def __cinit__(self, *ignorelist, **kwargs):
|
|
self.ignorelist = <int*>calloc((len(ignorelist)+1), sizeof(int))
|
|
self.ignorelist_python = ignorelist
|
|
for i in range(len(ignorelist)):
|
|
self.ignorelist[i] = ignorelist[i]
|
|
if "final_check" in kwargs.keys() and not kwargs["final_check"]:
|
|
final_check = True
|
|
final_check = False
|
|
|
|
def __enter__(self):
|
|
global pgsl_error_handler
|
|
self.ignorelist_old = pgsl_set_ignorelist(self.ignorelist)
|
|
self.old_handler = gsl_set_error_handler(pgsl_error_handler)
|
|
return
|
|
|
|
def __exit__(self, type, value, traceback):
|
|
global pgsl_errno_ignorelist, pgsl_error_handler
|
|
pgsl_set_ignorelist(self.ignorelist_old)
|
|
gsl_set_error_handler(self.old_handler)
|
|
if self.final_check:
|
|
pgsl_check_err(retval = None, ignore = self.ignorelist_python)
|
|
|
|
def __dealloc__(self):
|
|
free(self.ignorelist)
|
|
|
|
def pgsl_check_err(retval = None, ignorelist = None):
|
|
global pgsl_err_reason, pgsl_err_file, pgsl_err_line, pgsl_errno
|
|
'''Check for possible errors encountered by pgsl_error_handler.
|
|
Takes return value of a function as an optional argument, which is now ignored.
|
|
'''
|
|
cdef int errno_was
|
|
if (pgsl_errno != 0):
|
|
errno_was = pgsl_errno
|
|
pgsl_errno = 0
|
|
raise RuntimeError("Error %d in GSL calculation in %s:%d: %s" % (errno_was,
|
|
string_c2py(pgsl_err_file), pgsl_err_line, string_c2py(pgsl_err_reason)))
|
|
if (retval is not None and retval != 0 and ignorelist is not None and retval not in ignorelist):
|
|
warnings.warn("Got non-zero return value %d" % retval)
|
|
if retval is not None:
|
|
return retval
|
|
else:
|
|
return 0
|
|
|
|
def set_gsl_pythonic_error_handling():
|
|
'''
|
|
Sets pgsl_error_handler as the GSL error handler to avoid crashing.
|
|
'''
|
|
gsl_set_error_handler(pgsl_error_handler)
|
|
|
|
cdef class PointGroup:
|
|
cdef readonly qpms_pointgroup_t G
|
|
|
|
def __init__(self, cls, qpms_gmi_t n = 0, IRot3 orientation = IRot3()):
|
|
cls = PointGroupClass(cls)
|
|
self.G.c = cls
|
|
if n <= 0 and qpms_pg_is_finite_axial(cls):
|
|
raise ValueError("For finite axial groups, n argument must be positive")
|
|
self.G.n = n
|
|
self.G.orientation = orientation.qd
|
|
|
|
def __len__(self):
|
|
return qpms_pg_order(self.G.c, self.G.n);
|
|
|
|
def __le__(PointGroup self, PointGroup other):
|
|
if qpms_pg_is_subgroup(self.G, other.G):
|
|
return True
|
|
else:
|
|
return False
|
|
def __ge__(PointGroup self, PointGroup other):
|
|
if qpms_pg_is_subgroup(other.G, self.G):
|
|
return True
|
|
else:
|
|
return False
|
|
def __lt__(PointGroup self, PointGroup other):
|
|
return qpms_pg_is_subgroup(self.G, other.G) and not qpms_pg_is_subgroup(other.G, self.G)
|
|
def __eq__(PointGroup self, PointGroup other):
|
|
return qpms_pg_is_subgroup(self.G, other.G) and qpms_pg_is_subgroup(other.G, self.G)
|
|
def __gt__(PointGroup self, PointGroup other):
|
|
return not qpms_pg_is_subgroup(self.G, other.G) and qpms_pg_is_subgroup(other.G, self.G)
|
|
|
|
def elems(self):
|
|
els = list()
|
|
cdef qpms_irot3_t *arr
|
|
arr = qpms_pg_elems(NULL, self.G)
|
|
cdef IRot3 q
|
|
for i in range(len(self)):
|
|
q = IRot3()
|
|
q.cset(arr[i])
|
|
els.append(q)
|
|
free(arr)
|
|
return els
|
|
|
|
cdef class FinitePointGroup:
|
|
'''
|
|
Wrapper over the qpms_finite_group_t structure.
|
|
|
|
TODO more functionality to make it better usable in Python
|
|
(group element class at least)
|
|
'''
|
|
|
|
def __cinit__(self, info):
|
|
'''Constructs a FinitePointGroup from PointGroupInfo'''
|
|
# TODO maybe I might use a try..finally statement to avoid leaks
|
|
# First, generate all basic data from info
|
|
permlist = info.deterministic_elemlist()
|
|
cdef int order = len(permlist)
|
|
permindices = {perm: i for i, perm in enumerate(permlist)} # 'invert' permlist
|
|
identity = info.permgroup.identity
|
|
# We use calloc to avoid calling free to unitialized pointers
|
|
self.G = <qpms_finite_group_t *>calloc(1,sizeof(qpms_finite_group_t))
|
|
if not self.G: raise MemoryError
|
|
self.G[0].name = make_c_string(info.name)
|
|
self.G[0].order = order
|
|
self.G[0].idi = permindices[identity]
|
|
self.G[0].mt = <qpms_gmi_t *>malloc(sizeof(qpms_gmi_t) * order * order)
|
|
if not self.G[0].mt: raise MemoryError
|
|
for i in range(order):
|
|
for j in range(order):
|
|
self.G[0].mt[i*order + j] = permindices[permlist[i] * permlist[j]]
|
|
self.G[0].invi = <qpms_gmi_t *>malloc(sizeof(qpms_gmi_t) * order)
|
|
if not self.G[0].invi: raise MemoryError
|
|
for i in range(order):
|
|
self.G[0].invi[i] = permindices[permlist[i]**-1]
|
|
self.G[0].ngens = len(info.permgroupgens)
|
|
self.G[0].gens = <qpms_gmi_t *>malloc(sizeof(qpms_gmi_t) * self.G[0].ngens)
|
|
if not self.G[0].gens: raise MemoryError
|
|
for i in range(self.G[0].ngens):
|
|
self.G[0].gens[i] = permindices[info.permgroupgens[i]]
|
|
self.G[0].permrep = <char **>calloc(order, sizeof(char *))
|
|
if not self.G[0].permrep: raise MemoryError
|
|
for i in range(order):
|
|
self.G[0].permrep[i] = make_c_string(str(permlist[i]))
|
|
if not self.G[0].permrep[i]: raise MemoryError
|
|
self.G[0].permrep_nelem = info.permgroup.degree
|
|
if info.rep3d is not None:
|
|
self.G[0].rep3d = <qpms_irot3_t *>malloc(order * sizeof(qpms_irot3_t))
|
|
for i in range(order):
|
|
self.G[0].rep3d[i] = info.rep3d[permlist[i]].qd
|
|
self.G[0].nirreps = len(info.irreps)
|
|
self.G[0].irreps = <qpms_finite_group_irrep_t *>calloc(self.G[0].nirreps, sizeof(qpms_finite_group_irrep_t))
|
|
if not self.G[0].irreps: raise MemoryError
|
|
cdef int dim
|
|
for iri, irname in enumerate(sorted(info.irreps.keys())):
|
|
irrep = info.irreps[irname]
|
|
is1d = isinstance(irrep[identity], (int, float, complex))
|
|
dim = 1 if is1d else irrep[identity].shape[0]
|
|
self.G[0].irreps[iri].dim = dim
|
|
self.G[0].irreps[iri].name = <char *>make_c_string(irname)
|
|
if not self.G[0].irreps[iri].name: raise MemoryError
|
|
self.G[0].irreps[iri].m = <cdouble *>malloc(dim*dim*sizeof(cdouble)*order)
|
|
if not self.G[0].irreps[iri].m: raise MemoryError
|
|
if is1d:
|
|
for i in range(order):
|
|
self.G[0].irreps[iri].m[i] = irrep[permlist[i]]
|
|
else:
|
|
for i in range(order):
|
|
for row in range(dim):
|
|
for col in range(dim):
|
|
self.G[0].irreps[iri].m[i*dim*dim + row*dim + col] = irrep[permlist[i]][row,col]
|
|
self.G[0].elemlabels = <char **> 0 # Elem labels not yet implemented
|
|
self.owns_data = True
|
|
|
|
def __dealloc__(self):
|
|
cdef qpms_gmi_t order
|
|
if self.owns_data:
|
|
if self.G:
|
|
order = self.G[0].order
|
|
free(self.G[0].name)
|
|
free(self.G[0].mt)
|
|
free(self.G[0].invi)
|
|
free(self.G[0].gens)
|
|
if self.G[0].permrep:
|
|
for i in range(order): free(self.G[0].permrep[i])
|
|
free(self.G[0].permrep)
|
|
if self.G[0].elemlabels: # this is not even contructed right now
|
|
for i in range(order): free(self.G[0].elemlabels[i])
|
|
if self.G[0].irreps:
|
|
for iri in range(self.G[0].nirreps):
|
|
free(self.G[0].irreps[iri].name)
|
|
free(self.G[0].irreps[iri].m)
|
|
free(self.G[0].irreps)
|
|
free(self.G)
|
|
self.G = <qpms_finite_group_t *>0
|
|
self.owns_data = False
|
|
|
|
property order: # LPTODO might instead be __len__() if iterable at some point
|
|
def __get__(self):
|
|
return self.G[0].order
|
|
|
|
@staticmethod
|
|
def TRIVIAL(): # quite ugly
|
|
from .symmetries import point_group_info
|
|
return FinitePointGroup(point_group_info['trivial_g'])
|
|
|
|
|
|
cdef class FinitePointGroupElement:
|
|
'''TODO'''
|
|
cdef readonly FinitePointGroup G
|
|
cdef readonly qpms_gmi_t gmi
|
|
def __cinit__(self, FinitePointGroup G, qpms_gmi_t gmi):
|
|
self.G = G
|
|
self.gmi = gmi
|
|
|
|
cdef class Particle:
|
|
'''
|
|
Wrapper over the qpms_particle_t structure.
|
|
'''
|
|
cdef qpms_particle_t p
|
|
cdef readonly TMatrixFunction f # Reference to ensure correct reference counting
|
|
|
|
|
|
def __cinit__(Particle self, pos, t, bspec = None):
|
|
"""Particle object constructor.
|
|
|
|
Parameters
|
|
----------
|
|
pos : (x, y, z)
|
|
Particle position in cartesian coordinates
|
|
|
|
t : TMatrixGenerator or CTMatrix or TMatrixInterpolator
|
|
T-matrix specification for the particle.
|
|
|
|
bspec : BaseSpec
|
|
WSWF basis specification for the particle. Might be omitted if
|
|
t is a CTMatrix.
|
|
"""
|
|
cdef TMatrixGenerator tgen
|
|
cdef BaseSpec spec
|
|
if(len(pos)>=2 and len(pos) < 4):
|
|
self.p.pos.x = pos[0]
|
|
self.p.pos.y = pos[1]
|
|
self.p.pos.z = pos[2] if len(pos)==3 else 0
|
|
else:
|
|
raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates")
|
|
if isinstance(t, CTMatrix):
|
|
tgen = TMatrixGenerator(t)
|
|
elif isinstance(t, TMatrixInterpolator):
|
|
tgen = TMatrixGenerator(t)
|
|
warnings.warn("Initialising a particle with interpolated T-matrix values. Imaginary frequencies will be discarded and mode search algorithm will yield nonsense (just saying).")
|
|
elif isinstance(t, TMatrixGenerator):
|
|
tgen = <TMatrixGenerator>t
|
|
else: raise TypeError('t must be either CTMatrix or TMatrixGenerator, was %s' % str(type(t)))
|
|
if bspec is not None:
|
|
spec = bspec
|
|
else:
|
|
if isinstance(tgen.holder, CTMatrix):
|
|
spec = (<CTMatrix>tgen.holder).spec
|
|
else:
|
|
raise ValueError("bspec argument must be specified separately for %s" % str(type(t)))
|
|
self.f = TMatrixFunction(tgen, spec)
|
|
self.p.tmg = self.f.rawpointer()
|
|
# TODO non-trivial transformations later; if modified, do not forget to update ScatteringSystem constructor
|
|
self.p.op = qpms_tmatrix_operation_noop
|
|
|
|
def __dealloc__(self):
|
|
qpms_tmatrix_operation_clear(&self.p.op)
|
|
|
|
cdef qpms_particle_t *rawpointer(Particle self):
|
|
'''Pointer to the qpms_particle_p structure.
|
|
'''
|
|
return &(self.p)
|
|
property rawpointer:
|
|
def __get__(self):
|
|
return <uintptr_t> &(self.p)
|
|
|
|
cdef qpms_particle_t cval(Particle self):
|
|
'''Provides a copy for assigning in cython code'''
|
|
return self.p
|
|
|
|
property x:
|
|
def __get__(self):
|
|
return self.p.pos.x
|
|
def __set__(self,x):
|
|
self.p.pos.x = x
|
|
property y:
|
|
def __get__(self):
|
|
return self.p.pos.y
|
|
def __set__(self,y):
|
|
self.p.pos.y = y
|
|
property z:
|
|
def __get__(self):
|
|
return self.p.pos.z
|
|
def __set__(self,z):
|
|
self.p.pos.z = z
|
|
property pos:
|
|
def __get__(self):
|
|
return (self.p.pos.x, self.p.pos.y, self.p.pos.z)
|
|
def __set__(self, pos):
|
|
if(len(pos)>=2 and len(pos) < 4):
|
|
self.p.pos.x = pos[0]
|
|
self.p.pos.y = pos[1]
|
|
self.p.pos.z = pos[2] if len(pos)==3 else 0
|
|
else:
|
|
raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates")
|
|
|
|
cpdef void scatsystem_set_nthreads(long n):
|
|
qpms_scatsystem_set_nthreads(n)
|
|
return
|
|
|
|
|
|
cdef class ScatteringSystem:
|
|
'''
|
|
Wrapper over the C qpms_scatsys_t structure, representing a collection
|
|
of scatterers, finite or periodic.
|
|
|
|
Currently, it does not have a standard constructor. Use the
|
|
ScatteringSystem.create() method instead.
|
|
'''
|
|
cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators)
|
|
#cdef list Tmatrices # Here we keep the references to occuring T-matrices
|
|
cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator
|
|
cdef qpms_scatsys_t *s
|
|
cdef FinitePointGroup sym
|
|
|
|
cdef qpms_iri_t iri_py2c(self, iri, allow_None = True):
|
|
if iri is None and allow_None:
|
|
return QPMS_NO_IRREP
|
|
cdef qpms_iri_t nir = self.nirreps
|
|
cdef qpms_iri_t ciri = iri
|
|
if ciri < 0 or ciri > nir:
|
|
raise ValueError("Invalid irrep index %s (of %d irreps)", str(iri), self.nirreps)
|
|
return ciri
|
|
|
|
def check_s(self): # cdef instead?
|
|
if self.s == <qpms_scatsys_t *>NULL:
|
|
raise ValueError("ScatteringSystem's s-pointer not set. You must not use the default constructor; use the create() method instead")
|
|
#TODO is there a way to disable the constructor outside this module?
|
|
|
|
@staticmethod # We don't have any "standard" constructor for this right now
|
|
def create(particles, medium, cdouble omega, FinitePointGroup sym = FinitePointGroup.TRIVIAL(),
|
|
latticebasis = None): # TODO tolerances
|
|
"""(Non-standard) constructor of ScatteringSystem
|
|
|
|
Parameters
|
|
----------
|
|
particles : list of Particles objects
|
|
Scatterers to be included in the system. These scatterers are then
|
|
copied around by the constructor using the point group symmetry defined in sym.
|
|
medium : EpsMu or EpsMuGenerator
|
|
Material properties of the background medium.
|
|
omega : complex
|
|
Any valid angular frequency for the initial T-matrix evaluation.
|
|
It must be a value which all the T-matrix generators and interpolators
|
|
referenced by particles can evaluate.
|
|
sym : FinitePointGroup
|
|
Symmetry group for a finite system.
|
|
Defaults to the trivial point group FinitePointGroup.TRIVIAL().
|
|
Currently, this must be left trivial for periodic systems.
|
|
latticebasis : None or array_like of float type.
|
|
Lattice base vectors of a periodic system in cartesian coordinates.
|
|
The shape must be [d][3], where d = 1, 2 or 3 determines the lattice dimension.
|
|
Defaults to None, i.e. a finite system.
|
|
|
|
Returns
|
|
-------
|
|
ss : ScatteringSystem
|
|
Object representing the system of compact scatterers.
|
|
ssw : _ScatteringSystemAtOmega
|
|
Object representing ss evaluated at omega (angular frequency used for the
|
|
initial evaluation).
|
|
|
|
Note
|
|
----
|
|
Currently, the ScatterinSystem's T-matrices need to be evaluated for at least
|
|
one frequency in order to initialise the symmetry structure. For this reason,
|
|
this "constructor" creates an instance of ScatteringSystem and an instance
|
|
of _ScatteringSystemAtOmega at the same time, so that the evaluated T-matrices
|
|
can be used right away if needed. This is why this approach is used instead
|
|
of the usual __init__ method.
|
|
"""
|
|
|
|
if latticebasis is not None and sym.order != 1:
|
|
raise NotImplementedError("Periodic systems don't currently support nontrivial point group symmetries")
|
|
|
|
# These we are going to construct
|
|
cdef ScatteringSystem self
|
|
cdef _ScatteringSystemAtOmega pyssw
|
|
|
|
cdef qpms_scatsys_t orig # This should be automatically init'd to 0 (CHECKME)
|
|
cdef qpms_ss_pi_t pi, p_count = len(particles)
|
|
cdef qpms_ss_tmi_t tmi, tm_count = 0
|
|
cdef qpms_ss_tmgi_t tmgi, tmg_count = 0
|
|
|
|
cdef qpms_scatsys_at_omega_t *ssw
|
|
cdef qpms_scatsys_t *ss
|
|
|
|
cdef Particle p
|
|
|
|
tmgindices = dict()
|
|
tmgobjs = list()
|
|
tmindices = dict()
|
|
tmlist = list()
|
|
for p in particles: # find and enumerate unique t-matrix generators
|
|
if p.p.op.typ != QPMS_TMATRIX_OPERATION_NOOP:
|
|
raise NotImplementedError("currently, only no-op T-matrix operations are allowed in ScatteringSystem constructor")
|
|
#tmg_key = id(p.f) # This causes a different generator for each particle -> SUPER SLOW
|
|
tmg_key = (id(p.f.generator), id(p.f.spec))
|
|
if tmg_key not in tmgindices:
|
|
tmgindices[tmg_key] = tmg_count
|
|
tmgobjs.append(p.f) # Save the references on BaseSpecs and TMatrixGenerators (via TMatrixFunctions)
|
|
tmg_count += 1
|
|
# Following lines have to be adjusted when nontrivial operations allowed:
|
|
tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None
|
|
if tm_derived_key not in tmindices:
|
|
tmindices[tm_derived_key] = tm_count
|
|
tmlist.append(tm_derived_key)
|
|
tm_count += 1
|
|
cdef EpsMuGenerator mediumgen = EpsMuGenerator(medium)
|
|
orig.medium = mediumgen.g
|
|
orig.tmg_count = tmg_count
|
|
orig.tm_count = tm_count
|
|
orig.p_count = p_count
|
|
try:
|
|
orig.tmg = <qpms_tmatrix_function_t *>malloc(orig.tmg_count * sizeof(orig.tmg[0]))
|
|
if not orig.tmg: raise MemoryError
|
|
orig.tm = <qpms_ss_derived_tmatrix_t *>malloc(orig.tm_count * sizeof(orig.tm[0]))
|
|
if not orig.tm: raise MemoryError
|
|
orig.p = <qpms_particle_tid_t *>malloc(orig.p_count * sizeof(orig.p[0]))
|
|
if not orig.p: raise MemoryError
|
|
for tmgi in range(orig.tmg_count):
|
|
orig.tmg[tmgi] = (<TMatrixFunction?>tmgobjs[tmgi]).raw()
|
|
for tmi in range(tm_count):
|
|
tm_derived_key = tmlist[tmi]
|
|
tmgi = tmgindices[tm_derived_key[0]]
|
|
orig.tm[tmi].tmgi = tmgi
|
|
orig.tm[tmi].op = qpms_tmatrix_operation_noop # TODO adjust when notrivial operations allowed
|
|
for pi in range(p_count):
|
|
p = particles[pi]
|
|
tmg_key = (id(p.f.generator), id(p.f.spec))
|
|
tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None
|
|
orig.p[pi].pos = p.cval().pos
|
|
orig.p[pi].tmatrix_id = tmindices[tm_derived_key]
|
|
if latticebasis is not None: # periodic system
|
|
assert(len(latticebasis) <= 3 and len(latticebasis) > 0)
|
|
orig.lattice_dimension = len(latticebasis)
|
|
for d in range(len(latticebasis)):
|
|
orig.per.lattice_basis[d] = {'x' : latticebasis[d][0], 'y' : latticebasis[d][1], 'z' : latticebasis[d][2] if len(latticebasis[d]) >= 3 else 0}
|
|
else: orig.lattice_dimension = 0
|
|
ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT)
|
|
ss = ssw[0].ss
|
|
finally:
|
|
free(orig.tmg)
|
|
free(orig.tm)
|
|
free(orig.p)
|
|
self = ScatteringSystem()
|
|
self.medium_holder = mediumgen
|
|
self.s = ss
|
|
self.tmgobjs = tmgobjs
|
|
self.sym = sym
|
|
pyssw = _ScatteringSystemAtOmega()
|
|
pyssw.ssw = ssw
|
|
pyssw.ss_pyref = self
|
|
return self, pyssw
|
|
|
|
def __call__(self, cdouble omega):
|
|
self.check_s()
|
|
cdef _ScatteringSystemAtOmega pyssw = _ScatteringSystemAtOmega()
|
|
pyssw.ssw = qpms_scatsys_at_omega(self.s, omega)
|
|
pyssw.ss_pyref = self
|
|
return pyssw
|
|
|
|
def __dealloc__(self):
|
|
if(self.s):
|
|
qpms_scatsys_free(self.s)
|
|
|
|
property particles_tmi:
|
|
def __get__(self):
|
|
self.check_s()
|
|
r = list()
|
|
cdef qpms_ss_pi_t pi
|
|
for pi in range(self.s[0].p_count):
|
|
r.append(self.s[0].p[pi])
|
|
return r
|
|
|
|
property particle_count:
|
|
"""Number of particles in the system (or unit cell)"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
return self.s[0].p_count
|
|
property fecv_size:
|
|
"""Length of the full excitation coefficient vector"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
return self.s[0].fecv_size
|
|
property saecv_sizes:
|
|
"""Lengths of the partial symmetry-adapted excitation coefficient vectors"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
return [self.s[0].saecv_sizes[i]
|
|
for i in range(self.s[0].sym[0].nirreps)]
|
|
property irrep_names:
|
|
def __get__(self):
|
|
self.check_s()
|
|
return [string_c2py(self.s[0].sym[0].irreps[iri].name)
|
|
if (self.s[0].sym[0].irreps[iri].name) else None
|
|
for iri in range(self.s[0].sym[0].nirreps)]
|
|
property nirreps:
|
|
"""Number of irreducible representations of the scattering system point symmetry group"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
return self.s[0].sym[0].nirreps
|
|
property lattice_dimension:
|
|
def __get__(self):
|
|
return self.s[0].lattice_dimension
|
|
|
|
property lattice_basis:
|
|
"""Direct lattice basis vectors for periodic system in 3D cartesian coordinates"""
|
|
def __get__(self):
|
|
cdef int d = self.lattice_dimension
|
|
if d == 0:
|
|
return None
|
|
else:
|
|
return np.array([[self.s[0].per.lattice_basis[i].x,
|
|
self.s[0].per.lattice_basis[i].y,
|
|
self.s[0].per.lattice_basis[i].z,
|
|
] for i in range(d)])
|
|
|
|
property reciprocal_basis:
|
|
"""Reciprocal lattice basis vectors for periodic system in 3D cartesian coordinates"""
|
|
def __get__(self):
|
|
cdef int d = self.lattice_dimension
|
|
if d == 0:
|
|
return None
|
|
else:
|
|
return np.array([[self.s[0].per.reciprocal_basis[i].x,
|
|
self.s[0].per.reciprocal_basis[i].y,
|
|
self.s[0].per.reciprocal_basis[i].z,
|
|
] for i in range(d)])
|
|
|
|
def bspec_pi(self, qpms_ss_pi_t pi):
|
|
"""Gives the BaseSpec of a particle defined by an index in the internal particles' order, in the form of an array that can be fed to BaseSpec constructor"""
|
|
self.check_s()
|
|
cdef qpms_ss_pi_t pcount = self.s[0].p_count
|
|
if pi >= pcount or pi < 0: raise IndexError("Invalid particle index")
|
|
cdef const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(self.s, pi)
|
|
return np.array([ bspec[0].ilist[i] for i in range(bspec[0].n)])
|
|
property bspecs:
|
|
"""Gives the BaseSpecs of all particles in the internal particles' order, in the form of arrays that can be fed to BaseSpec constructor"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
cdef qpms_ss_pi_t pi, pcount = self.s[0].p_count
|
|
return [self.bspec_pi(pi) for pi in range(pcount)]
|
|
|
|
property unitcell_volume:
|
|
def __get__(self):
|
|
self.check_s()
|
|
if self.lattice_dimension:
|
|
return self.s[0].per.unitcell_volume
|
|
else:
|
|
return None
|
|
|
|
property eta:
|
|
"""Ewald parameter η (only relevant for periodic systems)
|
|
|
|
This is the default value that will be copied into
|
|
_ScatteringSystemAtOmega by __call__"""
|
|
def __get__(self):
|
|
self.check_s()
|
|
if self.lattice_dimension:
|
|
return self.s[0].per.eta
|
|
else:
|
|
return None
|
|
|
|
def __set__(self, eta):
|
|
self.check_s()
|
|
if self.lattice_dimension:
|
|
self.s[0].per.eta = eta
|
|
else:
|
|
raise AttributeError("Cannot set Ewald parameter for finite system") # different exception?
|
|
|
|
|
|
def pack_vector(self, vect, iri):
|
|
"""Converts (projects) a full excitation coefficient vector into an irrep subspace.
|
|
|
|
Parameters
|
|
----------
|
|
vect : array_like of shape (self.fecv_size,)
|
|
The full excitation coefficient vector to be converted.
|
|
iri : int
|
|
Index of the irreducible representation.
|
|
|
|
Returns
|
|
-------
|
|
packed : ndarray of shape (self.saecv_sizes[iri],)
|
|
"irrep-packed" excitation vector: the part of vect belonging to the
|
|
irrep indexed by iri, in terms of given irrep's internal coordinates.
|
|
|
|
See Also
|
|
--------
|
|
unpack_vector : The beckward (not exactly inverse) operation.
|
|
pack_matrix : Corresponding conversion for matrices.
|
|
"""
|
|
self.check_s()
|
|
if len(vect) != self.fecv_size:
|
|
raise ValueError("Length of a full vector has to be %d, not %d"
|
|
% (self.fecv_size, len(vect)))
|
|
vect = np.array(vect, dtype=complex, copy=False, order='C')
|
|
cdef cdouble[::1] vect_view = vect;
|
|
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
|
|
(self.saecv_sizes[iri],), dtype=complex, order='C')
|
|
cdef cdouble[::1] target_view = target_np
|
|
qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri)
|
|
return target_np
|
|
|
|
def unpack_vector(self, packed, iri):
|
|
"""Unpacks an "irrep-packed" excitation coefficient vector to full coordinates.
|
|
|
|
Parameters
|
|
----------
|
|
packed : array_like of shape (self.saecv_sizes[iri],)
|
|
Excitation coefficient vector part belonging to the irreducible representation
|
|
indexed by iri, in terms of given irrep's internal coordinates.
|
|
iri : int
|
|
Index of the irreducible representation.
|
|
|
|
Returns
|
|
-------
|
|
vect : ndarray of shape (self.fecv_size,)
|
|
The contribution to a full excitation coefficient vector from iri'th irrep.
|
|
|
|
See Also
|
|
--------
|
|
pack_vector : The inverse operation.
|
|
unpack_matrix : Corresponding conversion for matrices.
|
|
"""
|
|
self.check_s()
|
|
if len(packed) != self.saecv_sizes[iri]:
|
|
raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d"
|
|
% (iri, self.saecv_sizes, len(packed)))
|
|
packed = np.array(packed, dtype=complex, copy=False, order='C')
|
|
cdef cdouble[::1] packed_view = packed
|
|
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
|
|
(self.fecv_size,), dtype=complex)
|
|
cdef cdouble[::1] target_view = target_np
|
|
qpms_scatsys_irrep_unpack_vector(&target_view[0], &packed_view[0],
|
|
self.s, iri, 0)
|
|
return target_np
|
|
|
|
def pack_matrix(self, fullmatrix, iri):
|
|
"""Converts (projects) a matrix into an irrep subspace.
|
|
|
|
Parameters
|
|
----------
|
|
fullmatrix : array_like of shape (self.fecv_size, self.fecv_size)
|
|
The full matrix (operating on the exctitation coefficient vectors) to be converted.
|
|
iri : int
|
|
Index of the irreducible representation.
|
|
|
|
Returns
|
|
-------
|
|
packedmatrix : ndarray of shape (self.saecv_sizes[iri], self.saecv_sizes[iri])
|
|
"irrep-packed" matrix: the part of fullmatrix belonging to the
|
|
irrep indexed by iri, in terms of given irrep's internal coordinates.
|
|
|
|
See Also
|
|
--------
|
|
unpack_matrix : The beckward (not exactly inverse) operation.
|
|
pack_vector : Corresponding conversion for excitation coefficient vectors.
|
|
"""
|
|
self.check_s()
|
|
cdef size_t flen = self.s[0].fecv_size
|
|
cdef size_t rlen = self.saecv_sizes[iri]
|
|
fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C')
|
|
if fullmatrix.shape != (flen, flen):
|
|
raise ValueError("Full matrix shape should be (%d,%d), is %s."
|
|
% (flen, flen, repr(fullmatrix.shape)))
|
|
cdef cdouble[:,::1] fullmatrix_view = fullmatrix
|
|
cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
|
|
(rlen, rlen), dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target_np
|
|
qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
|
|
self.s, iri)
|
|
return target_np
|
|
def unpack_matrix(self, packedmatrix, iri):
|
|
"""Unpacks an "irrep-packed" excitation coefficient vector to full coordinates.
|
|
|
|
Parameters
|
|
----------
|
|
packedmatrix : array_like of shape (self.saecv_sizes[iri], self.saecv_sizes[iri])
|
|
A matrix (operating on the excitation coefficient vectors) part belonging to the
|
|
irreducible representation indexed by iri, in terms of given irrep's internal
|
|
coordinates.
|
|
iri : int
|
|
Index of the irreducible representation.
|
|
|
|
Returns
|
|
-------
|
|
fullmatrix : ndarray of shape (self.fecv_size, self.fecv_size)
|
|
The iri'th irrep contribution to a full matrix.
|
|
|
|
See Also
|
|
--------
|
|
pack_matrix : The inverse operation.
|
|
unpack_vector : Corresponding conversion for excitation coefficient vectors.
|
|
"""
|
|
self.check_s()
|
|
cdef size_t flen = self.s[0].fecv_size
|
|
cdef size_t rlen = self.saecv_sizes[iri]
|
|
packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C')
|
|
if packedmatrix.shape != (rlen, rlen):
|
|
raise ValueError("Packed matrix shape should be (%d,%d), is %s."
|
|
% (rlen, rlen, repr(packedmatrix.shape)))
|
|
cdef cdouble[:,::1] packedmatrix_view = packedmatrix
|
|
cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
|
|
(flen, flen), dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target_np
|
|
qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
|
|
self.s, iri, 0)
|
|
return target_np
|
|
|
|
def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS, double eta = 0):
|
|
"""Constructs full translation matrix of a scattering system.
|
|
|
|
This method enables to use any wave number for the background medium ignoring the
|
|
background EpsMuGenerator), using only system's particle positions (and lattice
|
|
basis for infinite system).
|
|
|
|
Parameters
|
|
----------
|
|
wavenumber : complex
|
|
Wave number of the medium
|
|
|
|
blochvector : array_like of shape (3,)
|
|
Bloch vector (only for periodic systems)
|
|
|
|
J : BesselType
|
|
Optionally, one can replace Hankel functions of the first kind with different
|
|
Bessel functions.
|
|
|
|
eta : float
|
|
Ewald parameter η, only relevant for periodic lattices;
|
|
if set to 0 (default), the actual value is determined
|
|
automatically.
|
|
|
|
See Also
|
|
--------
|
|
ScatteringSystemAtOmega.translation_matrix_full : Translation matrix at a given frequency.
|
|
"""
|
|
self.check_s()
|
|
cdef size_t flen = self.s[0].fecv_size
|
|
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
|
(flen,flen),dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target
|
|
cdef cart3_t blochvector_c
|
|
if self.lattice_dimension == 0:
|
|
if blochvector is None:
|
|
qpms_scatsys_build_translation_matrix_e_full(&target_view[0][0], self.s, wavenumber, J)
|
|
else: raise ValueError("Can't use blochvector with non-periodic system")
|
|
else:
|
|
if blochvector is None: raise ValueError("Valid blochvector must be specified for periodic system")
|
|
else:
|
|
if J != QPMS_HANKEL_PLUS:
|
|
raise NotImplementedError("Translation operators based on other than Hankel+ functions not supperted in periodic systems")
|
|
blochvector_c = {'x': blochvector[0], 'y': blochvector[1], 'z': blochvector[2]}
|
|
with pgsl_ignore_error(15):
|
|
qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c, eta)
|
|
return target
|
|
|
|
def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
|
|
self.check_s()
|
|
cdef size_t rlen = self.saecv_sizes[iri]
|
|
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
|
(rlen,rlen),dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target
|
|
qpms_scatsys_build_translation_matrix_e_irrep_packed(&target_view[0][0],
|
|
self.s, iri, wavenumber, J)
|
|
return target
|
|
|
|
property fullvec_psizes:
|
|
def __get__(self):
|
|
self.check_s()
|
|
cdef np.ndarray[int32_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.int32)
|
|
cdef int32_t[::1] ar_view = ar
|
|
for pi in range(self.s[0].p_count):
|
|
ar_view[pi] = qpms_ss_bspec_pi(self.s, pi)[0].n
|
|
return ar
|
|
|
|
|
|
property fullvec_poffsets:
|
|
def __get__(self):
|
|
self.check_s()
|
|
cdef np.ndarray[intptr_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.intp)
|
|
cdef intptr_t[::1] ar_view = ar
|
|
cdef intptr_t offset = 0
|
|
for pi in range(self.s[0].p_count):
|
|
ar_view[pi] = offset
|
|
offset += qpms_ss_bspec_pi(self.s, pi)[0].n
|
|
return ar
|
|
|
|
property positions:
|
|
def __get__(self):
|
|
self.check_s()
|
|
cdef np.ndarray[np.double_t, ndim=2] ar = np.empty((self.s[0].p_count, 3), dtype=float)
|
|
cdef np.double_t[:,::1] ar_view = ar
|
|
for pi in range(self.s[0].p_count):
|
|
ar_view[pi,0] = self.s[0].p[pi].pos.x
|
|
ar_view[pi,1] = self.s[0].p[pi].pos.y
|
|
ar_view[pi,2] = self.s[0].p[pi].pos.z
|
|
return ar
|
|
|
|
def planewave_full(self, k_cart, E_cart):
|
|
self.check_s()
|
|
k_cart = np.array(k_cart)
|
|
E_cart = np.array(E_cart)
|
|
if k_cart.shape != (3,) or E_cart.shape != (3,):
|
|
raise ValueError("k_cart and E_cart must be ndarrays of shape (3,)")
|
|
cdef qpms_incfield_planewave_params_t p
|
|
p.use_cartesian = 1
|
|
p.k.cart.x = <cdouble>k_cart[0]
|
|
p.k.cart.y = <cdouble>k_cart[1]
|
|
p.k.cart.z = <cdouble>k_cart[2]
|
|
p.E.cart.x = <cdouble>E_cart[0]
|
|
p.E.cart.y = <cdouble>E_cart[1]
|
|
p.E.cart.z = <cdouble>E_cart[2]
|
|
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
|
|
(self.fecv_size,), dtype=complex)
|
|
cdef cdouble[::1] target_view = target_np
|
|
qpms_scatsys_incident_field_vector_full(&target_view[0],
|
|
self.s, qpms_incfield_planewave, <void *>&p, 0)
|
|
return target_np
|
|
|
|
def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri, iri = None,
|
|
blochvector = None,
|
|
size_t contour_points = 20, double rank_tol = 1e-4, size_t rank_min_sel=1,
|
|
double res_tol = 0):
|
|
"""
|
|
Attempts to find the eigenvalues and eigenvectors using Beyn's algorithm.
|
|
"""
|
|
|
|
cdef beyn_result_t *res
|
|
cdef double blochvector_c[3]
|
|
if self.lattice_dimension == 0:
|
|
assert(blochvector is None)
|
|
res = qpms_scatsys_finite_find_eigenmodes(self.s, self.iri_py2c(iri),
|
|
omega_centre, omega_rr, omega_ri, contour_points,
|
|
rank_tol, rank_min_sel, res_tol)
|
|
else:
|
|
if iri is not None: raise NotImplementedError("Irreps decomposition not yet supported for periodic systems")
|
|
blochvector_c[0] = blochvector[0]
|
|
blochvector_c[1] = blochvector[1]
|
|
blochvector_c[2] = blochvector[2] if len(blochvector) > 2 else 0
|
|
res = qpms_scatsys_periodic_find_eigenmodes(self.s, blochvector_c,
|
|
omega_centre, omega_rr, omega_ri, contour_points, rank_tol, rank_min_sel, res_tol)
|
|
if res == NULL: raise RuntimeError
|
|
|
|
cdef size_t neig = res[0].neig, i, j
|
|
cdef size_t vlen = res[0].vlen # should be equal to self.s.fecv_size
|
|
|
|
cdef np.ndarray[complex, ndim=1] eigval = np.empty((neig,), dtype=complex)
|
|
cdef cdouble[::1] eigval_v = eigval
|
|
cdef np.ndarray[complex, ndim=1] eigval_err = np.empty((neig,), dtype=complex)
|
|
cdef cdouble[::1] eigval_err_v = eigval_err
|
|
cdef np.ndarray[double, ndim=1] residuals = np.empty((neig,), dtype=np.double)
|
|
cdef double[::1] residuals_v = residuals
|
|
cdef np.ndarray[complex, ndim=2] eigvec = np.empty((neig,vlen),dtype=complex)
|
|
cdef cdouble[:,::1] eigvec_v = eigvec
|
|
cdef np.ndarray[double, ndim=1] ranktest_SV = np.empty((vlen), dtype=np.double)
|
|
cdef double[::1] ranktest_SV_v = ranktest_SV
|
|
|
|
for i in range(neig):
|
|
eigval_v[i] = res[0].eigval[i]
|
|
eigval_err_v[i] = res[0].eigval_err[i]
|
|
residuals_v[i] = res[0].residuals[i]
|
|
for j in range(vlen):
|
|
eigvec_v[i,j] = res[0].eigvec[i*vlen + j]
|
|
for i in range(vlen):
|
|
ranktest_SV_v[i] = res[0].ranktest_SV[i]
|
|
|
|
zdist = eigval - omega_centre
|
|
eigval_inside_metric = np.hypot(zdist.real / omega_rr, zdist.imag / omega_ri)
|
|
|
|
beyn_result_free(res)
|
|
retdict = {
|
|
'eigval':eigval,
|
|
'eigval_inside_metric':eigval_inside_metric,
|
|
'eigvec':eigvec,
|
|
'residuals':residuals,
|
|
'eigval_err':eigval_err,
|
|
'ranktest_SV':ranktest_SV,
|
|
'iri': iri,
|
|
'blochvector': blochvector
|
|
}
|
|
|
|
return retdict
|
|
|
|
@boundscheck(False)
|
|
def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, btyp=BesselType.HANKEL_PLUS, bint alt=False):
|
|
# TODO DOC, periodic systems
|
|
if self.lattice_dimension != 0: # TODO enable also periodic systems
|
|
raise NotImplementedError("For periodic system, use _ScatteringSystemAtOmegaK.scattered_E")
|
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
|
if evalpos.shape[-1] != 3:
|
|
raise ValueError("Last dimension of evalpos has to be 3")
|
|
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
|
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False, dtype=complex)
|
|
cdef cdouble[::1] scv_view = scv
|
|
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
|
|
cdef ccart3_t res
|
|
cdef cart3_t pos
|
|
cdef Py_ssize_t i
|
|
with nogil, parallel(), boundscheck(False), wraparound(False):
|
|
for i in prange(evalpos_a.shape[0]):
|
|
pos.x = evalpos_a[i,0]
|
|
pos.y = evalpos_a[i,1]
|
|
pos.z = evalpos_a[i,2]
|
|
if alt:
|
|
res = qpms_scatsys_scattered_E__alt(self.s, btyp_c, wavenumber, &scv_view[0], pos)
|
|
else:
|
|
res = qpms_scatsys_scattered_E(self.s, btyp_c, wavenumber, &scv_view[0], pos)
|
|
results[i,0] = res.x
|
|
results[i,1] = res.y
|
|
results[i,2] = res.z
|
|
return results.reshape(evalpos.shape)
|
|
|
|
def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double maxomega):
|
|
'''Empty (2D, xy-plane) lattice mode (diffraction order) frequencies of a non-dispersive medium.
|
|
|
|
reciprocal_basis is of the "mutliplied by 2п" type.
|
|
'''
|
|
cdef double *omegas_c
|
|
cdef size_t n
|
|
cdef cart2_t k_, b1, b2
|
|
k_.x = wavevector[0]
|
|
k_.y = wavevector[1]
|
|
b1.x = reciprocal_basis[0][0]
|
|
b1.y = reciprocal_basis[0][1]
|
|
b2.x = reciprocal_basis[1][0]
|
|
b2.y = reciprocal_basis[1][1]
|
|
if(epsmu.n.imag != 0):
|
|
warnings.warn("Got complex refractive index", epsmu.n, "ignoring the imaginary part")
|
|
refindex = epsmu.n.real
|
|
n = qpms_emptylattice2_modes_maxfreq(&omegas_c, b1, b2, BASIS_RTOL,
|
|
k_, GSL_CONST_MKSA_SPEED_OF_LIGHT / refindex, maxomega)
|
|
cdef np.ndarray[double, ndim=1] omegas = np.empty((n,), dtype=np.double)
|
|
cdef double[::1] omegas_v = omegas
|
|
cdef size_t i
|
|
for i in range(n):
|
|
omegas_v[i] = omegas_c[i]
|
|
free(omegas_c)
|
|
return omegas
|
|
|
|
cdef class _ScatteringSystemAtOmegaK:
|
|
'''
|
|
Wrapper over the C qpms_scatsys_at_omega_k_t structure
|
|
|
|
This represents an infinite periodic system with a given (fixed) frequency and Bloch vector.
|
|
'''
|
|
cdef qpms_scatsys_at_omega_k_t sswk
|
|
cdef _ScatteringSystemAtOmega ssw_pyref
|
|
|
|
cdef qpms_scatsys_at_omega_k_t *rawpointer(self):
|
|
return &self.sswk
|
|
|
|
property eta:
|
|
"""Ewald parameter η"""
|
|
def __get__(self):
|
|
return self.sswk.eta
|
|
def __set__(self, double eta):
|
|
self.sswk.eta = eta
|
|
|
|
@boundscheck(False)
|
|
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
|
|
"""Evaluate electric field for a given excitation coefficient vector (periodic system)
|
|
|
|
Parameters
|
|
----------
|
|
scatcoeffvector_full: array_like of shape (self.fecv_size,)
|
|
Onedimensional excitation coefficient vector, describing SVWFs leaving
|
|
the scatterers inside the representative unit cell.
|
|
evalpos: array_like of floats and shape (..., 3)
|
|
Evaluation points in cartesian coordinates.
|
|
|
|
Returns
|
|
-------
|
|
ndarray of complex, with the same shape as `evalpos`
|
|
Electric field at the positions given in `evalpos`, in cartesian coordinates.
|
|
"""
|
|
if(btyp != QPMS_HANKEL_PLUS):
|
|
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
|
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
|
if evalpos.shape[-1] != 3:
|
|
raise ValueError("Last dimension of evalpos has to be 3")
|
|
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
|
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, dtype=complex, copy=False)
|
|
cdef cdouble[::1] scv_view = scv
|
|
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
|
|
cdef ccart3_t res
|
|
cdef cart3_t pos
|
|
cdef Py_ssize_t i
|
|
with nogil, wraparound(False), parallel():
|
|
for i in prange(evalpos_a.shape[0]):
|
|
pos.x = evalpos_a[i,0]
|
|
pos.y = evalpos_a[i,1]
|
|
pos.z = evalpos_a[i,2]
|
|
res = qpms_scatsyswk_scattered_E(&self.sswk, btyp_c, &scv_view[0], pos)
|
|
results[i,0] = res.x
|
|
results[i,1] = res.y
|
|
results[i,2] = res.z
|
|
return results.reshape(evalpos.shape)
|
|
|
|
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS):
|
|
# TODO examples
|
|
"""Evaluate scattered field "basis" (periodic system)
|
|
|
|
This function enables the evaluation of "scattered" fields
|
|
generated by the system for many different excitation
|
|
coefficients vectors, without the expensive re-evaluation of the respective
|
|
translation operators for each excitation coefficient vector.
|
|
|
|
Parameters
|
|
----------
|
|
evalpos: array_like of floats and shape (..., 3)
|
|
Evaluation points in cartesian coordinates.
|
|
|
|
Returns
|
|
-------
|
|
ndarray of complex, with the shape `evalpos.shape[:-1] + (self.fecv_size, 3)`
|
|
"Basis" fields at the positions given in `evalpos`, in cartesian coordinates.
|
|
"""
|
|
if(btyp != QPMS_HANKEL_PLUS):
|
|
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
|
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
|
cdef Py_ssize_t fecv_size = self.fecv_size
|
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
|
if evalpos.shape[-1] != 3:
|
|
raise ValueError("Last dimension of evalpos has to be 3")
|
|
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
|
cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex)
|
|
cdef ccart3_t *res
|
|
res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
|
|
cdef cart3_t pos
|
|
cdef Py_ssize_t i, j
|
|
with nogil, wraparound(False), parallel():
|
|
for i in prange(evalpos_a.shape[0]):
|
|
pos.x = evalpos_a[i,0]
|
|
pos.y = evalpos_a[i,1]
|
|
pos.z = evalpos_a[i,2]
|
|
qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos)
|
|
for j in range(fecv_size):
|
|
results[i,j,0] = res[j].x
|
|
results[i,j,1] = res[j].y
|
|
results[i,j,2] = res[j].z
|
|
free(res)
|
|
return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
|
|
|
|
property fecv_size:
|
|
def __get__(self): return self.ssw_pyref.ss_pyref.fecv_size
|
|
|
|
cdef class _ScatteringSystemAtOmega:
|
|
'''
|
|
Wrapper over the C qpms_scatsys_at_omega_t structure
|
|
that keeps the T-matrix and background data evaluated
|
|
at specific frequency.
|
|
'''
|
|
cdef qpms_scatsys_at_omega_t *ssw
|
|
cdef ScatteringSystem ss_pyref
|
|
|
|
def check(self): # cdef instead?
|
|
if not self.ssw:
|
|
raise ValueError("_ScatteringSystemAtOmega's ssw-pointer not set. You must not use the default constructor; ScatteringSystem.create() instead")
|
|
self.ss_pyref.check_s()
|
|
#TODO is there a way to disable the constructor outside this module?
|
|
|
|
|
|
def ensure_finite(self):
|
|
if self.ssw[0].ss[0].lattice_dimension != 0:
|
|
raise NotImplementedError("Operation not supported for periodic systems")
|
|
|
|
def __dealloc__(self):
|
|
if (self.ssw):
|
|
qpms_scatsys_at_omega_free(self.ssw)
|
|
|
|
def apply_Tmatrices_full(self, a):
|
|
self.check()
|
|
if len(a) != self.fecv_size:
|
|
raise ValueError("Length of a full vector has to be %d, not %d"
|
|
% (self.fecv_size, len(a)))
|
|
a = np.array(a, dtype=complex, copy=False, order='C')
|
|
cdef cdouble[::1] a_view = a;
|
|
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
|
|
(self.fecv_size,), dtype=complex, order='C')
|
|
cdef cdouble[::1] target_view = target_np
|
|
qpms_scatsysw_apply_Tmatrices_full(&target_view[0], &a_view[0], self.ssw)
|
|
return target_np
|
|
|
|
cdef qpms_scatsys_at_omega_t *rawpointer(self):
|
|
return self.ssw
|
|
|
|
def scatter_solver(self, iri=None, k=None):
|
|
self.check()
|
|
cdef _ScatteringSystemAtOmegaK sswk # used only for periodic systems
|
|
if(self.ssw[0].ss[0].lattice_dimension == 0):
|
|
return ScatteringMatrix(self, iri=iri)
|
|
else:
|
|
if iri is not None:
|
|
raise NotImplementedError("Irrep decomposition not (yet) supported for periodic systems")
|
|
sswk = self._sswk(k)
|
|
return ScatteringMatrix(ssw=self, sswk=sswk, iri=None)
|
|
|
|
def _sswk(self, k):
|
|
cdef _ScatteringSystemAtOmegaK sswk = _ScatteringSystemAtOmegaK()
|
|
sswk.sswk.ssw = self.ssw
|
|
sswk.ssw_pyref=self
|
|
sswk.sswk.k[0] = k[0]
|
|
sswk.sswk.k[1] = k[1]
|
|
sswk.sswk.k[2] = k[2]
|
|
sswk.eta = qpms_ss_adjusted_eta(self.ssw[0].ss, self.ssw[0].wavenumber, sswk.sswk.k)
|
|
return sswk
|
|
|
|
property fecv_size:
|
|
def __get__(self): return self.ss_pyref.fecv_size
|
|
property saecv_sizes:
|
|
def __get__(self): return self.ss_pyref.saecv_sizes
|
|
property irrep_names:
|
|
def __get__(self): return self.ss_pyref.irrep_names
|
|
property nirreps:
|
|
def __get__(self): return self.ss_pyref.nirreps
|
|
property wavenumber:
|
|
def __get__(self): return self.ssw[0].wavenumber
|
|
|
|
|
|
def modeproblem_matrix_full(self, k=None):
|
|
self.check()
|
|
cdef size_t flen = self.ss_pyref.s[0].fecv_size
|
|
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
|
(flen,flen),dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target
|
|
cdef _ScatteringSystemAtOmegaK sswk
|
|
if k is not None:
|
|
sswk = self._sswk(k)
|
|
qpms_scatsyswk_build_modeproblem_matrix_full(&target_view[0][0], &sswk.sswk)
|
|
else:
|
|
qpms_scatsysw_build_modeproblem_matrix_full(&target_view[0][0], self.ssw)
|
|
return target
|
|
|
|
def modeproblem_matrix_packed(self, qpms_iri_t iri, version='pR'):
|
|
self.check()
|
|
self.ensure_finite()
|
|
cdef size_t rlen = self.saecv_sizes[iri]
|
|
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
|
(rlen,rlen),dtype=complex, order='C')
|
|
cdef cdouble[:,::1] target_view = target
|
|
if (version == 'R'):
|
|
qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(&target_view[0][0], self.ssw, iri)
|
|
elif (version == 'pR'):
|
|
with nogil:
|
|
qpms_scatsysw_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.ssw, iri)
|
|
else:
|
|
qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.ssw, iri)
|
|
return target
|
|
|
|
def translation_matrix_full(self, blochvector = None):
|
|
"""Constructs full translation matrix of a scattering system at given frequency.
|
|
|
|
Parameters
|
|
----------
|
|
blochvector : array_like of shape (3,)
|
|
Bloch vector (only for periodic systems)
|
|
|
|
See Also
|
|
--------
|
|
ScatteringSystem.translation_matrix_full: translation matrix for any wavenumber
|
|
|
|
"""
|
|
return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector)
|
|
|
|
def translation_matrix_packed(self, iri, J = QPMS_HANKEL_PLUS):
|
|
return self.ss_pyref.translation_matrix_packed(wavenumber=self.wavenumber, iri=iri, J=J)
|
|
|
|
@boundscheck(False)
|
|
def scattered_E(self, scatcoeffvector_full, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS, bint alt=False):
|
|
# FIXME TODO this obviously does not work for periodic systems
|
|
"""Evaluate electric field for a given excitation coefficient vector
|
|
|
|
Parameters
|
|
----------
|
|
scatcoeffvector_full: array_like of shape (self.fecv_size,)
|
|
Onedimensional excitation coefficient vector, describing SVWFs leaving
|
|
the scatterers.
|
|
evalpos: array_like of floats and shape (..., 3)
|
|
Evaluation points in cartesian coordinates.
|
|
blochvector: array_like or None
|
|
Bloch vector, must be supplied (non-None) for periodic systems, else None.
|
|
|
|
Returns
|
|
-------
|
|
array_like of complex, with the same shape as `evalpos`
|
|
Electric field at the positions given in `evalpos`.
|
|
|
|
See Also
|
|
--------
|
|
_ScatteringSystemAtOmegaK.scattered_E: Evaluate scattered field in an infinite periodic system.
|
|
"""
|
|
if(self.ssw[0].ss[0].lattice_dimension != 0):
|
|
if blochvector is None:
|
|
raise ValueError("For a periodic system, blochvector must be specified.")
|
|
return self._sswk(blochvector).scattered_E(scatcoeffvector_full, evalpos=evalpos, btyp=btyp)
|
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
|
if evalpos.shape[-1] != 3:
|
|
raise ValueError("Last dimension of evalpos has to be 3")
|
|
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
|
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, dtype=complex, copy=False)
|
|
cdef cdouble[::1] scv_view = scv
|
|
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
|
|
cdef ccart3_t res
|
|
cdef cart3_t pos
|
|
cdef Py_ssize_t i
|
|
with wraparound(False), nogil, parallel():
|
|
for i in prange(evalpos_a.shape[0]):
|
|
pos.x = evalpos_a[i,0]
|
|
pos.y = evalpos_a[i,1]
|
|
pos.z = evalpos_a[i,2]
|
|
if alt:
|
|
res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos)
|
|
else:
|
|
res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos)
|
|
results[i,0] = res.x
|
|
results[i,1] = res.y
|
|
results[i,2] = res.z
|
|
return results.reshape(evalpos.shape)
|
|
|
|
|
|
cdef class ScatteringMatrix:
|
|
'''
|
|
Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix.
|
|
'''
|
|
cdef _ScatteringSystemAtOmega ssw # Here we keep the reference to the parent scattering system
|
|
cdef _ScatteringSystemAtOmegaK sswk
|
|
cdef qpms_ss_LU lu
|
|
|
|
def __cinit__(self, _ScatteringSystemAtOmega ssw, sswk=None, iri=None):
|
|
ssw.check()
|
|
self.ssw = ssw
|
|
if sswk is None:
|
|
ssw.ensure_finite()
|
|
# TODO? pre-allocate the matrix with numpy to make it transparent?
|
|
if iri is None:
|
|
self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU(
|
|
NULL, NULL, ssw.rawpointer())
|
|
else:
|
|
self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
|
|
NULL, NULL, ssw.rawpointer(), iri)
|
|
else:
|
|
# TODO check sswk validity
|
|
self.sswk = sswk
|
|
self.lu = qpms_scatsyswk_build_modeproblem_matrix_full_LU(NULL, NULL, self.sswk.rawpointer())
|
|
|
|
def __dealloc__(self):
|
|
qpms_ss_LU_free(self.lu)
|
|
|
|
property iri:
|
|
def __get__(self):
|
|
return None if self.lu.full else self.lu.iri
|
|
|
|
def __call__(self, a_inc):
|
|
cdef size_t vlen
|
|
cdef qpms_iri_t iri = -1;
|
|
if self.lu.full:
|
|
vlen = self.lu.ssw[0].ss[0].fecv_size
|
|
if len(a_inc) != vlen:
|
|
raise ValueError("Length of a full coefficient vector has to be %d, not %d"
|
|
% (vlen, len(a_inc)))
|
|
else:
|
|
iri = self.lu.iri
|
|
vlen = self.lu.ssw[0].ss[0].saecv_sizes[iri]
|
|
if len(a_inc) != vlen:
|
|
raise ValueError("Length of a %d. irrep packed coefficient vector has to be %d, not %d"
|
|
% (iri, vlen, len(a_inc)))
|
|
a_inc = np.array(a_inc, dtype=complex, copy=False, order='C')
|
|
cdef const cdouble[::1] a_view = a_inc;
|
|
cdef np.ndarray f = np.empty((vlen,), dtype=complex, order='C')
|
|
cdef cdouble[::1] f_view = f
|
|
qpms_scatsys_scatter_solve(&f_view[0], &a_view[0], self.lu)
|
|
return f
|
|
|
|
def pitau(double theta, qpms_l_t lMax, double csphase = -1):
|
|
if(abs(csphase) != 1):
|
|
raise ValueError("csphase must be 1 or -1, is %g" % csphase)
|
|
cdef size_t nelem = qpms_lMax2nelem(lMax)
|
|
cdef np.ndarray[np.float_t, ndim=1] lega = np.empty((nelem,), dtype=float)
|
|
cdef np.ndarray[np.float_t, ndim=1] pia = np.empty((nelem,), dtype=float)
|
|
cdef np.ndarray[np.float_t, ndim=1] taua = np.empty((nelem,), dtype=float)
|
|
cdef double[::1] leg = lega
|
|
cdef double[::1] pi = pia
|
|
cdef double[::1] tau = taua
|
|
qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase)
|
|
return (lega, pia, taua)
|
|
|
|
def linton_gamma(cdouble x):
|
|
return clilgamma(x)
|
|
|
|
def linton_gamma_real(double x):
|
|
return lilgamma(x)
|
|
|
|
def gamma_inc(double a, cdouble x, int m = 0):
|
|
cdef qpms_csf_result res
|
|
with pgsl_ignore_error(15): #15 is underflow
|
|
complex_gamma_inc_e(a, x, m, &res)
|
|
return (res.val, res.err)
|
|
|
|
def gamma_inc_series(double a, cdouble x):
|
|
cdef qpms_csf_result res
|
|
with pgsl_ignore_error(15): #15 is underflow
|
|
cx_gamma_inc_series_e(a, x, &res)
|
|
return (res.val, res.err)
|
|
|
|
def gamma_inc_CF(double a, cdouble x):
|
|
cdef qpms_csf_result res
|
|
with pgsl_ignore_error(15): #15 is underflow
|
|
cx_gamma_inc_CF_e(a, x, &res)
|
|
return (res.val, res.err)
|
|
|
|
def lll_reduce(basis, double delta=0.75):
|
|
"""
|
|
Lattice basis reduction with the Lenstra-Lenstra-Lovász algorithm.
|
|
|
|
basis is array_like with dimensions (n, d), where
|
|
n is the size of the basis (dimensionality of the lattice)
|
|
and d is the dimensionality of the space into which the lattice
|
|
is embedded.
|
|
"""
|
|
basis = np.array(basis, copy=True, order='C', dtype=np.double)
|
|
if len(basis.shape) != 2:
|
|
raise ValueError("Expected two-dimensional array (got %d-dimensional)"
|
|
% len(basis.shape))
|
|
cdef size_t n, d
|
|
n, d = basis.shape
|
|
if n > d:
|
|
raise ValueError("Real space dimensionality (%d) cannot be smaller than"
|
|
"the dimensionality of the lattice (%d) embedded into it."
|
|
% (d, n))
|
|
cdef double [:,:] basis_view = basis
|
|
if 0 != qpms_reduce_lattice_basis(&basis_view[0,0], n, d, delta):
|
|
raise RuntimeError("Something weird happened")
|
|
return basis
|
|
|
|
|
|
cdef PGen get_PGen_direct(direct_basis, bint include_origin=False, double layers=30):
|
|
dba = np.array(direct_basis)
|
|
if not (dba.shape == (2,2)):
|
|
raise NotImplementedError
|
|
cdef cart2_t b1, b2
|
|
b1.x = dba[0,0]
|
|
b1.y = dba[0,1]
|
|
b2.x = dba[1,0]
|
|
b2.y = dba[0,1]
|
|
cdef double maxR = layers*max(cart2norm(b1), cart2norm(b2))
|
|
return PGen_xyWeb_new(b1, b2, BASIS_RTOL, CART2_ZERO, 0, include_origin, maxR, False)
|
|
|
|
cdef double get_unitcell_volume(direct_basis):
|
|
dba = np.array(direct_basis)
|
|
if not (dba.shape == (2,2)):
|
|
raise NotImplementedError
|
|
cdef cart2_t b1, b2
|
|
b1.x = dba[0,0]
|
|
b1.y = dba[0,1]
|
|
b2.x = dba[1,0]
|
|
b2.y = dba[0,1]
|
|
return l2d_unitcell_area(b1, b2)
|
|
|
|
cdef PGen get_PGen_reciprocal2pi(direct_basis, double layers = 30):
|
|
dba = np.array(direct_basis)
|
|
if not (dba.shape == (2,2)):
|
|
raise NotImplementedError
|
|
cdef cart2_t b1, b2, rb1, rb2
|
|
b1.x = dba[0,0]
|
|
b1.y = dba[0,1]
|
|
b2.x = dba[1,0]
|
|
b2.y = dba[0,1]
|
|
if(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2) != 0):
|
|
raise RuntimeError
|
|
cdef double maxK = layers*max(cart2norm(rb1), cart2norm(rb2))
|
|
return PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, CART2_ZERO,
|
|
0, True, maxK, False)
|
|
|
|
cdef class Ewald3Calculator:
|
|
'''Wrapper class over qpms_ewald3_constants_t.
|
|
|
|
Mainly for testing low-level scalar Ewald summation functionality.'''
|
|
cdef qpms_ewald3_constants_t *c
|
|
|
|
def __cinit__(self, qpms_l_t lMax, int csphase = -1):
|
|
if (csphase != -1 and csphase != 1):
|
|
raise ValueError("csphase must be +1 or -1, not %d" % csphase)
|
|
self.c = qpms_ewald3_constants_init(lMax, csphase)
|
|
|
|
def __dealloc__(self):
|
|
qpms_ewald3_constants_free(self.c)
|
|
|
|
def sigma0(self, double eta, cdouble wavenumber, do_err = False):
|
|
cdef int retval
|
|
cdef double err
|
|
cdef cdouble result
|
|
retval = ewald3_sigma0(&result, &err, self.c, eta, wavenumber)
|
|
if retval:
|
|
raise RuntimeError("ewald3_sigma0 returned non-zero value (%d)" % retval)
|
|
if do_err:
|
|
return (result, err)
|
|
else:
|
|
return result
|
|
|
|
def sigma_short(self, double eta, cdouble wavenumber, direct_basis, wavevector, particle_shift, do_err=False):
|
|
# FIXME now only 2d XY lattice in 3D is implemented here, we don't even do proper dimensionality checks.
|
|
cdef cart3_t beta, pshift
|
|
beta.x = wavevector[0]
|
|
beta.y = wavevector[1]
|
|
beta.z = 0
|
|
pshift.x = particle_shift[0]
|
|
pshift.y = particle_shift[1]
|
|
pshift.z = 0
|
|
cdef qpms_l_t n = self.c[0].nelem_sc
|
|
cdef np.ndarray[complex, ndim=1] result = np.empty((n,), dtype=complex)
|
|
cdef cdouble[::1] result_v = result
|
|
cdef np.ndarray[double, ndim=1] err
|
|
cdef double[::1] err_v
|
|
if do_err:
|
|
err = np.empty((n,), dtype=np.double)
|
|
err_v = err
|
|
cdef bint include_origin = not (particle_shift[0] == 0 and particle_shift[1] == 0)
|
|
cdef PGen rgen = get_PGen_direct(direct_basis, include_origin)
|
|
cdef int retval = ewald3_sigma_short(&result_v[0], &err_v[0] if do_err else NULL,
|
|
self.c, eta, wavenumber, LAT_2D_IN_3D_XYONLY, &rgen, False, beta, pshift)
|
|
if rgen.stateData: PGen_destroy(&rgen)
|
|
if retval: raise RuntimeError("ewald3_sigma_short returned %d" % retval)
|
|
if do_err:
|
|
return (result, err)
|
|
else:
|
|
return result
|
|
|
|
def sigma_long(self, double eta, cdouble wavenumber, direct_basis, wavevector, particle_shift, do_err=False):
|
|
# FIXME now only 2d XY lattice in 3D is implemented here, we don't even do proper dimensionality checks.
|
|
cdef cart3_t beta, pshift
|
|
beta.x = wavevector[0]
|
|
beta.y = wavevector[1]
|
|
beta.z = 0
|
|
pshift.x = particle_shift[0]
|
|
pshift.y = particle_shift[1]
|
|
pshift.z = 0
|
|
cdef qpms_l_t n = self.c[0].nelem_sc
|
|
cdef np.ndarray[complex, ndim=1] result = np.empty((n,), dtype=complex)
|
|
cdef cdouble[::1] result_v = result
|
|
cdef np.ndarray[double, ndim=1] err
|
|
cdef double[::1] err_v
|
|
if do_err:
|
|
err = np.empty((n,), dtype=np.double)
|
|
err_v = err
|
|
cdef PGen kgen = get_PGen_reciprocal2pi(direct_basis)
|
|
cdef double unitcell_volume = get_unitcell_volume(direct_basis)
|
|
cdef int retval = ewald3_sigma_long(&result_v[0], &err_v[0] if do_err else NULL,
|
|
self.c, eta, wavenumber, unitcell_volume, LAT_2D_IN_3D_XYONLY, &kgen, False, beta, pshift)
|
|
|
|
if kgen.stateData: PGen_destroy(&kgen)
|
|
if retval: raise RuntimeError("ewald3_sigma_long returned %d" % retval)
|
|
if do_err:
|
|
return (result, err)
|
|
else:
|
|
return result
|
|
|
|
|
|
|
|
|