Compare commits

...

55 Commits

Author SHA1 Message Date
Marek Nečada 24c3ff386c Fix obvious errors caught by compiler 2022-10-22 17:28:40 +03:00
Marek Nečada 3a2fc6685c An ugly hack to convince newer version of numpy to save list of arrays 2022-06-24 08:42:56 +03:00
Marek Nečada e5270b54fb finiterect-constant-driving: sort slice labels for deterministic output 2022-06-24 02:52:39 +03:00
Marek Nečada ecac798529 Handle "empty" irreps in finiterectlat-scatter.py (not tested!) 2022-06-24 01:13:44 +03:00
Marek Nečada bb15a2b035 Handle "empty" irreps in finiterectlat-constant-driving.py 2022-06-23 18:35:33 +03:00
Marek Nečada 6b89feed70 Named constant for failed vswf index remap 2022-06-23 11:38:23 +03:00
Marek Nečada 1659588196 Fix qpms_vswf_set_reindex().
Stupid typo with possibly serious consequences...
2022-06-23 11:30:07 +03:00
Marek Nečada 0313d3b4ad vswf set specification support for "single particle type" scripts 2022-06-23 09:12:27 +03:00
Marek Nečada cfd25c6218 Constant T-matrix example add symmetries
continuous-integration/drone/push Build is failing Details
This should work, but it douesn't. This is obviously a bug.
2022-06-22 12:36:12 +03:00
Marek Nečada 8df02429d5 Initial example with manually specified constant T-matrix 2022-06-22 12:36:12 +03:00
Marek Nečada 5086e42353 argproc don't convert --vswf-set arguments directly to BaseSpec
BaseSpec is currently not picklable, and we want to be able to save
the command line arguments used.

Hence we first save the ilist as tuple instead.
2022-06-22 12:36:12 +03:00
Marek Nečada 74d72100d4 Fix typos in argproc.py, metavars to "labeled" parameters 2022-06-22 12:36:12 +03:00
Marek Nečada ccfb13a15c argproc.py WIP constant t-matrix, mutual exclusivity checks 2022-06-22 12:36:12 +03:00
Marek Nečada 68d894ec96 WIP argproc.py constant T-matrix 2022-06-22 12:36:12 +03:00
Marek Nečada 8a2d4a2969 Forgotten source catch_aux.C
continuous-integration/drone/push Build is failing Details
2022-06-22 12:33:48 +03:00
Marek Nečada d062ec4da9 Forgotten header lattice_types.h 2022-06-22 12:33:48 +03:00
Marek Nečada cd1406d23f argproc.py: enable manual specification of BaseSpec 2022-06-22 12:33:48 +03:00
Marek Nečada 4c0453a865 Nicer BaseSpec.__repr__() for BaseSpec(lMax=...) 2022-06-22 12:33:48 +03:00
Marek Nečada f0d0e48876 BaseSpec.__repr__() 2022-06-22 12:33:48 +03:00
Marek Nečada 6ddcc4cfcf BaseSpec constructor docs 2022-06-22 12:33:48 +03:00
Marek Nečada cf4b9cb52b Comments, doxygen
Remove old (not valid anymore) comments from vswf.h

Doxygenise translations_dbg.h
2022-06-22 12:33:48 +03:00
Marek Nečada ffb4bab365 Dummy (zero) T-matrix generator 2022-06-22 12:33:48 +03:00
Marek Nečada 487bc25aaf VSWF fill test, fix memory issue 2022-06-22 12:33:48 +03:00
Marek Nečada 87079dd09b Missing include 2022-06-22 12:33:48 +03:00
Marek Nečada 5dba7acd40 Indexing.h C++ compatibility; flag C++ incompatible headers. 2022-06-22 12:33:48 +03:00
Marek Nečada d43ca9026b Restructure headers to ensure ewald.h C++ compatibility 2022-06-22 12:33:48 +03:00
Marek Nečada 31a1094d08 vectors.h explicit type conversions for C++ compatibility 2022-06-22 12:33:48 +03:00
Marek Nečada 0aab207ab0 extern "C" in headers 2022-06-22 12:33:48 +03:00
Marek Nečada 45dba25b3c Introduce testing witch catch2 (no runtime tests yet) 2022-06-22 12:33:48 +03:00
Marek Nečada 5af1fd7734 Replace "complex" with "_Complex" in headers
Needed for compatibility with C++ (in which "complex" is
a template).
2022-06-22 12:33:48 +03:00
Marek Nečada d640176d62 Expressing VSWF in terms of SSWF 2022-06-22 12:33:48 +03:00
Marek Nečada bac3c3276d Indexing documentation, pedantic use of qpms_y_sc_t 2022-06-22 12:33:48 +03:00
Marek Nečada 7e318fc6cd Remove duplicit function declaration.
qpms_scatsysw_scattered_field_basis()
2022-06-22 12:33:48 +03:00
Marek Nečada 9e9f70cf65 analysis: fix symmetry group info 2022-06-22 12:33:48 +03:00
Marek Nečada 030a96d56f Analysis script also for finiterectlat constant driving results 2022-06-22 12:33:48 +03:00
Marek Nečada 2c12d28c70 Expose more testing functions via cython api. 2022-06-22 12:33:48 +03:00
Marek Nečada fed391b20c analysis.py: fix duplicities in position mapping due to rounding errors 2022-06-22 12:33:48 +03:00
Marek Nečada 679ee72bee Helper class to analyse results of finiterectlat-modes.py 2022-06-22 12:33:48 +03:00
Marek Nečada dc0d030ef0 WIP fixing some example scripts 2022-06-22 12:33:48 +03:00
Marek Nečada 798d6e7943 Move old cython translation_calculator wrapper 2022-06-22 12:33:48 +03:00
Marek Nečada 7599c22d06 Attempt to fix the Ewald sum SSWF phase bug 2022-06-22 12:33:48 +03:00
Marek Nečada 788a990a29 Revert scatsystem cython parallelisation.
It didn't work correctly on my laptop.
This must be reviewed later, as the waiting is annoying.
2022-06-22 12:33:48 +03:00
Marek Nečada 1e5418e7d0 Dirty fix of the otherwise non-literal arrays. 2022-06-22 12:33:48 +03:00
Marek Nečada a66dd24f89 Cython evaluation of "normal" and Ewald-summed SSWF 2022-06-22 12:33:48 +03:00
Marek Nečada 94e218c459 ScatteringSystem, ... definitions to pxd file 2022-06-22 12:33:48 +03:00
Marek Nečada 3a8e5b99cb SSWF evaluation functions for debugging Ewald sum. 2022-06-22 12:33:48 +03:00
Marek Nečada 6ccd785164 Additional phase factors in Ewald sum for debugging 2022-06-22 12:33:48 +03:00
Marek Nečada e486e04e2f WIP documentation on phase/normalisation implementation. 2022-06-22 12:33:48 +03:00
Marek Nečada fd588b7c02 Ewald summation consistency test for lattice basis fields.
In-plane test currently passes, out-of-plane does not.
2022-06-22 12:33:48 +03:00
Marek Nečada 2cef16834a finiterect-constant-driving more metadata output 2022-06-22 12:33:48 +03:00
Marek Nečada c0ee750f5d Remove bugs of removed code from BUGS.rst.
Some of the other bugs might have been already fixed
as well, this has to be reviewed, and also doxygenised.
2022-06-22 12:33:48 +03:00
Marek Nečada 8abe7e3357 Enable separation of long- and short-range of Ewald sum in python 2022-06-22 12:33:48 +03:00
Marek Nečada 5fdbb362a3 Fix prange race conditions. 2022-06-22 12:33:48 +03:00
Marek Nečada 478b010a15 "Basis fields" for finite systems.
There seem to be race conditions in the prange'd cython parts.
2022-06-22 12:33:48 +03:00
Marek Nečada 9e350cdac6 "Basis fields" for finite systems 2022-06-22 12:33:48 +03:00
75 changed files with 2384 additions and 638 deletions

View File

@ -23,28 +23,6 @@ Plane wave decompositions gives wrong value on the longitudinal part.
The implementation of the L coefficients OR the longitudinal waves
is thus probably wrong.
Scattering result asymmetry
---------------------------
The lattice scattering code (such as finitesqlatzsym-scattery.py) produces
asymmetric results where one should not get them, due to the system symmetry.
It seems that the asymmetry appears mostly in the y-direction (i.e.
for example the scattering/absorption cross section at k = (kx, ky, kz)
is not exactly the same as k = (kx, -ky, kz).
What has been checked (hopefully):
- The flip operators for electric waves
- Some weird kind of rounding or other numerical error depending on
the position ordering of the matrix (randomized indices give
the same asymmetry).
What has not been checked (so well):
- The x,y,z-flip operators for magnetic waves (i.e. are they really
supposet to bejust the
same as the operators for electric waves, just with opposite sign?)
- zplane_pq_y
- the translation operators
Singular value asymmetry
------------------------
@ -60,6 +38,3 @@ Moreover, the non-normalized legendre functions that are used in translations.c
are likely to overflow for even smaller values of l.
Therefore TODO: modify translations.c in a way that uses normalised legendre functions everywhere.
Memory management
-----------------
When compiled with optimizations, I get stack smashing error in qpms_trans_calculator_get_2DFT_longrange_AB_arrays().

View File

@ -48,6 +48,7 @@ add_subdirectory(faddeeva)
add_subdirectory (qpms)
add_subdirectory (tests/catch EXCLUDE_FROM_ALL)
enable_testing()
add_subdirectory (tests EXCLUDE_FROM_ALL)

View File

@ -9,10 +9,10 @@ ${MISCDIR}/lat2d_modes.py \
-b s${A1X_nm}e-9 s${A1Y_nm}e-9 \
-b s${A2X_nm}e-9 s${A2Y_nm}e-9 \
-p s${P1X_nm}e-9 s${P1Y_nm}e-9 \
-p s${P2X_nm}e-9 s${P2Y_nm}e-9 \
-L 3 -m $METAL -r ${RADIUS_nm}e-9 -H ${HEIGHT_nm}e-9 \
-k s${KPOINTX_nmi}e9 s${KPOINTY_nmi}e9 \
-k 0 0 \
-d -3 \
-t 0.01 \
-c 250 \
-P
#-k s${B1X_nmi}e9 s${B1Y_nmi}e9 \

View File

@ -0,0 +1,7 @@
An example with arectangular array of a very simplistic
model of scatterer with completely fabricated constant
T-matrix and manually specified set of VSWFs the scatterer
interacts with (i.e. instead of constant lMax cutoff).
This might be useful e.g. in symmetry-breaking thought
experiments.

View File

@ -0,0 +1,53 @@
#!/bin/bash
# Common parameters for a rectangular array
# N.B. We put those into bc, which does not understant exponential notation
export PX_nm=375
export PY_nm=375
export BG_REFINDEX=1.52
export BSPEC='[6,14]' # l = 1, m = +- 1 (i.e. xy) electric dipoles, nothing more
# Input parameters for T-matrix. It shall have the following form:
# '[[ T_mm, T_mp],
# [ T_mp, T_pp]]'
#
# The following lines are used to specify parameters for the most
# simplistic T-matrix that b having only diagonal elements, but
T_PHASE_MM=0.1
T_PHASE_PP=0.05
# "gain" factors; values larger than 1 are unphysical, 1 conserves energy
T_GAINFACTOR_MM=1
T_GAINFACTOR_PP=1
# Setup bc
echo 'scale=20;pi=3.14159265358979323846;' > bc_env
export BC_ENV_ARGS="bc_env"
# We have only one particle per unit cell here
export P1X_nm=0
export P1Y_nm=0
# Lattice vectors (for the general scripts)
export A1X_nm=${PX_nm}
export A1Y_nm=0
export A2X_nm=0
export A2Y_nm=${PY_nm}
# Reciprocal lattice vectors
export B1X_nmi=$(bc <<< '1/'${PX_nm})
export B1Y_nmi=0
export B2X_nmi=0
export B2Y_nmi=$(bc <<< '1/'${PY_nm})
# Constant T-matrix
T_MM="$(bc -l <<< -.5+.5*${T_GAINFACTOR_MM}*c\(${T_PHASE_MM}\))+$(bc -l <<< .5*${T_GAINFACTOR_MM}*s\(${T_PHASE_MM}\))j"
T_PP="$(bc -l <<< -.5+.5*${T_GAINFACTOR_PP}*c\(${T_PHASE_PP}\))+$(bc -l <<< .5*${T_GAINFACTOR_PP}*s\(${T_PHASE_PP}\))j"
export TMATRIX="\
[[${T_MM},0],\
[0,${T_PP}]]" # note no unescaped whitespaces allowed here

View File

@ -0,0 +1,17 @@
#!/bin/bash
SCRIPTDIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
MISCDIR=../../../misc
source ${SCRIPTDIR}/00_params.sh
# try several lMaxes
${MISCDIR}/lat2d_realfreqsvd.py \
-B $BG_REFINDEX \
-b s${A1X_nm}e-9 s${A1Y_nm}e-9 \
-b s${A2X_nm}e-9 s${A2Y_nm}e-9 \
-p s${P1X_nm}e-9 s${P1Y_nm}e-9 \
-w ${BSPEC} -T ${TMATRIX} \
-k 0 0 \
-F 2.001 0.001 2.250 \
-P

View File

@ -0,0 +1,18 @@
#!/bin/bash
SCRIPTDIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
MISCDIR=../../../misc
source ${SCRIPTDIR}/00_params.sh
# try several lMaxes
${MISCDIR}/lat2d_realfreqsvd.py \
-B $BG_REFINDEX \
-b s${A1X_nm}e-9 s${A1Y_nm}e-9 \
-b s${A2X_nm}e-9 s${A2Y_nm}e-9 \
-p s${P1X_nm}e-9 s${P1Y_nm}e-9 \
-w ${BSPEC} -T ${TMATRIX} \
-k 0 0 \
-F 2.001 0.001 2.250 \
-g C4 \
-P

View File

@ -4,7 +4,7 @@ import math
from qpms.argproc import ArgParser, make_dict_action, sslice, annotate_pdf_metadata
figscale=3
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'single_omega'])
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax_or_vswfset', 'single_omega'])
ap.add_argument("-k", '--wavevector', nargs=2, type=float, required=True, help='"Bloch" vector, modulating phase of the driving', metavar=('KX', 'KY'), default=(0., 0.))
# ap.add_argument("--kpi", action='store_true', help="Indicates that the k vector is given in natural units instead of SI, i.e. the arguments given by -k shall be automatically multiplied by pi / period (given by -p argument)")
ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
@ -36,8 +36,8 @@ px, py = a.period
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
defaultprefix = "cd_%s_p%gnmx%gnm_%dx%d_m%s_n%s_k_%g_%g_f%geV_L%d_micro-%s" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.wavevector[0], a.wavevector[1], a.eV, a.lMax, "SO3" if a.symmetry_adapted is None else a.symmetry_adapted)
defaultprefix = "cd_%s_p%gnmx%gnm_%dx%d_m%s_n%s_k_%g_%g_f%geV_%s_micro-%s" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.wavevector[0], a.wavevector[1], a.eV, ap.bspecstr, "SO3" if a.symmetry_adapted is None else a.symmetry_adapted)
logging.info("Default file prefix: %s" % defaultprefix)
import numpy as np
@ -54,6 +54,8 @@ eh = eV/hbar
# Check slice ranges and generate all corresponding combinations
slicepairs = []
slicelabels = set(a.xslice.keys()) | set(a.yslice.keys())
slicelabels = sorted(slicelabels, key= lambda l:
'.None' if l is None else str(l))
for label in slicelabels:
rowslices = a.xslice.get(label, None)
colslices = a.yslice.get(label, None)
@ -71,13 +73,17 @@ def realdipfieldlabels(yp):
if yp == 1: return 'y'
if yp == 2: return 'z'
raise ValueError
def realdipfields(vecgrid, yp):
def realdipfields(vecgrid, yp, bspec=BaseSpec(lMax=1)):
if yp == 0 or yp == 1:
im = np.where(bspec.ilist == 6)[0][0]
ip = np.where(bspec.ilist == 14)[0][0]
if yp == 1:
return vecgrid[...,0] + vecgrid[...,2]
return vecgrid[...,ip] + vecgrid[...,ip]
if yp == 0:
return -1j*(vecgrid[...,0] - vecgrid[...,2])
return -1j*(vecgrid[...,im] - vecgrid[...,ip])
if yp == 2:
return vecgrid[...,1]
i0 = np.where(bspec.ilist == 10)[0][0]
return vecgrid[...,i0]
raise ValueError
def float_nicestr(x, tol=1e-5):
@ -126,7 +132,7 @@ orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
omega = ap.omega
bspec = BaseSpec(lMax = a.lMax)
bspec = ap.bspec
medium = EpsMuGenerator(ap.background_epsmu)
particles= [Particle(orig_xy[i], ap.tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]
@ -199,6 +205,9 @@ scattered_ir = [None for iri in range(ss.nirreps)]
ir_contained = np.ones((nsp, nelem, ss.nirreps), dtype=bool)
for iri in range(ss.nirreps):
if ss.saecv_sizes[iri] == 0:
logging.info('irrep %d/%d has an empty VSWF set, skipping' % (iri, ss.nirreps))
continue
logging.info("processing irrep %d/%d" % (iri, ss.nirreps))
LU = None # to trigger garbage collection before the next call
translation_matrix = None
@ -251,10 +260,15 @@ if not math.isnan(a.ccd_distance):
logging.info("Far fields done")
outfile = defaultprefix + ".npz" if a.output is None else a.output
for sir in scattered_ir:# An ugly hack to convince newer versions of numpy that we indeed want to save a list of arrays
if sir is None:
break
else:
scattered_ir.append(None)
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()},
omega=omega, wavenumber=wavenumber, nelem=nelem, wavevector=np.array(a.wavevector), phases=phases,
positions = ss.positions[:,:2],
scattered_ir_packed = np.array(scattered_ir, dtype=np.object),
scattered_ir_packed = np.array(scattered_ir, dtype=object),
scattered_full = scattered_full,
ir_contained = ir_contained,
t=t, l=l, m=m,
@ -264,6 +278,7 @@ np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()},
#ccd_size = ccd_size if not math.isnan(a.ccd_distance) else None,
ccd_points = ccd_points if not math.isnan(a.ccd_distance) else None,
ccd_fields = ccd_fields if not math.isnan(a.ccd_distance) else None,
fullvec_poffsets = ss.fullvec_poffsets,
)
logging.info("Saved to %s" % outfile)
@ -317,15 +332,25 @@ if a.plot or (a.plot_out is not None):
vecgrid_ff = np.fft.fftshift(np.fft.fft2(vecgrid, axes=(0,1)),axes=(0,1))
lemax = np.amax(abs(vecgrid))
for yp in range(0,3):
if(np.amax(abs(realdipfields(vecgrid,yp))) > lemax*1e-5):
axes[y,yp*4].imshow(abs(realdipfields(vecgrid,yp)), vmin=0, interpolation='none')
axes[y,yp*4].text(0.5, 0.5, '%g' % np.amax(abs(realdipfields(vecgrid,yp))), horizontalalignment='center', verticalalignment='center', transform=axes[y,yp*4].transAxes)
axes[y,yp*4+1].imshow(np.angle(realdipfields(vecgrid,yp)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none')
axes[y,yp*4+2].imshow(abs(realdipfields(vecgrid_ff,yp)), vmin=0, interpolation='none')
axes[y,yp*4+3].imshow(np.angle(realdipfields(vecgrid_ff,yp)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none')
try:
if(np.amax(abs(realdipfields(vecgrid,yp,bspec))) > lemax*1e-5):
axes[y,yp*4].imshow(abs(realdipfields(vecgrid,yp,bspec)), vmin=0, interpolation='none')
axes[y,yp*4].text(0.5, 0.5, '%g' % np.amax(abs(realdipfields(vecgrid,yp,bspec))), horizontalalignment='center', verticalalignment='center', transform=axes[y,yp*4].transAxes)
axes[y,yp*4+1].imshow(np.angle(realdipfields(vecgrid,yp,bspec)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none')
axes[y,yp*4+2].imshow(abs(realdipfields(vecgrid_ff,yp,bspec)), vmin=0, interpolation='none')
axes[y,yp*4+3].imshow(np.angle(realdipfields(vecgrid_ff,yp,bspec)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none')
else:
for c in range(0,4):
axes[y,yp*4+c].tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)
except (KeyError, IndexError) as e:
logging.info("skipping dipole plot #%d for y=%d, spi=%d (dipole probably not included in the VSWF set) %s" %
(yp, y, spi, e))
for c in range(0,4):
axx= axes[y,yp*4+c]
axx.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)
axx.text(0.5, 0.5, 'skipped', horizontalalignment='center',verticalalignment='center', transform=axx.transAxes)
if not math.isnan(a.ccd_distance):
fxye=(-ccd_size/2, ccd_size/2, -ccd_size/2, ccd_size/2)
e2vmax = np.amax(np.linalg.norm(ccd_fields[spi,y], axis=-1)**2)

View File

@ -4,7 +4,7 @@ import math
from qpms.argproc import ArgParser, annotate_pdf_metadata
ap = ArgParser(['rectlattice2d_finite', 'background_analytical', 'single_particle', 'single_lMax', ])
ap = ArgParser(['rectlattice2d_finite', 'background_analytical', 'single_particle', 'single_lMax_or_vswfset', ])
ap.add_argument("-t", "--rank-tolerance", type=float, default=1e11)
ap.add_argument("-c", "--min-candidates", type=int, default=1, help='always try at least this many eigenvalue candidates, even if their SVs in the rank tests are lower than rank_tolerance')
@ -35,8 +35,8 @@ thegroup = 'D4h' if px == py and Nx == Ny and not a.D2 else 'D2h'
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_B%s_L%d_c(%s±%g±%gj)eV_cn%d_%s" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.lMax,
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_B%s_%s_c(%s±%g±%gj)eV_cn%d_%s" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), ap.bspecstr,
str(a.centre), a.ar, a.ai, a.N,
thegroup,
)
@ -70,7 +70,7 @@ orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
bspec = BaseSpec(lMax = a.lMax)
bspec = ap.bspec
medium = EpsMuGenerator(ap.background_epsmu)
particles= [Particle(orig_xy[i], ap.tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]

View File

@ -5,7 +5,7 @@ pi = math.pi
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave'])
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax_or_vswfset', 'omega_seq_real_ng', 'planewave'])
ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)")
ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path")
@ -36,8 +36,8 @@ px, py = a.period
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_bg%s%gπ_θ(%g_%g)π_ψ%gπ_χ%gπ_f%s_L%d" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.phi/pi, np.amin(a.theta)/pi, np.amax(a.theta)/pi, a.psi/pi, a.chi/pi, ap.omega_descr, a.lMax, )
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_bg%s%gπ_θ(%g_%g)π_ψ%gπ_χ%gπ_f%s_%s" % (
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.phi/pi, np.amin(a.theta)/pi, np.amax(a.theta)/pi, a.psi/pi, a.chi/pi, ap.omega_descr, ap.bspecstr, )
logging.info("Default file prefix: %s" % defaultprefix)
@ -47,7 +47,7 @@ orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
bspec = BaseSpec(lMax = a.lMax)
bspec = ap.bspec
particles= [Particle(orig_xy[i], ap.tmgen, bspec=bspec) for i in np.ndindex(orig_xy.shape[:-1])]
sym = FinitePointGroup(point_group_info['D2h'])
@ -70,8 +70,8 @@ ndir = a.theta.shape[0]
k_cart_arr = np.empty((nfreq, ndir, 3), dtype=float)
wavenumbers = np.empty((nfreq,), dtype=float)
σ_ext_arr_ir = np.empty((nfreq, ndir, ss.nirreps), dtype=float)
σ_scat_arr_ir = np.empty((nfreq, ndir, ss.nirreps), dtype=float)
σ_ext_arr_ir = np.zeros((nfreq, ndir, ss.nirreps), dtype=float)
σ_scat_arr_ir = np.zeros((nfreq, ndir, ss.nirreps), dtype=float)
outfile_tmp = defaultprefix + ".tmp" if a.output is None else a.output + ".tmp"
@ -90,6 +90,9 @@ for i, omega in enumerate(ap.allomegas):
k_cart_arr[i] = sph2cart(k_sph_list)
for iri in range(ss.nirreps):
if ss.saecv_sizes[iri] == 0:
logging.info('irrep %d/%d has an empty VSWF set, skipping' % (iri, ss.nirreps))
continue
logging.info("processing irrep %d/%d" % (iri, ss.nirreps))
LU = None # to trigger garbage collection before the next call
translation_matrix = None

View File

@ -3,7 +3,7 @@
import math
from qpms.argproc import ArgParser, annotate_pdf_metadata
ap = ArgParser(['rectlattice2d', 'const_real_background', 'single_particle', 'single_lMax']) # const_real_background needed for calculation of the diffracted orders
ap = ArgParser(['rectlattice2d', 'const_real_background', 'single_particle', 'single_lMax_or_vswfset']) # const_real_background needed for calculation of the diffracted orders
ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
ap.add_argument("--kpi", action='store_true', help="Indicates that the k vector is given in natural units instead of SI, i.e. the arguments given by -k shall be automatically multiplied by pi / period (given by -p argument)")
ap.add_argument("--rank-tol", type=float, required=False)
@ -28,8 +28,8 @@ if a.kpi:
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
defaultprefix = "%s_p%gnmx%gnm_m%s_n%g_b%+d_k(%g_%g)um-1_L%d_cn%d" % (
particlestr, px*1e9, py*1e9, str(a.material), a.refractive_index, a.band_index, a.k[0]*1e-6, a.k[1]*1e-6, a.lMax, a.N)
defaultprefix = "%s_p%gnmx%gnm_m%s_n%g_b%+d_k(%g_%g)um-1_%s_cn%d" % (
particlestr, px*1e9, py*1e9, str(a.material), a.refractive_index, a.band_index, a.k[0]*1e-6, a.k[1]*1e-6, ap.bspecstr, a.N)
import logging
logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
@ -91,7 +91,7 @@ freqradius = .5 * (top - bottom) * a.interval_factor
centfreq = bottom + freqradius if a.band_index > 0 else top - freqradius
bspec = BaseSpec(lMax = a.lMax)
bspec = ap.bspec
pp = Particle(orig_xy[0][0], t = ap.tmgen, bspec=bspec)
ss, ssw = ScatteringSystem.create([pp], ap.background_emg, centfreq, latticebasis = ap.direct_basis)

33
notes/VSWF_from_SSWF.md Normal file
View File

@ -0,0 +1,33 @@
VSWF expansions in terms of SSWF
================================
From
\cite necada_multiple-scattering_2021, eq. (2.19)
\f[
\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
\sum_{\tau'l'm'} \tropSr{\kappa(\vect r_2 - \vect r_1)}_{\tau l m;\tau'l'm} \wfkcreg_{\tau'l'm'}(\vect r -\vect r_2),
\qquad |\vect r -\vect r_2| < |\vect r_1 - \vect r_2|,
\f]
setting \f$ \vect r = \vect r_2\f$ and considering that
\f$ \wfkcreg_{\tau'l'm'}(\vect 0) \ne \vect 0 \f$ only for electric dipole waves (\f$ \tau = \mathrm{E}, l=1 \f$),
we have
\f[
\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
\sum_{m'} \tropSr{\kappa(\vect r - \vect r_1)}_{\tau l m;\mathrm{E}1m} \wfkcreg_{\mathrm{E}1m'}(\vect 0),
\qquad \vect r \ne \vect r_1 .
\f]
Combining this with
\cite necada_multiple-scattering_2021, eq. (2.25)
\f[
\tropSr{\vect d}_{\tau l m; \tau' l' m'} = \sum_{\lambda =|l-l'|+|\tau-\tau'|}^{l+l'}
C^{\lambda}_{\tau l m;\tau' l'm'} \underbrace{ \spharm{\lambda}{m-m'}(\uvec d) h_\lambda^{(1)}(d)}_{\sswfout_\lambda^{m-m'}(\vect d)},
\f]
we get
\f[
\wfkcout_{\tau lm}(\vect d) = \sum_{m'=-1}^1 \wfkcreg_{\mathrm{E}1m'}(\vect 0)
\sum_{\lambda=l-1+|\tau-\tau'|}^{l+1}
C^\lambda_{\tau l m;\mathrm{E}1m'} \sswfout_\lambda^{m-m'}(\vect d).
\f]
Note that the VSWF components in this expression are given in global "cartesian" basis,
*not* the local orthonormal basis derived from spherical coordinates.
(This is mostly desirable, because in lattices we need to work with flat coordinates anyway.)

View File

@ -1,6 +1,10 @@
VSWF conventions {#vswf_conventions}
====================================
*This page provides reference about the VSWF conventions used in the literature.
For VSWF convention specification in QPMS API, see
[SWF Conventions in QPMS](@ref swf_conventions_qpms).*
In general, the (transversal) VSWFs can be defined using (some) vector spherical harmonics
as follows: \f[
\wfm\pr{k\vect r}_{lm} = \sphbes_l(kr) \vshrot_{lm} (\uvec r),\\
@ -50,7 +54,7 @@ where the connection to negative orders is
\dlmfFer{\nu}{m}(x) = (-1)^m \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfFer{\nu}{m}(x),\\
%\dlmfLeg{\nu}{m}(x) = \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfLeg{\nu}{m}(x).\\
\f]
Note that there are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
Note that they are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
different meaning / conventions (Ferrers functions being defined for \f$ \abs{x} \le 1 \f$, whereas
Legendre for \f$ \abs{x} \ge 1 \f$. We will not use the DLMF "Legendre" functions here.

114
notes/conventions_qpms.md Normal file
View File

@ -0,0 +1,114 @@
SWF conventions in QPMS {#swf_conventions_qpms}
=================================================
*This page describes how (V)SWF conventions are specified
internally and in QPMS API. For a general overview of VSWF
conventions in the literature, see [VSWF Conventions](@ref vswf_conventions).*
Convention enumerator
---------------------
Most of the meaningful phase and normalisation conventions for spherical waves
can be specified by the enum type @ref qpms_normalisation_t.
The type can be also used to specify conventions that are currently not fully
supported in QPMS (such as those based on real spherical harmonics).
As an enum type, it does not cover all the conventions possibly imaginable,
but it does cover the most meaningful ones and most of those that can be found
in the literature.
(Most notably, it does not cover the “anti-normalisation”
that does appear in certain expressions in some literature where the spherical
harmonics contain unnormalised Legendre functions, so that the basis set of
of spherical harmonics has different norms for different signs of *m*
for the same *l*. This is a bad idea overall and an absolutely atrocious
approach for numerics. Do not use that.)
VSWF evaluation
---------------
Evaluation of VSWFs using qpms_vswf_fill(), qpms_eval_vswf(),
qpms_uvswf_fill() and other functions from vswf.h are evaluated as follows.
These fuctions take a @ref qpms_normalisation_t as an argument.
The Ferrers-Legendre functions and the π, τ functions are evaluated
by qpms_pitau_get(), which internally uses gsl_sf_legendre_deriv_array_e().
Note only the information about the Condon-Shortley
phase is passed to qpms_pitau_get() the result of this function
uses always the GSL_SF_LEGENDRE_SPHARM normalisation,
and possible normalisation and other phase factors are evaluated afterwards
using the inline functions
qpms_normalisation_factor_L_noCS(),
qpms_normalisation_factor_M_noCS(),
qpms_normalisation_factor_N_noCS().
Evaluation of vector spherical harmonics only with qpms_vecspharm_fill()
works similarly but TODO.
TODO reference to pi, tau.
VSWF translation operator evaluation
------------------------------------
In practice, translation operators are calculated by first creating
an instance of the qpms_trans_calculator structure, which contains
a table of constant normalisation factors for a given phase/normalisation
convention (it is assumed that the phase/normalisation conventions do not
change with the translation), and then calling
qpms_trans_calculator_get_AB_arrays()
(or others).
The precomputed factor table in qpms_trans_calculator_t contains a CS phase
related factor (via qpms_trans_normfac()).
Function qpms_trans_calculator_get_AB_arrays_buf() then calculates
the unnormalised (GSL_SF_LEGENDRE_NONE)
associated Legendre functions always with CS phase -1;
and the A, B arrays are filled (via qpms_trans_calculator_get_AB_arrays_precalcbuf())
by individual calls of qpms_trans_calculator_get_A_precalcbuf()
and qpms_trans_calculator_B_precalcbuf(), which basically just multiply
and sum the precalculated constant factors with the radial (Bessel),
polar (Legendre) and azimuthal (exponential/trigonometric) functions.
**This means that the normalisation and phase convention is fully embedded
in the constant factor tables, and nothing is calculated during "runtime".**
The "higher-level" qpms_trans_calculator_get_trans_array() currently
just calls qpms_trans_calculator_get_AB_arrays() and then reorders
the elements (using qpms_trans_array_from_AB()), asserting that
the normalisation conventions remain the same.
There seems to be an inconsistency between
qpms_trans_calculator_get_B_buf() and
qpms_trans_calculator_get_A_buf() on one hand, and
qpms_trans_calculator_get_AB_buf_p() and
qpms_trans_calculator_get_AB_arrays_buf() on the other.
While the latter two functions use always -1 as the CS phase,
the former two take it from the normalisation enumerator.
**Although the former two are probably used nowhere in the production,
this needs to be fixed.**
Lattice sums
------------
### Scalar SWFs
### Translation operators
Function qpms_trans_calculator_get_AB_arrays_e32_e()
first compute the scalar lattice sums (using ewald3_sigma_short(),
ewald3_sigma_long() and ewald3_sigma0() calls).
These are then transformed into the VSWF translation operator
elements in a similar manner as in
qpms_trans_calculator_get_A_precalcbuf() and
qpms_trans_calculator_get_B_precalcbuf(), although there some optical
differences (CHECK!).
### VSWFs

View File

@ -31,8 +31,15 @@ MathJax.Hub.Config({
spharm: ["{{Y_{\\mathrm{#1}}}_{#2}^{#3}}", 3, ""], // Spherical harmonics
spharmR: ["{{Y_{\\mathrm{#1}}}_{\\mathrm{#1}{#2}{#3}}", 4, ""], // Spherical harmonics
csphase: "\\mathsf{C_{CS}}", // Condon-Shortley phase
tropSrr: ["{{S^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
tropRrr: ["{{R^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
tropS: "{\\mathcal{S}}", // Translation operator singular
tropR: "{\\mathcal{R}}", // Translation operator regular
tropSr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator singular
tropRr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator regular
tropSrr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
tropRrr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
sswfout: "{\\psi}", // outgoing SSWF
sswfreg: "{\\phi}", // regular SSWF
// Kristensson's VSWFs, complex version (2014 notes)
wfkc: "{\\vect{y}}", // any wave

View File

@ -24,6 +24,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
lll.c beyn.c trivialgroup.c version.c
translations_dbg.c scatsystem_dbg.c
)
use_c99()

View File

@ -17,7 +17,7 @@ except ImportError as ex:
from .cyquaternions import CQuat, IRot3
from .cybspec import VSWFNorm, BaseSpec, default_bspec
from .cytmatrices import CTMatrix, TMatrixInterpolator, TMatrixGenerator
from .cytranslations import trans_calculator
#from .cytranslations import trans_calculator
from .cymaterials import MaterialInterpolator, EpsMu, LorentzDrudeModel, lorentz_drude, EpsMuGenerator
from .cycommon import dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType, VSWFType
from .cywaves import vswf_single

150
qpms/analysis.py Normal file
View File

@ -0,0 +1,150 @@
'''
Auxiliary functions to analyse results obtained using scripts in
the ../misc directory.
This needs to be kept in sync with the scripts and qpms/argproc.py.
'''
import math
import numpy as np
import qpms
from qpms.qpms_c import Particle
from qpms.cytmatrices import CTMatrix, TMatrixGenerator
from qpms.cybspec import BaseSpec
from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
from qpms.symmetries import point_group_info
from qpms import FinitePointGroup, ScatteringSystem, eV, hbar
eh = eV / hbar
def cleanarray(a, atol=1e-8, copy=True, fullNaN=True):
a = np.array(a, copy=copy)
sieve = abs(a.real) < atol
a[sieve] = 1j * a[sieve].imag
sieve = abs(a.imag) < atol
a[sieve] = a[sieve].real
if fullNaN:
if np.all(a == 0):
a[...] = 0
return a
def nicerot(a, atol=1e-10, copy=True):
'''
Gives an array a "nice" phase.
'''
a = np.array(a, copy=copy)
i = np.argmax(abs(a))
a = a / a[i] * abs(a[i])
return a
def _meta_matspec2emg(matspec):
if isinstance(matspec, (float, complex)): # number is interpreted as refractive index
return EpsMuGenerator(EpsMu(matspec**2))
else:
return EpsMuGenerator(lorentz_drude[matspec])
class FiniteRectLatAnalysis:
'''Recreates the scattering system structure from finiterectlat-modes.py or
finiterectlat-constant-driving.py'''
def __init__(self, data):
meta = data['meta'][()]
scatter = 'symmetry_adapted' in meta.keys()# finiterectlat-constant-driving ; TODO unify
if scatter:
#thegroup = meta['symmetry_adapted'] # always D2h, this option does something else!
thegroup = 'D2h'
else: # finiterectlat-modes
thegroup = 'D2h' if meta['D2'] else 'D4h'
Nx, Ny = meta['size']
px, py = meta['period']
orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
bspec = BaseSpec(lMax = meta['lMax'])
nelem = len(bspec)
bg = _meta_matspec2emg(meta['background'])
fg = _meta_matspec2emg(meta['material'])
if meta.get('height', None) is None:
tmgen = TMatrixGenerator.sphere(bg, fg, meta['radius'])
else:
tmgen = TMatrixGenerator.cylinder(bg, fg, meta['radius'], meta['height'], lMax_extend=meta['lMax_extend'])
particles = [Particle(orig_xy[i], tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]
sym = FinitePointGroup(point_group_info[thegroup])
omega = data['omega'][()] if scatter else meta['centre'] * eh
self.ss, self.ssw = ScatteringSystem.create(particles, bg, omega, sym=sym)
ss1, ssw1 = ScatteringSystem.create([Particle([0,0,0], tmgen, bspec)], bg, omega, sym=sym)
self.ss1, self.ssw1 = ss1, ssw1
# per-particle projectors/transformation to symmetry-adapted bases
fvcs1 = np.empty((nelem, nelem), dtype=complex)
y = 0
iris1 = []
for iri1 in range(ss1.nirreps):
for j in range(ss1.saecv_sizes[iri1]):
pvc1 = np.zeros((ss1.saecv_sizes[iri1],), dtype=complex)
pvc1[j] = 1
fvcs1[y] = ss1.unpack_vector(pvc1, iri1)
fvcs1[y] = cleanarray(nicerot(fvcs1[y], copy=False), copy=False)
y += 1
iris1.append(iri1)
# Mapping between ss particles and grid positions
positions = self.ss.positions.astype(np.float32)
# round to get rid of duplicities due to rounding errors
positions = positions.round(7-int(math.log(np.amax(abs(positions)),10)))
xpositions = np.unique(positions[:,0])
ypositions = np.unique(positions[:,1])
assert(len(xpositions) == Nx)
assert(len(ypositions == Ny))
# particle positions as integer indices
posmap = np.empty((positions.shape[0],2), dtype=int)
#invposmap = np.empty((Nx, Ny), dtype=int)
for i, pos in enumerate(positions):
posmap[i,0] = np.searchsorted(xpositions, positions[i,0])
posmap[i,1] = np.searchsorted(ypositions, positions[i,1])
#invposmap[posmap[i,0], posmap[i, 1]] = i
self.meta = meta
self.npzfile = data
self.fvcs1 = fvcs1
self.iris1 = np.array(iris1)
self.bspec = bspec
self.nelem = nelem
self.Nx, self.Ny, self.px, self.py = Nx, Ny, px, py
self.posmap = posmap
self.fullvec_poffsets = nelem * np.arange(Nx*Ny)
self.group = thegroup
if scatter:
self.typ = 'scatter'
self.fvcs1_fromdata = data['fvcs1']
self.iris1_fromdata = data['iris1']
self.positions_fromdata = data['positions']
self.driving_symmetry_adapted = data['symmetry_adapted']
else: # modes
self.typ = 'modes'
self.eigval = data['eigval']
self.neig = len(self.eigval)
self.eigvec = data['eigvec']
self.residuals = data['residuals']
self.ranktest_SV = data['ranktest_SV']
self.iri = data['iri'][()]
def fullvec2grid(self, fullvec, swapxy=False):
arr = np.empty((self.Nx, self.Ny, self.nelem), dtype=complex)
for pi, offset in enumerate(self.fullvec_poffsets):
ix, iy = self.posmap[pi]
arr[ix, iy] = fullvec[offset:offset+self.nelem]
return np.swapaxes(arr, 0, 1) if swapxy else arr
def saecv2grid(self, saecv, swapxy=False):
fullvec = self.ss.unpack_vector(saecv, iri=self.iri)
return self.fullvec2grid(fullvec, swapxy)
def eigvec_asgrid(self, i, swapxy = False):
'''
Returns i-th eigenmode as a (Nx, Ny, nelem)-shaped grid of excitation
coefficients.
'''
if self.iri is not None:
return self.saecv2grid(self.eigvec[i], swapxy)
else:
return self.fullvec2grid(self.eigvec[i], swapxy)

View File

@ -5,6 +5,8 @@ Common snippets for argument processing in command line scripts.
import argparse
import sys
import warnings
import ast
from collections import namedtuple
def flatten(S):
if S == []:
@ -174,6 +176,33 @@ def sint(string):
else: raise exc
return res
def _resolve_exclusive_groups(namespace, poslabel, *groups):
"""
namespace: arguments namespace whose elements
are dicts that possibly
each element of groups shall be a sequence of strings,
each of the strings should correspond
if mutual exclusivity is respected, returns index
of the "selected" group.
"""
selected = None
example = None
for i, group in enumerate(groups):
for name in group:
ledict = getattr(namespace, name, None)
if ledict is not None: # TODO Maybe it shuld always be not none?
res = ledict.get(poslabel, None)
if res is not None:
if selected is None:
selected = i
example = name
elif selected != i:
raise ArgumentProcessingError("Mutually exclusive parameters (%s and %s) specified for label %s." %
(name, example, poslabel))
if selected is None and poslabel is not None: # try default params
return _resolve_exclusive_groups(namespace, None, *groups)
return selected
def material_spec(string):
"""Tries to parse a string as a material specification, i.e. a
real or complex number or one of the string in built-in Lorentz-Drude models.
@ -192,6 +221,32 @@ def material_spec(string):
raise argparse.ArgumentTypeError("Material specification must be a supported material name %s, or a number" % (str(lorentz_drude.keys()),)) from ve
return lemat
class hashable_complex_matrix(tuple):
"""
Converts string to a hashable equivalent of two-dimensional
numpy.ndarray(..., dtype=complex) (which by itself is not hashable)
"""
def __new__(self, matrix):
import numpy as np
if isinstance(matrix, str):
matrix = ast.literal_eval(matrix)
matrix = np.array(matrix, dtype=complex, copy=False)
return super().__new__(hashable_complex_matrix, (tuple(row) for row in matrix))
# auxiliary structures for t-matrix generator specifications
constant_tmatrix_spec = namedtuple("constant_tmatrix_spec", ("bspec", "matrix"))
cyl_sph_dimensions = namedtuple("cyl_sph_dimensions", ("radius", "height", "lMax_extend"))
tmgen_spec = namedtuple("tmgen_spec", ("bgspec", "fgspec", "dims"))
def string2bspec_tuple(string):
"""
Converts string representation of list to BaseSpec...
And then to tuple, because we want the arguments to be picklable :(
"""
from .cybspec import BaseSpec
bspec = BaseSpec(ast.literal_eval(string))
return tuple(i for i in bspec.ilist)
class ArgParser:
''' Common argument parsing engine for QPMS python CLI scripts. '''
@ -223,31 +278,56 @@ class ArgParser:
mpgrp.add_argument("-p", "--position", nargs='+', action=make_dict_action(argtype=sfloat, postaction='append',
first_is_key=False), help="Particle positions, cartesion coordinates (default particle properties)")
mpgrp.add_argument("+p", "++position", nargs='+', action=make_dict_action(argtype=sfloat, postaction='append',
first_is_key=True), help="Particle positions, cartesian coordinates (labeled)")
first_is_key=True), metavar=('LABEL', 'POSITION'), help="Particle positions, cartesian coordinates (labeled)")
mpgrp.add_argument("-L", "--lMax", nargs=1, default={},
action=make_dict_action(argtype=int, postaction='store', first_is_key=False,),
help="Cutoff multipole degree (default)")
mpgrp.add_argument("+L", "++lMax", nargs=2,
action=make_dict_action(argtype=int, postaction='store', first_is_key=True,),
help="Cutoff multipole degree (labeled)")
metavar=('LABEL', 'LMAX'),
help="Cutoff multipole degree (labeled)")
mpgrp.add_argument("-m", "--material", nargs=1, default={},
action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=False,),
help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (default)')
mpgrp.add_argument("+m", "++material", nargs=2,
action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=True,),
metavar=('LABEL', 'MATERIAL'),
help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (labeled)')
mpgrp.add_argument("-r", "--radius", nargs=1, default={},
action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
help='particle radius (sphere or cylinder; default)')
mpgrp.add_argument("+r", "++radius", nargs=2,
action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
metavar=('LABEL', 'RADIUS'),
help='particle radius (sphere or cylinder; labeled)')
mpgrp.add_argument("-H", "--height", nargs=1, default={},
action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
help='particle radius (cylinder; default)')
mpgrp.add_argument("+H", "++height", nargs=2,
action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
metavar=('LABEL', 'HEIGȞT'),
help='particle radius (cylinder; labeled)')
# Alternatively, add a constant T-matrix
mpgrp.add_argument("-w", "--vswf-set", nargs=1, default={},
action=make_dict_action(argtype=string2bspec_tuple, postaction='store', first_is_key=False),
help='Manual specification of VSWF set codes (format as a python list of integers); see docs on qpms_uvswfi_t for valid codes or simply use --lMax instead. Overrides --lMax.')
mpgrp.add_argument("+w", "++vswf-set", nargs=2, default={},
action=make_dict_action(argtype=string2bspec_tuple, postaction='store', first_is_key=True),
metavar=('LABEL', "VSWF_SET"),
help='Manual specification of VSWF set codes (format as a python list of integers); see docs on qpms_uvswfi_t for valid codes or simply use ++lMax instead. Overrides ++lMax and --lMax.')
mpgrp.add_argument("-T", "--constant-tmatrix", nargs=1, default={},
action=make_dict_action(argtype=hashable_complex_matrix, postaction='store', first_is_key=False),
help='constant T-matrix (elements must correspond to --vswf-set)')
mpgrp.add_argument("+T", "++constant-tmatrix", nargs=2, default={},
metavar=('LABEL', 'TMATRIX'),
action=make_dict_action(argtype=hashable_complex_matrix, postaction='store', first_is_key=True),
help='constant T-matrix (elements must correspond to ++vswf-set)')
def __add_single_lMax_or_vswfset_group(ap):
grp = ap.add_mutually_exclusive_group(required=True)
grp.add_argument("-w", "--vswf-set", type=string2bspec_tuple,
help='Manual specification of VSWF set codes (format as a python list of integers); see docs on qpms_uvswfi_t for valid codes or simply use --lMax instead. Overrides --lMax.')
grp.add_argument("-L", "--lMax", type=int, default=3, help='multipole degree cutoff')
atomic_arguments = {
'rectlattice2d_periods': lambda ap: ap.add_argument("-p", "--period", type=float, nargs='+', required=True, help='square/rectangular lattice periods', metavar=('px','[py]')),
@ -264,6 +344,7 @@ class ArgParser:
'bg_real_refractive_index': lambda ap: ap.add_argument("-n", "--refractive-index", type=float, default=1., help='background medium strictly real refractive index'),
'bg_analytical': lambda ap: ap.add_argument("-B", "--background", type=material_spec, default=1., help="Background medium specification (constant real or complex refractive index, or supported material label)"),
'single_lMax': lambda ap: ap.add_argument("-L", "--lMax", type=int, required=True, default=3, help='multipole degree cutoff'),
'single_lMax_or_vswfset': __add_single_lMax_or_vswfset_group,
'single_lMax_extend': lambda ap: ap.add_argument("--lMax-extend", type=int, required=False, default=6, help='multipole degree cutoff for T-matrix calculation (cylindrical particles only'),
'outfile': lambda ap: ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)'), # TODO consider type=argparse.FileType('w')
'plot_out': lambda ap: ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)"),
@ -280,6 +361,7 @@ class ArgParser:
'single_particle': ("Single particle definition (shape [currently spherical or cylindrical]) and materials, incl. background)", ('background',), ('single_material', 'single_radius', 'single_height', 'single_lMax_extend'), ('_eval_single_tmgen',), ()),
'multi_particle': ("One or more particle definition (shape [curently spherical or cylindrical]), materials, and positions)", ('background',), ('multi_particle',), ('_process_multi_particle',), ()),
'single_lMax': ("Single particle lMax definition", (), ('single_lMax',), (), ()),
'single_lMax_or_vswfset': ("Single particle lMax definition", (), ('single_lMax_or_vswfset',), ('_eval_single_lMax_or_vswfset',), ()),
'single_omega': ("Single angular frequency", (), ('single_frequency_eV',), ('_eval_single_omega',), ()),
'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',), ()),
'omega_seq_real_ng': ("Equidistant real frequency ranges or individual frequencies (new syntax)", (), ('real_frequencies_eV_ng',), ('_eval_omega_seq_real_ng',), ()),
@ -325,18 +407,29 @@ class ArgParser:
if tmgspec in self._tmg_register.keys():
return self._tmg_register[tmgspec]
else:
from .cytmatrices import TMatrixGenerator
bgspec, fgspec, (radius, height, lMax_extend) = tmgspec
bg = self._add_emg(bgspec)
fg = self._add_emg(fgspec)
if height is None:
tmgen = TMatrixGenerator.sphere(bg, fg, radius)
from .cytmatrices import TMatrixGenerator, CTMatrix
if isinstance(tmgspec, constant_tmatrix_spec):
tmgen = TMatrixGenerator(CTMatrix(tmgspec.bspec, tmgspec.matrix))
else:
tmgen = TMatrixGenerator.cylinder(bg, fg, radius, height, lMax_extend=lMax_extend)
bgspec, fgspec, (radius, height, lMax_extend) = tmgspec
bg = self._add_emg(bgspec)
fg = self._add_emg(fgspec)
if height is None:
tmgen = TMatrixGenerator.sphere(bg, fg, radius)
else:
tmgen = TMatrixGenerator.cylinder(bg, fg, radius, height, lMax_extend=lMax_extend)
self._tmg_register[tmgspec] = tmgen
return tmgen
def _add_bspec(self, key):
"""
Transforms a given key into a BaseSpec and registers
the BaseSpec with the key.
Always returns a BaseSpec (unless an exception occurs).
If an equivalent instance of BaseSpec is registered, returns that one
(therefore, all equivalent BaseSpecs in the register values should
be the same instance).
"""
if key in self._bspec_register.keys():
return self._bspec_register[key]
else:
@ -345,7 +438,9 @@ class ArgParser:
bspec = key
elif isinstance(key, int):
bspec = self._add_bspec(BaseSpec(lMax=key))
else: raise TypeError("Can't register this as a BaseSpec")
#else: raise TypeError("Can't register this as a BaseSpec")
else:
bspec = self._add_bspec(BaseSpec(key))
self._bspec_register[key] = bspec
return bspec
@ -406,13 +501,24 @@ class ArgParser:
from .cymaterials import EpsMuGenerator, lorentz_drude
from .cytmatrices import TMatrixGenerator
self.foreground_emg = self._add_emg(a.material)
self.tmgen = self._add_tmg((a.background, a.material, (a.radius, a.height, a.lMax_extend)))
self.tmgen = self._add_tmg(tmgen_spec(a.background, a.material,
cyl_sph_dimensions(a.radius, a.height, a.lMax_extend)))
self.bspec = self._add_bspec(a.lMax)
def _eval_single_omega(self): # feature: single_omega
from .constants import eV, hbar
self.omega = self.args.eV * eV / hbar
def _eval_single_lMax_or_vswfset(self): #feature: single_lMax_or_vswfset
from .cybspec import BaseSpec
a = self.args
if a.vswf_set is not None:
self.bspec = BaseSpec(a.vswf_set)
self.bspecstr = "WF%s" % (str(a.vswf_set),)
else:
self.bspec = BaseSpec(lMax=a.lMax)
self.bspecstr = "L%d" % (a.lMax,)
def _eval_omega_seq(self): # feature: omega_seq
import numpy as np
from .constants import eV, hbar
@ -485,30 +591,65 @@ class ArgParser:
self.positions = {}
pos13, pos23, pos33 = False, False, False # used to
if len(a.position.keys()) == 0:
warnings.warn("No particle position (-p or +p) specified, assuming single particle in the origin / single particle per unit cell!")
warnings.warn("No particle position (-p or +p) specified, assuming single particle in the "
"origin / single particle per unit cell!")
a.position[None] = [(0.,0.,0.)]
for poslabel in a.position.keys():
_resolve_exclusive_groups(a, poslabel, ("lMax",), ("vswf_set",))
try:
lMax = a.lMax.get(poslabel, False) or a.lMax[None]
radius = a.radius.get(poslabel, False) or a.radius[None]
# Height is "inherited" only together with radius
height = a.height.get(poslabel, None) if poslabel in a.radius.keys() else a.height.get(None, None)
if hasattr(a, 'lMax_extend'):
lMax_extend = a.lMax_extend.get(poslabel, False) or a.lMax_extend.get(None, False) or None
else:
lMax_extend = None
material = a.material.get(poslabel, False) or a.material[None]
lMax_or_bspec = ( a.vswf_set.get(poslabel, False)
or a.lMax.get(poslabel, False)
or a.vswf_set.get(None, False)
or a.lMax[None] )
bspec = self._add_bspec(lMax_or_bspec)
self.bspecs[poslabel] = bspec
except (TypeError, KeyError) as exc:
if poslabel is None:
raise ArgumentProcessingError("Unlabeled particles' positions (-p) specified, but some default particle properties are missing (--lMax, --radius, and --material have to be specified)") from exc
raise ArgumentProcessingError("Unlabeled particles' positions (-p) specified, "
"but neither --lMax nor --vswf-set is specified") from exc
else:
raise ArgumentProcessingError(("Incomplete specification of '%s'-labeled particles: you must"
"provide at least ++lMax, ++radius, ++material arguments with the label, or the fallback arguments"
"--lMax, --radius, --material.")%(str(poslabel),)) from exc
tmspec = (a.background, material, (radius, height, lMax_extend))
"provide either ++lMax or ++vswf-set argument with the label, "
"or one of the fallback arguments --lMax or --vswf-set."
)%(str(poslabel),)) from exc
agi = _resolve_exclusive_groups(a, poslabel,
("constant_tmatrix",),
("radius", "height", "material", "lMax_extend"),
)
if agi == 0: # constant T-matrix
tmspec = constant_tmatrix_spec(bspec,
hashable_complex_matrix(a.constant_tmatrix.get(poslabel, None)
or a.constant_tmatrix[None]))
elif agi == 1: #
try:
radius = a.radius.get(poslabel, False) or a.radius[None]
# Height is "inherited" only together with radius
height = a.height.get(poslabel, None) if poslabel in a.radius.keys() else a.height.get(None, None)
if hasattr(a, 'lMax_extend'):
lMax_extend = a.lMax_extend.get(poslabel, False) or a.lMax_extend.get(None, False) or None
else:
lMax_extend = None
material = a.material.get(poslabel, False) or a.material[None]
except (TypeError, KeyError) as exc:
if poslabel is None:
raise ArgumentProcessingError("Unlabeled particles' positions (-p) specified, "
"but some default particle properties are missing "
"(either --lMax or --vswf-set, and --constant-tmatrix or "
"--radius and --material have to be specified)") from exc
else:
raise ArgumentProcessingError(("Incomplete specification of '%s'-labeled particles: you must"
"provide at least ++lMax (or ++vswf-set), "
"++radius and ++material (or ++constant-tmatrix) arguments with the label, "
"or the fallback arguments --lMax (or --vswf-sit), --radius, --material "
"(or --constant-tmatrix)."
)%(str(poslabel),)) from exc
tmspec = tmgen_spec(a.background, material,
cyl_sph_dimensions(radius, height, lMax_extend))
else: raise AssertionErrorw
self.tmspecs[poslabel] = tmspec
self.tmgens[poslabel] = self._add_tmg(tmspec)
self.bspecs[poslabel] = self._add_bspec(lMax)
poslist_cured = []
for pos in a.position[poslabel]:
if len(pos) == 1:

View File

@ -3,18 +3,20 @@
*/
#ifndef BEYN_H
#define BEYN_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
#include <complex.h>
/// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
/** Pure C array version */
typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
typedef int (*beyn_function_M_t)(_Complex double *target_M, size_t m, _Complex double z, void *params);
/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
/** Pure C array version */
typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
const complex double *Vhat, complex double z, void *params);
typedef int (*beyn_function_M_inv_Vhat_t)(_Complex double *target_M_inv_Vhat, size_t m, size_t l,
const _Complex double *Vhat, _Complex double z, void *params);
/// Complex plane integration contour structure.
typedef struct beyn_contour_t {
@ -26,15 +28,15 @@ typedef struct beyn_contour_t {
* It does not have to be a centre in some strictly defined sense,
* but it should be "somewhere around" where the contour is.
*/
complex double centre;
_Complex double centre;
/// Function testing that a point \a z lies inside the contour (optional).
_Bool (*inside_test)(struct beyn_contour_t *, complex double z);
complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
_Bool (*inside_test)(struct beyn_contour_t *, _Complex double z);
_Complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
} beyn_contour_t;
/// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */
beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
beyn_contour_t *beyn_contour_ellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints);
typedef enum {
@ -47,23 +49,23 @@ typedef enum {
/// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */
beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
beyn_contour_halfellipse_orientation or);
beyn_contour_t *beyn_contour_halfellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints,
beyn_contour_halfellipse_orientation o);
/// Similar to halfellipse but with rounded corners.
beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im,
beyn_contour_t *beyn_contour_kidney(_Complex double centre, double halfax_re, double halfax_im,
double rounding, ///< Must be in interval [0, 0.5)
size_t n, beyn_contour_halfellipse_orientation or);
size_t n, beyn_contour_halfellipse_orientation o);
/// Beyn algorithm result structure (pure C array version).
typedef struct beyn_result_t {
size_t neig; ///< Number of eigenvalues found.
size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
complex double *eigval;
complex double *eigval_err;
_Complex double *eigval;
_Complex double *eigval_err;
double *residuals;
complex double *eigvec; // Rows are the eigenvectors
_Complex double *eigvec; // Rows are the eigenvectors
double *ranktest_SV;
} beyn_result_t;
@ -83,4 +85,7 @@ beyn_result_t *beyn_solve(
double res_tol ///< (default: `0.0`) TODO DOC.
);
#ifdef __cplusplus
}
#endif
#endif // BEYN_H

View File

@ -1,4 +1,5 @@
import numpy as np
cimport numpy as np
import enum
from .cycommon import get_mn_y, tlm2uvswfi
__VSWF_norm_dict = {
@ -35,6 +36,13 @@ cdef class BaseSpec:
whenever used in other wrapper classes that need the pointer
to qpms_vswf_set_spec_t, remember to set a (private, probably immutable) reference to qpms.basespec to ensure
correct reference counting and garbage collection.
BaseSpec can be constructed either using a single keyword argument `lMax` to include
all the transversal vector spherical wavefunction to a specific degree, or
using an explicit list of integer basis VSWF codes, \see qpms_uvswfi_t
BaseSpec(lMax = <positive_integer>)
or
BaseSpec([ uvswfi0, ... ])
'''
#cdef qpms_vswf_set_spec_t s # in pxd
#cdef np.ndarray __ilist # in pxd
@ -132,6 +140,14 @@ cdef class BaseSpec:
# TODO raise correct errors (TypeError on bad type of key, IndexError on exceeding index)
return self.__ilist[key]
def __repr__(self):
cdef qpms_l_t lMax_candidate
lMax_candidate = qpms_nelem2lMax(len(self.__ilist) // 2)
if lMax_candidate >= 0:
if np.array_equal(self.__ilist, BaseSpec(lMax=lMax_candidate).ilist):
return 'BaseSpec(lMax=%d)' % lMax_candidate
return 'BaseSpec([' + ','.join('%d' % i for i in self.__ilist) + '])'
property ilist:
def __get__(self):
return self.__ilist

View File

@ -20,6 +20,12 @@ class BesselType(enum.IntEnum):
HANKEL_PLUS = QPMS_HANKEL_PLUS
HANKEL_MINUS = QPMS_HANKEL_MINUS
class EwaldPart(enum.IntEnum):
LONG_RANGE = QPMS_EWALD_LONG_RANGE
SHORT_RANGE = QPMS_EWALD_SHORT_RANGE
FULL = QPMS_EWALD_FULL
ZEROTERM = QPMS_EWALD_0TERM
class PointGroupClass(enum.IntEnum):
CN = QPMS_PGS_CN
S2N = QPMS_PGS_S2N
@ -138,6 +144,46 @@ def get_y_mn_unsigned(int nmax):
i = i + 1
return(ymn_plus, ymn_minus)
def get_nelem_scalar(lMax):
# TODO DOC, exceptions etc.
#return qpms_lMax2nelem(lMax)
return lMax * (lMax + 2) + 1 # = (lMax + 1)**2
## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax
@cython.boundscheck(False)
def get_mn_y_scalar(int nmax):
"""
Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing;
inc. nmax.
('y to mn' conversion)
Parameters
----------
nmax : int
The maximum order to which the SSWFs / Legendre functions etc. will be evaluated.
Returns
-------
output : (m, n)
Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax + 1), dtype=np.int),
where [(m[y],n[y]) for y in range(nmax*nmax + 2*nmax + 1)] covers all possible
integer pairs n >= 0, -n <= m <= n.
"""
cdef Py_ssize_t nelems = (nmax + 1)**2
cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int)
cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int)
cdef Py_ssize_t i = 0
cdef np.int_t n, m
for n in range(0,nmax+1):
for m in range(-n,n+1):
m_arr[i] = m
n_arr[i] = n
i = i + 1
return (m_arr, n_arr)
def tlm2uvswfi(t, l, m):
''' TODO doc

View File

@ -1,12 +1,36 @@
from .qpms_cdefs cimport *
from .qpms_c cimport ScatteringSystem, _ScatteringSystemAtOmega, _ScatteringSystemAtOmegaK
from .cycommon import EwaldPart, BesselType
from libc.stdlib cimport malloc, free, calloc
import numpy as np
from cython.parallel cimport prange, parallel
from cython cimport boundscheck, wraparound
cdef extern from "ewald.h":
void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z)
void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z)
void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z, bint bigimz)
int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result)
extern int ewald_factor_ipow_l
extern int ewald_factor_ipow_m
cdef extern from "translations_dbg.h":
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
cdouble *dest, const csph_t kdlj, const qpms_bessel_t J) nogil
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
cdouble *dest, const cdouble k, const cart3_t r, const qpms_bessel_t J) nogil
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
cdouble * const sigmas_total, double * serr_total,
const double eta, const cdouble k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts) nogil
cdef extern from "scatsystem_dbg.h":
int qpms_scatsyswk_test_sswf_basis_e(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, cart3_t evalpoint, qpms_ewald_part parts) nogil
def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'):
cdef np.ndarray[double, ndim=1] err_np
@ -29,9 +53,117 @@ def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, met
else:
return target_np
def scatsyswk_sswf_basis(_ScatteringSystemAtOmegaK sswk, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
"""Evaluate Ewald-summed scalar SWFs that are used to compose the tranlation operators.
Parameters
----------
evalpos: array_like of floats and shape (..., 3)
Evaluation points in cartesian coordinates.
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (particle_count, pnelem)`
"""
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 qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
cdef ScatteringSystem ss = sswk.ssw_pyref.ss_pyref
cdef qpms_scatsys_t *ss_c = ss.s
cdef qpms_scatsys_at_omega_k_t *sswk_c = sswk.rawpointer()
cdef const qpms_trans_calculator *c = ss_c[0].c
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
cdef qpms_ss_pi_t particle_count = sswk.ssw_pyref.ss_pyref.particle_count
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], particle_count, pnelem), dtype=complex)
cdef cdouble *res
cdef cart3_t pos
cdef Py_ssize_t i, pi, j
with nogil, wraparound(False), parallel():
res = <cdouble *> malloc(particle_count * pnelem * sizeof(cdouble))
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_test_sswf_basis_e(res, sswk_c, btyp_c, pos, ewaldparts_c)
for pi in range(particle_count):
for j in range(pnelem):
results[i,pi,j] = res[pi * pnelem + j]
free(res)
return results.reshape(evalpos.shape[:-1] + (particle_count, pnelem))
def scatsysw_eval_sswf(_ScatteringSystemAtOmega ssw, evalpos, btyp=QPMS_HANKEL_PLUS):
"""Evaluate sswfs by internal qpms_trans_calculator at given locations.
Parameters
----------
ssw: _ScatteringSystemAtOmega
We take the wavenumber and qpms_trans_calculator from this, nothing else.
(In a way, this is an overkill, but there is no up-to-date standalone cython
wrapper for qpms_trans_calculator.)
evalpos: array_like of floats and shape (..., 3)
Evaluation points in cartesian coordinates.
btyp: BesselType
Kind of Bessel functions to be used
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (pnelem,)`
where `pnelem` is TODO
"""
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 ScatteringSystem ss = ssw.ss_pyref
cdef qpms_scatsys_t *ss_c = ss.s
cdef const qpms_trans_calculator *c = ss_c[0].c
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
cdef cdouble k = ssw.wavenumber
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[complex,ndim=2] results = np.empty((evalpos_a.shape[0], pnelem), dtype=complex)
cdef cdouble *res
cdef cart3_t pos
cdef csph_t kdlj
cdef Py_ssize_t i, j
with nogil, wraparound(False), parallel():
res = <cdouble *> malloc(pnelem * sizeof(cdouble))
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]
kdlj = cart2csph(pos)
kdlj.r *= k
qpms_trans_calculator_test_sswf(c, res, kdlj, btyp_c)
for j in range(pnelem):
results[i,j] = res[j]
free(res)
return results.reshape(evalpos.shape[:-1] + (pnelem,))
def gamma_inc(double a, cdouble x, int xbranch=0):
cdef qpms_csf_result res
complex_gamma_inc_e(a, x, xbranch, &res)
return res.val
def ipow_l_factor(n = None):
global ewald_factor_ipow_l
if n is not None:
ewald_factor_ipow_l = n
return ewald_factor_ipow_l
def ipow_m_factor(n = None):
global ewald_factor_ipow_m
if n is not None:
ewald_factor_ipow_m = n
return ewald_factor_ipow_m

View File

@ -286,6 +286,10 @@ cdef class TMatrixGenerator:
self.holder = what
self.g.function = qpms_tmatrix_generator_interpolator
self.g.params = <void*>(<TMatrixInterpolator?>self.holder).rawpointer()
elif what == 0:
self.holder = 0
self.g.function = qpms_tmatrix_generator_zero
self.g.params = <void*> 0
elif callable(what):
warnings.warn("Custom python T-matrix generators are an experimental feature. Also expect it to be slow.")
self.holder = what
@ -403,3 +407,15 @@ cdef class TMatrixGenerator:
EpsMuGenerator(outside), EpsMuGenerator(inside),
ArcFunction(__ArcCylinder(r, h)), *args, **kwargs))
@staticmethod
def dummy():
"""Returns a dummy T-matrix generator (returns a zero T-matrix).
Returns
-------
tmgen_dummy : TMatrixGenerator
"""
return tmgen_dummy
# pre-generate a dummy TMatrixGenerator (which returns zero T-matrix)
tmgen_dummy = TMatrixGenerator(0)

View File

@ -7,11 +7,23 @@
#include <complex.h>
#include "tiny_inlines.h"
#include "qpms_error.h"
#include "lattices.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_expint.h>
/*
* Control additional phase factors
* that could possibly affect the resulting VSWF
* shape.
* Currently not used in ewald3_1_z_sigma_long()
*/
int ewald_factor_ipow_l = 0;
int ewald_factor_ipow_m = 0;
// parameters for the quadrature of integral in (4.6)
#ifndef INTEGRATION_WORKSPACE_LIMIT
#define INTEGRATION_WORKSPACE_LIMIT 30000
@ -125,7 +137,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->s1_constfacs_base = malloc(s1_constfacs_sz * sizeof(complex double));
size_t s1_constfacs_sz_cumsum = 0;
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
if ((m + n) % 2 == 0) {
c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum;
@ -161,7 +173,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->S1_constfacs = malloc((1+c->nelem_sc) * sizeof(complex double *));
//determine sizes
size_t S1_constfacs_sz = 0;
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
const qpms_l_t L_M = n - abs(m);
for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
@ -176,7 +188,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->S1_constfacs_base = malloc(S1_constfacs_sz * sizeof(complex double));
size_t S1_constfacs_sz_cumsum = 0; // second count
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
const complex double yfactor = -2 * ipow(n+1) * M_SQRTPI
* factorial((n-m)/2) * factorial((n+m)/2);
@ -303,7 +315,7 @@ int ewald3_21_xy_sigma_long (
const bool k_is_real = (cimag(k) == 0);
assert((latdim & LAT_XYONLY) && (latdim & SPACE3D));
assert((latdim & LAT1D) || (latdim & LAT2D));
const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
assert(nelem_sc > 0);
const qpms_l_t lMax = c->lMax;
@ -435,7 +447,7 @@ int ewald3_21_xy_sigma_long (
for(qpms_m_t m = -n; m <= n; ++m) {
if((particle_shift.z == 0) && ((m+n) % 2 != 0)) // odd coefficients are zero.
continue;
const qpms_y_t y = qpms_mn2y_sc(m, n);
const qpms_y_sc_t y = qpms_mn2y_sc(m, n);
size_t constidx = 0; // constants offset
const complex double e_imalpha_pq = cexp(I*m*arg_pq);
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
@ -444,7 +456,7 @@ int ewald3_21_xy_sigma_long (
assert((n-abs(m))/2 == c->s1_jMaxes[y]);
for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
complex double summand = rbeta_pq_div_k_pow[n-2*j]
* e_imalpha_pq * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
* e_imalpha_pq * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow(m) // This line can actually go outside j-loop
* gamma_pq_powm1[2*j]// * Gamma_pq[j] bellow (GGG) after error computation
* c->s1_constfacs[y][j];
if(err) {
@ -455,6 +467,7 @@ int ewald3_21_xy_sigma_long (
ckahanadd(&jsum, &jsum_c, summand);
}
jsum *= phasefac * factor1d; // PFC
jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, jsum);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(jsum));
@ -483,6 +496,7 @@ int ewald3_21_xy_sigma_long (
ckahanadd(&jsum, &jsum_c, jsummand);
}
jsum *= phasefac; // factor1d not here, off-axis sums not implemented/allowed.
jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, jsum);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(jsum));
@ -496,7 +510,7 @@ int ewald3_21_xy_sigma_long (
#ifdef EWALD_AUTO_CUTOFF
++li;
double cursum_min = INFINITY;
for (qpms_y_t y = 0; y < nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
const double a = cabs(target[y]);
if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
cursum_min = MIN(cursum_min, a);
@ -516,10 +530,10 @@ ewald3_21_xy_sigma_long_end_point_loop:
free(err_c);
free(target_c);
for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
target[y] *= commonfac;
if(err)
for(qpms_y_t y = 0; y < nelem_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
err[y] *= commonfac;
return 0;
}
@ -548,8 +562,8 @@ int ewald3_1_z_sigma_long (
assert(beta.x == 0 && beta.y == 0);
assert(particle_shift.x == 0 && particle_shift.y == 0);
const double beta_z = beta.z;
const double particle_shift_z = particle_shift_z;
const qpms_y_t nelem_sc = c->nelem_sc;
const double particle_shift_z = particle_shift.z;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax;
// Manual init of the ewald summation targets
@ -600,7 +614,7 @@ int ewald3_1_z_sigma_long (
// TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
// and just fetched for each n
for(qpms_l_t n = 0; n <= lMax; ++n) {
const qpms_y_t y = qpms_mn2y_sc(0, n);
const qpms_y_sc_t y = qpms_mn2y_sc(0, n);
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
for(qpms_l_t j = 0; j <= n/2; ++j) {
@ -634,10 +648,10 @@ int ewald3_1_z_sigma_long (
free(err_c);
free(target_c);
for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
target[y] *= commonfac;
if(err)
for(qpms_y_t y = 0; y < nelem_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
err[y] *= commonfac;
return 0;
}
@ -780,7 +794,7 @@ int ewald3_sigma_short(
{
const bool k_is_real = (cimag(k) == 0); // TODO check how the compiler optimises the loops
const double kreal = creal(k);
const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax;
gsl_integration_workspace *workspace =
gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
@ -902,7 +916,8 @@ int ewald3_sigma_short(
if(err)
kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
* interr[n])); // TODO include also other errors
complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow_m_neg(m);
complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow(m);
thesummand *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, thesummand);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(thesummand));
@ -916,7 +931,7 @@ int ewald3_sigma_short(
#ifdef EWALD_AUTO_CUTOFF
++li;
double cursum_min = INFINITY;
for (qpms_y_t y = 0; y < nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
const double a = cabs(target[y]);
if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
cursum_min = MIN(cursum_min, a);

View File

@ -26,22 +26,18 @@
#ifndef EWALD_H
#define EWALD_H
#ifdef __cplusplus
extern "C" {
#endif
#include <gsl/gsl_sf_result.h>
#include <stdlib.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_errno.h>
#include <math.h> // for inlined lilgamma
#include <complex.h>
#include "qpms_types.h"
#include "lattices.h"
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
#include "lattices_types.h"
#include <complex.h>
/// Use this handler to ignore underflows of incomplete gamma.
@ -55,11 +51,11 @@ gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
*/
typedef struct qpms_ewald3_constants_t {
qpms_l_t lMax;
qpms_y_t nelem_sc;
qpms_y_sc_t nelem_sc;
/// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`.
qpms_l_t *s1_jMaxes;
/// The constant factors for the long range part of a 2D Ewald sum.
complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
_Complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN:
*
@ -73,7 +69,7 @@ typedef struct qpms_ewald3_constants_t {
*
* s1_constfacs[y(m,n)][j] = 0
*/
complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
_Complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
// similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0)
/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
/** If the summation points lie along a different direction, use the formula for
@ -86,7 +82,7 @@ typedef struct qpms_ewald3_constants_t {
// TODO indexing mechanisms
/// The constant factors for the long range part of a 2D Ewald sum.
complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
_Complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN:
*
@ -100,7 +96,7 @@ typedef struct qpms_ewald3_constants_t {
*
* S1_constfacs[y(m,n)][j] = 0
*/
complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
_Complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
/** If the summation points lie along a different direction, use the formula for
* 2D sum with additional factor of
@ -108,7 +104,7 @@ typedef struct qpms_ewald3_constants_t {
*/
complex double **s1_constfacs_1Dz;
_Complex double **s1_constfacs_1Dz;
/* These are the actual numbers now:
* s1_constfacs_1Dz[n][j] =
*
@ -116,7 +112,7 @@ typedef struct qpms_ewald3_constants_t {
* --------------------------
* j! * 2**(2*j) * (n - 2*j)!
*/
complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
_Complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
* what the multipliers from translations.c count with.
@ -139,13 +135,13 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
/// Structure for holding complex-valued result of computation and an error estimate.
/** Similar to gsl_sf_result, but with complex val. */
typedef struct qpms_csf_result {
complex double val; ///< Calculation result.
_Complex double val; ///< Calculation result.
double err; ///< Error estimate.
} qpms_csf_result;
// [1, (A.9)]
static inline complex double lilgamma(double t) {
static inline _Complex double lilgamma(double t) {
t = fabs(t);
if (t >= 1)
return sqrt(t*t - 1);
@ -154,8 +150,8 @@ static inline complex double lilgamma(double t) {
}
// [1, (A.8)], complex version of lilgamma()
static inline complex double clilgamma(complex double z) {
complex double a1 = z - 1, a2 = z + 1;
static inline _Complex double clilgamma(_Complex double z) {
_Complex double a1 = z - 1, a2 = z + 1;
// ensure -pi/2 < arg(z + 1) < 3*pi/2
if (creal(a2) < 0 && cimag(a2) <= 0)
a2 = -csqrt(a2);
@ -180,7 +176,7 @@ static inline complex double clilgamma(complex double z) {
* even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`.
*/
int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
int cx_gamma_inc_series_e(double a, _Complex double z, qpms_csf_result * result);
/// Incomplete Gamma function as continued fractions.
/**
@ -192,7 +188,7 @@ int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
* even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`.
*/
int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
int cx_gamma_inc_CF_e(double a, _Complex double z, qpms_csf_result * result);
/// Incomplete gamma for complex second argument.
/**
@ -209,7 +205,7 @@ int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
* Another than principal branch can be selected using non-zero \a m
* argument.
*/
int complex_gamma_inc_e(double a, complex double x,
int complex_gamma_inc_e(double a, _Complex double x,
/// Branch index.
/** If zero, the principal value is calculated.
* Other branches might be chosen using non-zero \a m.
@ -226,7 +222,7 @@ int complex_gamma_inc_e(double a, complex double x,
/// Exponential integral for complex second argument.
/** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */
int complex_expint_n_e(int n, complex double x, qpms_csf_result *result);
int complex_expint_n_e(int n, _Complex double x, qpms_csf_result *result);
/// Hypergeometric 2F2, used to calculate some errors.
@ -245,15 +241,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
* unsuitable especially for big values of \a maxn.
*
*/
void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z);
void ewald3_2_sigma_long_Delta(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
/** This function always uses Kambe's (corrected) recurrent formula.
* For production, use ewald3_2_sigma_long_Delta() instead.
*/
void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z, _Bool bigimz);
void ewald3_2_sigma_long_Delta_recurrent(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z, _Bool bigimz);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
/** This function always uses Taylor expansion in \a z.
@ -263,26 +259,26 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_
* parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
* of the error.
*/
void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z);
void ewald3_2_sigma_long_Delta_series(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z);
// General functions acc. to [2], sec. 4.6 currently valid for 2D and 1D lattices in 3D space
/// The Ewald sum "self-interaction" term that appears in the lattice sums with zero (direct-space) Bravais lattice displacement.
int ewald3_sigma0(complex double *result, ///< Pointer to save the result (single complex double).
int ewald3_sigma0(_Complex double *result, ///< Pointer to save the result (single _Complex double).
double *err, ///< Pointer to save the result error estimate (single double).
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber ///< Wavenumber of the background medium.
_Complex double wavenumber ///< Wavenumber of the background medium.
);
/// Short-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{S}(\vect k,\vect s)\f$.
int ewald3_sigma_short(
complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
_Complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
double *target_sigmasr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber, ///< Wavenumber of the background medium.
_Complex double wavenumber, ///< Wavenumber of the background medium.
/// Lattice dimensionality.
/** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
LatticeDimensionality latdim,
@ -293,7 +289,7 @@ int ewald3_sigma_short(
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_R,
struct PGen *pgen_R,
/// Indicates whether pgen_R already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(),
* so the function assumes that the generated points correspond to the unshifted Bravais lattice,
@ -310,11 +306,11 @@ int ewald3_sigma_short(
/// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$.
int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
_Complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
double *target_sigmalr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber, ///< Wavenumber of the background medium.
_Complex double wavenumber, ///< Wavenumber of the background medium.
double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality).
/// Lattice dimensionality.
LatticeDimensionality latdim,
@ -325,7 +321,7 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_K,
struct PGen *pgen_K,
/// Indicates whether pgen_K already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(),
* so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
@ -339,4 +335,13 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
cart3_t particle_shift
);
// If nonzero, adds an additional factor \f$ i^{nl} \f$ to the Ewald sum result (for debugging).
extern int ewald_factor_ipow_l;
// If nonzero, adds an additional factor \f$ i^{nm} \f$ to the Ewald sum result (for debubbing).
extern int ewald_factor_ipow_m;
#ifdef __cplusplus
}
#endif
#endif //EWALD_H

View File

@ -14,6 +14,7 @@
#include <stdbool.h>
#include <Faddeeva.h>
#include "tiny_inlines.h"
#include "qpms_error.h"
// Some magic constants

View File

@ -3,6 +3,10 @@
*/
#ifndef GAUNT_H
#define GAUNT_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stdlib.h>
#define _GAUNT_H_MIN(x,y) (((x) > (y)) ? (y) : (x))
@ -30,4 +34,7 @@ double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu);
int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu);
#endif //GAUNT_PRECOMPILED
#ifdef __cplusplus
}
#endif
#endif //GAUNT_H

View File

@ -23,6 +23,9 @@
*/
#ifndef QPMS_GROUPS_H
#define QPMS_GROUPS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <assert.h>
@ -35,7 +38,7 @@ struct qpms_finite_group_irrep_t {
/** The r-th row, c-th column of the representation of the i'th element is retrieved as
* m[i * dim * dim + r * dim + c]
*/
complex double *m;
_Complex double *m;
};
/// A point group with its irreducible representations and some metadata.
@ -93,4 +96,7 @@ qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *na
extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL;
extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL_G;
#ifdef __cplusplus
}
#endif
#endif // QPMS_GROUPS_H

View File

@ -1,57 +1,106 @@
/*! \file indexing.h
* \brief Various index conversion functions.
*
* Index conventions in QPMS
* -------------------------
*
* Depending on circumstances (e.g. whether scalar or vector spherical wavefunctions
* are being used), QPMS uses two mappings between the multipole degree/order index
* pair (\a l, \a m), due to historical reasons ocassionally also noted as (\a n, \a m),
* and a flat non-negative integer index.
*
* Both mappings map the (\a l, \a m) pair in lexicographic order (with
* the basic constraint \f$ |m| \le l \f$), the difference is that the "scalar"
* mapping starts from \a l = 0, whereas the "vector" mapping starts from \a l = 1,
* so that for the zeroth elements we have
* * (0, 0) 0, (1, -1) 1 in the "scalar" mapping,
* * (1, -1) 0 in the "vector" mapping.
*
* The vector SWF flat index is typically denoted by `y` and its type is \ref qpms_y_t,
* whereas * the scalar SWF flat index is typically denoted by `y_sc` and its
* type is \ref qpms_y_sc_t.
*
* Moreover, there is another mapping that encodes also the electric/magnetic/scalar
* WSWF type index, see \ref qpms_uvswfi_t.
*
*/
#ifndef QPMS_INDEXING_H
#define QPMS_INDEXING_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <math.h>
#include "optim.h"
static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) {
/// Mapping from multipole degree, order to flat index (degree ≥ 1)
static inline qpms_y_t qpms_mn2y(
qpms_m_t m, ///< multipole order
qpms_l_t n ///< multipole degree (a.k.a. \a l); must be greater than 0
) {
return n * (n + 1) + m - 1;
}
/// Mapping from flat index to multipole degree (degree ≥ 1)
static inline qpms_lm_t qpms_y2n(qpms_y_t y) {
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
return sqrt(y+1);
}
/// Mapping from flat index, multipole degree (≥ 1) to multipole order
static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) {
return y-qpms_mn2y(0,n);
}
/// Mapping from flat index to multipole degree, order (degree ≥ 1)
static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
*m=qpms_yn2m(y,*n=qpms_y2n(y));
}
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 1
static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2);
}
/// Inverse of qpms_lMax2nelem(); returns -1 if there is no exact match.
static inline qpms_l_t qpms_nelem2lMax(qpms_y_t nelem) {
qpms_l_t l = qpms_y2n(nelem);
if (l == -qpms_yn2m(nelem, l))
return l - 1;
else
return -1;
}
// Scalar versions: they have a place for the 0, 0 term in the beginning
static inline qpms_y_t qpms_mn2y_sc(qpms_m_t m, qpms_l_t n) {
/// Mapping from multipole degree, order to flat index (degree ≥ 0)
static inline qpms_y_sc_t qpms_mn2y_sc(
qpms_m_t m, ///< multipole order
qpms_l_t n ///< multipole degree (a.k.a. \a l); muste be greater or equal than 0
) {
return n * (n + 1) + m;
}
static inline qpms_lm_t qpms_y2n_sc(qpms_y_t y) {
/// Mapping from flat index to multipole degree (degree ≥ 0)
static inline qpms_lm_t qpms_y2n_sc(qpms_y_sc_t y) {
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
return sqrt(y);
}
static inline qpms_m_t qpms_yn2m_sc(qpms_y_t y, qpms_l_t n) {
/// Mapping from flat index, multipole degree (≥ 0) to multipole order
static inline qpms_m_t qpms_yn2m_sc(qpms_y_sc_t y, qpms_l_t n) {
return y-qpms_mn2y_sc(0,n);
}
static inline void qpms_y2mn_sc_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
/// Mapping from flat index to multipole degree, order (degree ≥ 0)
static inline void qpms_y2mn_sc_p(qpms_y_sc_t y, qpms_m_t *m, qpms_l_t *n){
*m=qpms_yn2m_sc(y,*n=qpms_y2n_sc(y));
}
static inline qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2) + 1;
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 0
static inline qpms_y_sc_t qpms_lMax2nelem_sc(qpms_l_t lmax){
return lmax * ((qpms_y_sc_t)lmax + 2) + 1;
}
// TODO maybe enable crashing / validity control by macro definitions...
@ -69,7 +118,7 @@ static const qpms_uvswfi_t QPMS_UI_L00 = 0;
/** Returns a non-zero value if the u value is invalid. */
static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
qpms_y_sc_t y_sc = u / 4;
qpms_y2mn_sc_p(y_sc, m, n);
// Test validity
@ -82,7 +131,7 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
/** Does *not* allow for longitudinal waves. */
static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_y_t *y) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
*y = u / 4 - 1;
if (QPMS_UNLIKELY(*t == 0 || *t == 3)) return QPMS_ERROR;
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@ -95,7 +144,7 @@ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
*/
static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_y_t *y) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
*y = u / 4 - 1;
if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR;
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@ -110,4 +159,7 @@ static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_INDEXING_H

View File

@ -3,8 +3,9 @@
*/
#ifndef KAHANSUM_H
#define KAHANSUM_H
#include <complex.h>
#ifdef __cplusplus
extern "C" {
#endif
static inline void kahaninit(double * const sum, double * const compensation) {
*sum = 0;
@ -19,16 +20,19 @@ static inline void kahanadd(double *sum, double *compensation, double input) {
}
static inline void ckahaninit(complex double * const sum, complex double * const compensation) {
static inline void ckahaninit(_Complex double * const sum, _Complex double * const compensation) {
*sum = 0;
*compensation = 0;
}
static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) {
complex double compensated_input = input - *compensation;
complex double nsum = *sum + compensated_input;
static inline void ckahanadd(_Complex double *sum, _Complex double *compensation, _Complex double input) {
_Complex double compensated_input = input - *compensation;
_Complex double nsum = *sum + compensated_input;
*compensation = (nsum - *sum) - compensated_input;
*sum = nsum;
}
#ifdef __cplusplus
}
#endif
#endif //KAHANSUM_H

View File

@ -621,10 +621,10 @@ typedef struct PGen_LatticeRadialHeap_StateData {
size_t heap_capacity;
double *r_heap;
int *coord_heap;
double b[0]; // basis vectors and offset
double offset_r;
double minR, maxR;
bool inc_minR, inc_maxR;
double b[0]; // basis vectors and offset
} PGen_LatticeRadialHeap_StateData;
static inline double nd2norm(const double a[], int d) {

View File

@ -1,9 +1,15 @@
/*! \file lattices.h
* \brief Lattice point generators and lattice vector analysis / transformation.
*
* \bug Header file not C++ compatible.
*/
#ifndef LATTICES_H
#define LATTICES_H
#ifdef __cplusplus // FIXME Not C++ compatible. Include "lattices_types.h" for minimal necessary enum decls.
extern "C" {
#endif
#include "lattices_types.h"
#include <math.h>
#include <stdbool.h>
@ -30,22 +36,6 @@
*
*/
typedef enum LatticeDimensionality {
LAT1D = 1,
LAT2D = 2,
LAT3D = 4,
SPACE1D = 8,
SPACE2D = 16,
SPACE3D = 32,
LAT_1D_IN_3D = 33,
LAT_2D_IN_3D = 34,
LAT_3D_IN_3D = 40,
// special coordinate arrangements (indicating possible optimisations)
LAT_ZONLY = 64,
LAT_XYONLY = 128,
LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
} LatticeDimensionality;
inline static bool LatticeDimensionality_checkflags(
LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) {
@ -945,4 +935,7 @@ int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxste
int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r);
void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g);
#ifdef __cplusplus
}
#endif
#endif // LATTICES_H

37
qpms/lattices_types.h Normal file
View File

@ -0,0 +1,37 @@
/*! \file lattices_types.h
* \brief Auxiliary lattice point generator and related enum declarations.
*
* This file contains necessary declarations needed in other header files.
* In contrast to lattices.h, this one is C++ compatible.
*
* \see lattices.h
*/
#ifndef LATTICES_TYPES_H
#define LATTICES_TYPES_H
#ifdef __cplusplus
extern "C" {
#endif
typedef enum LatticeDimensionality {
LAT1D = 1,
LAT2D = 2,
LAT3D = 4,
SPACE1D = 8,
SPACE2D = 16,
SPACE3D = 32,
LAT_1D_IN_3D = 33,
LAT_2D_IN_3D = 34,
LAT_3D_IN_3D = 40,
// special coordinate arrangements (indicating possible optimisations)
LAT_ZONLY = 64,
LAT_XYONLY = 128,
LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
} LatticeDimensionality;
struct PGen;
#ifdef __cplusplus
}
#endif
#endif // LATTICES_TYPES_H

View File

@ -5,6 +5,7 @@
#include <string.h>
#include "materials.h"
#include "qpms_error.h"
#include <stdbool.h>
#define SQ(x) ((x)*(x))

View File

@ -3,7 +3,12 @@
*/
#ifndef QPMS_MATERIALS_H
#define QPMS_MATERIALS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <complex.h>
#include <gsl/gsl_spline.h>
#ifndef SPEED_OF_LIGHT
@ -20,36 +25,36 @@ typedef struct qpms_epsmu_generator_t {
* qpms_permittivity_interpolator_epsmu_g(),
* qpms_lorentzdrude_epsmu_g().
*/
qpms_epsmu_t (*function) (complex double omega, const void *params);
qpms_epsmu_t (*function) (_Complex double omega, const void *params);
const void *params;
} qpms_epsmu_generator_t;
/// Convenience function for generating material properties at given frequency.
static inline qpms_epsmu_t qpms_epsmu_generator_eval(
qpms_epsmu_generator_t gen, complex double omega) {
qpms_epsmu_generator_t gen, _Complex double omega) {
return gen.function(omega, gen.params);
}
/// Constant optical property "generator" for qpms_epsmu_generator_t.
qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored.
qpms_epsmu_t qpms_epsmu_const_g(_Complex double omega, ///< Frequency ignored.
const void *epsmu ///< Points to the qpms_epsmu_t to be returned.
);
/// Gets refractive index of a material from its permeability and permittivity.
/** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */
static inline complex double qpms_refindex(qpms_epsmu_t em) {
static inline _Complex double qpms_refindex(qpms_epsmu_t em) {
return csqrt(em.eps * em.mu);
}
/// Gets wave number \a k from angular frequency and material permeability and permittivity.
/** \f[ k = \frac{n\omega}{c_0} = \frac{\omega\sqrt{\mu_r \varepsilon_r}}{c_0} \f] */
static inline complex double qpms_wavenumber(complex double omega, qpms_epsmu_t em) {
static inline _Complex double qpms_wavenumber(_Complex double omega, qpms_epsmu_t em) {
return qpms_refindex(em)*omega/SPEED_OF_LIGHT;
}
/// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity.
/** \eta_r = \sqrt{\mu_r / \varepsilon_r} \f] */
static inline complex double qpms_waveimpedance(qpms_epsmu_t em) {
static inline _Complex double qpms_waveimpedance(qpms_epsmu_t em) {
return csqrt(em.mu / em.eps);
}
@ -67,7 +72,7 @@ typedef struct qpms_ldparams_triple_t {
* \f]
*/
typedef struct qpms_ldparams_t {
complex double eps_inf; ///< Permittivity at infinity.
_Complex double eps_inf; ///< Permittivity at infinity.
double omega_p; ///< Plasma frequency.
size_t n; ///< Number of "oscillators".
qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters.
@ -86,14 +91,14 @@ extern const qpms_ldparams_t *const QPMS_LDPARAMS_PT; ///< Lorentz-Drude paramet
extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten.
/// Lorentz-Drude permittivity.
complex double qpms_lorentzdrude_eps(complex double omega, const qpms_ldparams_t *);
_Complex double qpms_lorentzdrude_eps(_Complex double omega, const qpms_ldparams_t *);
/// Lorentz-Drude optical properties, with relative permeability set always to one.
qpms_epsmu_t qpms_lorentzdrude_epsmu(complex double omega, const qpms_ldparams_t *);
qpms_epsmu_t qpms_lorentzdrude_epsmu(_Complex double omega, const qpms_ldparams_t *);
/// Lorentz-Drude optical properties, with relative permeability set always to one, compatible with qpms_epsmu_generator_t.
qpms_epsmu_t qpms_lorentzdrude_epsmu_g(
complex double omega,
_Complex double omega,
const void *ldparams ///< Lorentz-Drude parameters, in reality const qpms_ldparams_t *.
);
@ -125,14 +130,14 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
);
/// Evaluates interpolated material permittivity at a given angular frequency.
complex double qpms_permittivity_interpolator_eps_at_omega(
_Complex double qpms_permittivity_interpolator_eps_at_omega(
const qpms_permittivity_interpolator_t *interp, double omega_SI);
/// Evaluates interpolated material permittivity at a given angular frequency, qpms_epsmu_generator_t compatible version.
/** Permeability is always set to one. Imaginary part of omega is discarded.
*/
qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g(
complex double omega, ///< Angular frequency. The imaginary part is ignored!
_Complex double omega, ///< Angular frequency. The imaginary part is ignored!
const void * interpolator ///< Interpolator of type qpms_permittivity_interpolator_t
);
@ -148,14 +153,17 @@ double qpms_permittivity_interpolator_omega_max(
void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
/// Relative permittivity from the Drude model.
static inline complex double qpms_drude_epsilon(
complex double eps_inf, ///< Relative permittivity "at infinity".
complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
static inline _Complex double qpms_drude_epsilon(
_Complex double eps_inf, ///< Relative permittivity "at infinity".
_Complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
_Complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
_Complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
) {
return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_MATERIALS_H

View File

@ -2,14 +2,18 @@
* \brief Convention-dependent coefficients for VSWFs.
*
* See also @ref qpms_normalisation_t and @ref vswf_conventions.
*
* \bug Header file not C++ compatible.
*/
#ifndef NORMALISATION_H
#define NORMALISATION_H
#ifdef __cplusplus //FIXME not C++ compatible yet (enum bit operations)
extern "C" {
#endif
#include "qpms_types.h"
#include "qpms_error.h"
#include <math.h>
#include <complex.h>
#include "indexing.h"
#include "optim.h"
@ -36,8 +40,8 @@ static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qp
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -51,8 +55,8 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
@ -62,8 +66,8 @@ static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t no
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -77,14 +81,14 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
/// Returns the factors of a electric basis VSWF divided by the factor of a magnetic VWFS of a given convention, compared to the reference one.
static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
static inline _Complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
return qpms_normalisation_factor_N_noCS(norm, l, m)
/ qpms_normalisation_factor_M_noCS(norm, l, m);
}
@ -95,8 +99,8 @@ static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -109,13 +113,13 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
/// Returns the factors of a basis VSWF of a given convention compared to the reference convention.
static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
static inline _Complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l;
qpms_uvswfi2tmn(ui, &t, &m, &l);
switch(t) {
@ -156,7 +160,7 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t
* 0 \quad \mbox{if } m>0. \\
* \f]
*/
static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
static inline _Complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
case 0:
@ -194,7 +198,7 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor
*
*
*/
static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
static inline _Complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
if(m==0) return 0;
switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
@ -283,4 +287,7 @@ static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t
}
#endif
#ifdef __cplusplus
}
#endif
#endif //NORMALISATION_H

View File

@ -3,6 +3,9 @@
*/
#ifndef QPMS_PARSING_H
#define QPMS_PARSING_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
@ -57,4 +60,7 @@ size_t qpms_parse_doubles_fromfile(
const char *filepath //< File to read from, or NULL, "", "-" to read from stdin.
);
#ifdef __cplusplus
}
#endif
#endif // QPMS_PARSING_H

View File

@ -1,6 +1,7 @@
#include "pointgroups.h"
#include <search.h>
#include <stdlib.h>
#include <stdbool.h>
#include <assert.h>
double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;

View File

@ -3,6 +3,9 @@
*/
#ifndef POINTGROUPS_H
#define POINTGROUPS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_error.h"
#include "quaternions.h"
@ -18,9 +21,9 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
case QPMS_PGS_DN:
case QPMS_PGS_DND:
case QPMS_PGS_DNH:
return true;
return 1;
default:
return false;
return 0;
}
}
@ -81,4 +84,7 @@ _Bool qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b);
#ifdef __cplusplus
}
#endif
#endif //POINTGROUPS_H

View File

@ -1,4 +1,6 @@
from .qpms_cdefs cimport qpms_finite_group_t
from .qpms_cdefs cimport qpms_finite_group_t, qpms_scatsys_t, qpms_iri_t,\
qpms_scatsys_at_omega_t, qpms_scatsys_at_omega_k_t
from .cymaterials cimport EpsMuGenerator
cdef class FinitePointGroup:
cdef readonly bint owns_data
@ -6,3 +8,30 @@ cdef class FinitePointGroup:
cdef inline qpms_finite_group_t *rawpointer(self):
return self.G
cdef class ScatteringSystem:
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=?)
cdef class _ScatteringSystemAtOmegaK:
cdef qpms_scatsys_at_omega_k_t sswk
cdef _ScatteringSystemAtOmega ssw_pyref
cdef inline qpms_scatsys_at_omega_k_t *rawpointer(self):
return &self.sswk
cdef class _ScatteringSystemAtOmega:
cdef qpms_scatsys_at_omega_t *ssw
cdef ScatteringSystem ss_pyref
cdef inline qpms_scatsys_at_omega_t *rawpointer(self):
return self.ssw

0
qpms/qpms_c.py Normal file
View File

View File

@ -11,7 +11,7 @@ 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 .cycommon import string_c2py, PointGroupClass, BesselType, EwaldPart
from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
from .cymaterials cimport EpsMuGenerator, EpsMu
from libc.stdlib cimport malloc, free, calloc
@ -372,12 +372,6 @@ cdef class ScatteringSystem:
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
@ -702,7 +696,7 @@ cdef class ScatteringSystem:
self.s, iri, 0)
return target_np
def pack_matrix(self, fullmatrix, iri):
def pack_matrix(self, fullmatrix, iri, version='normal'):
"""Converts (projects) a matrix into an irrep subspace.
Parameters
@ -734,10 +728,14 @@ cdef class ScatteringSystem:
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],
if version == 'stupid':
qpms_scatsys_irrep_pack_matrix_stupid(&target_view[0][0], &fullmatrix_view[0][0],
self.s, iri)
else:
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):
def unpack_matrix(self, packedmatrix, iri, version='normal'):
"""Unpacks an "irrep-packed" excitation coefficient vector to full coordinates.
Parameters
@ -770,7 +768,11 @@ cdef class ScatteringSystem:
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],
if version == 'stupid':
qpms_scatsys_irrep_unpack_matrix_stupid(&target_view[0][0], &packedmatrix_view[0][0],
self.s, iri, 0)
else:
qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
self.s, iri, 0)
return target_np
@ -822,6 +824,16 @@ cdef class ScatteringSystem:
qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c, eta)
return target
def irrep_transform_matrix(self, qpms_iri_t iri):
self.check_s()
cdef size_t rlen = self.saecv_sizes[iri]
cdef size_t fullen = self.fecv_size
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
(rlen, fullen), dtype=complex, order='C')
cdef cdouble[:,::1] target_view = target
qpms_scatsys_irrep_transform_matrix(&target_view[0][0], self.s, iri)
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]
@ -965,8 +977,8 @@ cdef class ScatteringSystem:
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]):
with nogil:
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
@ -1006,18 +1018,13 @@ def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double max
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):
@ -1026,7 +1033,8 @@ cdef class _ScatteringSystemAtOmegaK:
self.sswk.eta = eta
@boundscheck(False)
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS,
ewaldparts = EwaldPart.FULL):
"""Evaluate electric field for a given excitation coefficient vector (periodic system)
Parameters
@ -1045,6 +1053,7 @@ cdef class _ScatteringSystemAtOmegaK:
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 qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
@ -1055,18 +1064,19 @@ cdef class _ScatteringSystemAtOmegaK:
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]):
with nogil:
for i in range(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)
res = qpms_scatsyswk_scattered_E_e(&self.sswk, btyp_c, &scv_view[0], pos, ewaldparts_c)
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):
@boundscheck(False)
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
# TODO examples
"""Evaluate scattered field "basis" (periodic system)
@ -1088,6 +1098,7 @@ cdef class _ScatteringSystemAtOmegaK:
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 qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
cdef Py_ssize_t fecv_size = self.fecv_size
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
@ -1095,20 +1106,20 @@ cdef class _ScatteringSystemAtOmegaK:
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]):
with nogil:
res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
for i in range(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)
qpms_scatsyswk_scattered_field_basis_e(res, &self.sswk, btyp_c, pos, ewaldparts_c)
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)
free(res)
return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
property fecv_size:
@ -1120,9 +1131,6 @@ cdef class _ScatteringSystemAtOmega:
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")
@ -1151,9 +1159,6 @@ cdef class _ScatteringSystemAtOmega:
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
@ -1237,7 +1242,6 @@ cdef class _ScatteringSystemAtOmega:
@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
@ -1274,19 +1278,85 @@ cdef class _ScatteringSystemAtOmega:
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]):
for i in range(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)
@boundscheck(False)
def scattered_field_basis(self, evalpos, blochvector=None, particle=None, btyp=QPMS_HANKEL_PLUS):
# TODO examples
# FIXME periodic case not implemented
"""Evaluate scattered field "basis"
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.
blochvector: array_like or None
Bloch vector, must be supplied (non-None) for periodic systems, else None.
particle_index: int or None (default), optional
A valid particle index; if specified, only the contribution of the given particle
is evaluated. Not applicable for periodic arrays.
btyp: BesselType, optional
Kind of the waves. Defaults to BesselType.HANKEL_PLUS.
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (n, 3)`
"Basis" fields at the positions given in `evalpos`, in cartesian coordinates.
`n` here stays either for `self.fecv_size` (if `particle==None`)
or `len(self.bspec_pi(particle))`.
"""
#if(btyp != QPMS_HANKEL_PLUS): # #TODO, IIRC not supported only for periodic systems
# raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef Py_ssize_t basissize = self.fecv_size if particle is None else len(self.bspec_pi(particle))
cdef qpms_ss_pi_t pi = particle
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], basissize, 3), dtype=complex)
cdef ccart3_t *res
cdef cart3_t pos
cdef Py_ssize_t i, j
with nogil:
res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t)) # thread-local
for i in range(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)
if particle is None:
qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, 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)
qpms_scatsysw_scattered_field_basis_pi(res, self.ssw, pi, btyp_c, pos)
for j in range(basissize):
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] + (basissize, 3))
def bspec_pi(self, pi):
return self.ss_pyref.bspec_pi(pi)
property bspecs:
def __get__(self):
return self.ss_pyref.bspecs
cdef class ScatteringMatrix:

View File

@ -150,6 +150,11 @@ cdef extern from "qpms_types.h":
cdouble mu
ctypedef enum qpms_coord_system_t:
pass
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
# maybe more if needed
cdef extern from "qpms_error.h":
@ -188,10 +193,13 @@ cdef extern from "vswf.h":
csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t n, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm)
cdef extern from "indexing.h":
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax)
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n)
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n)
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u)
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax) nogil
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) nogil
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n) nogil
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) nogil
qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax) nogil
void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n)
qpms_l_t qpms_nelem2lMax(qpms_y_t nelem)
# maybe more if needed
# Point generators from lattices.h
@ -254,95 +262,95 @@ cdef extern from "lattices.h":
cdef extern from "vectors.h":
cart2_t cart2_add(const cart2_t a, const cart2_t b)
cart2_t cart2_substract(const cart2_t a, const cart2_t b)
cart2_t cart2_scale(const double c, const cart2_t v)
double cart2_dot(const cart2_t a, const cart2_t b)
double cart2_normsq(const cart2_t a)
double cart2norm(const cart2_t v)
pol_t cart2pol(const cart2_t cart)
sph_t pol2sph_equator(const pol_t pol)
sph_t cart22sph(const cart2_t cart)
sph_t cart12sph_zaxis(double z)
cart3_t cart12cart3z(double z)
cart3_t cart22cart3xy(const cart2_t a)
cart2_t cart3xy2cart2(const cart3_t a)
double cart3_dot(const cart3_t a, const cart3_t b)
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b)
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c)
double cart3_normsq(const cart3_t a)
double cart3norm(const cart3_t v)
sph_t cart2sph(const cart3_t cart)
cart3_t sph2cart(const sph_t sph)
cart2_t pol2cart(const pol_t pol)
cart3_t pol2cart3_equator(const pol_t pol)
cart3_t cart3_add(const cart3_t a, const cart3_t b)
cart3_t cart3_substract(const cart3_t a, const cart3_t b)
cart3_t cart3_scale(const double c, const cart3_t v)
cart3_t cart3_divscale( const cart3_t v, const double c)
double cart3_dist(const cart3_t a, const cart3_t b)
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol)
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v)
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b)
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b)
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b)
ccart3_t cart32ccart3(cart3_t c)
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b)
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
csphvec_t csphvec_scale(cdouble c, const csphvec_t v)
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b)
sph_t sph_scale(double c, const sph_t s)
csph_t sph_cscale(cdouble c, const sph_t s)
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at)
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at)
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at)
csph_t sph2csph(sph_t s)
sph_t csph2sph(csph_t s)
csph_t ccart2csph(const ccart3_t cart)
csph_t cart2csph(const cart3_t cart)
ccart3_t csph2ccart(const csph_t sph)
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation)
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation)
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input)
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input)
double csphvec_norm(const csphvec_t a)
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance)
double csphvec_reldiff(const csphvec_t a, const csphvec_t b)
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t)
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t)
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t)
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t)
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t)
cart2_t cart2_add(const cart2_t a, const cart2_t b) nogil
cart2_t cart2_substract(const cart2_t a, const cart2_t b) nogil
cart2_t cart2_scale(const double c, const cart2_t v) nogil
double cart2_dot(const cart2_t a, const cart2_t b) nogil
double cart2_normsq(const cart2_t a) nogil
double cart2norm(const cart2_t v) nogil
pol_t cart2pol(const cart2_t cart) nogil
sph_t pol2sph_equator(const pol_t pol) nogil
sph_t cart22sph(const cart2_t cart) nogil
sph_t cart12sph_zaxis(double z) nogil
cart3_t cart12cart3z(double z) nogil
cart3_t cart22cart3xy(const cart2_t a) nogil
cart2_t cart3xy2cart2(const cart3_t a) nogil
double cart3_dot(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b) nogil
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c) nogil
double cart3_normsq(const cart3_t a) nogil
double cart3norm(const cart3_t v) nogil
sph_t cart2sph(const cart3_t cart) nogil
cart3_t sph2cart(const sph_t sph) nogil
cart2_t pol2cart(const pol_t pol) nogil
cart3_t pol2cart3_equator(const pol_t pol) nogil
cart3_t cart3_add(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_substract(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_scale(const double c, const cart3_t v) nogil
cart3_t cart3_divscale( const cart3_t v, const double c) nogil
double cart3_dist(const cart3_t a, const cart3_t b) nogil
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) nogil
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v) nogil
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) nogil
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) nogil
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b) nogil
ccart3_t cart32ccart3(cart3_t c) nogil
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) nogil
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) nogil
csphvec_t csphvec_scale(cdouble c, const csphvec_t v) nogil
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b) nogil
sph_t sph_scale(double c, const sph_t s) nogil
csph_t sph_cscale(cdouble c, const sph_t s) nogil
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) nogil
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at) nogil
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at) nogil
csph_t sph2csph(sph_t s) nogil
sph_t csph2sph(csph_t s) nogil
csph_t ccart2csph(const ccart3_t cart) nogil
csph_t cart2csph(const cart3_t cart) nogil
ccart3_t csph2ccart(const csph_t sph) nogil
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation) nogil
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) nogil
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) nogil
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) nogil
double csphvec_norm(const csphvec_t a) nogil
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance) nogil
double csphvec_reldiff(const csphvec_t a, const csphvec_t b) nogil
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t) nogil
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t) nogil
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t) nogil
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t) nogil
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t) nogil
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) nogil
void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
void cart3_to_double_array(double *a, cart3_t b)
cart3_t cart3_from_double_array(const double *a)
void cart2_to_double_array(double *a, cart2_t b)
cart2_t cart2_from_double_array(const double *a)
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
void cart3_to_double_array(double *a, cart3_t b) nogil
cart3_t cart3_from_double_array(const double *a) nogil
void cart2_to_double_array(double *a, cart2_t b) nogil
cart2_t cart2_from_double_array(const double *a) nogil
const cart2_t CART2_ZERO
cdef extern from "quaternions.h":
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q)
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q)
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q)
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q)
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q)
qpms_quat_t qpms_quat_conj(qpms_quat_t q)
double qpms_quat_norm(qpms_quat_t q)
double qpms_quat_imnorm(qpms_quat_t q)
qpms_quat_t qpms_quat_normalise(qpms_quat_t q)
qpms_quat_t qpms_quat_log(qpms_quat_t q)
qpms_quat_t qpms_quat_exp(qpms_quat_t q)
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent)
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q) nogil
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_conj(qpms_quat_t q) nogil
double qpms_quat_norm(qpms_quat_t q) nogil
double qpms_quat_imnorm(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_normalise(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_log(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_exp(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent) nogil
cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
qpms_m_t mp, qpms_m_t m)
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q)
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n)
qpms_quat_t qpms_quat_from_rotvector(cart3_t v)
qpms_m_t mp, qpms_m_t m) nogil
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q) nogil
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n) nogil
qpms_quat_t qpms_quat_from_rotvector(cart3_t v) nogil
cdef extern from "groups.h":
struct qpms_finite_group_irrep_t:
@ -490,6 +498,7 @@ cdef extern from "tmatrices.h":
qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, cdouble omega, const void *params)
struct qpms_tmatrix_generator_sphere_param_t:
qpms_epsmu_generator_t outside
qpms_epsmu_generator_t inside
@ -623,6 +632,7 @@ cdef extern from "scatsystem.h":
size_t *saecv_sizes
const qpms_finite_group_t *sym
qpms_scatsys_periodic_info_t per
qpms_trans_calculator *c
# per[] and other stuff not currently needed in cython
void qpms_scatsys_free(qpms_scatsys_t *s)
@ -638,6 +648,11 @@ cdef extern from "scatsystem.h":
cdouble omega, const qpms_tolerance_spec_t *tol)
qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega)
void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw)
cdouble *qpms_scatsys_irrep_transform_matrix(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_pack_matrix_stupid(cdouble *target_packed,
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_matrix_stupid(cdouble *target_full,
const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed,
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full,
@ -708,8 +723,20 @@ cdef extern from "scatsystem.h":
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where, qpms_ewald_part parts) nogil
qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss,
qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_t *ss,
qpms_ss_pi_t pi, qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where, qpms_ewald_part parts) nogil
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil
@ -718,12 +745,6 @@ cdef extern from "ewald.h":
cdouble val
double err
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
struct qpms_ewald3_constants_t:
qpms_l_t lMax
qpms_y_t nelem_sc

View File

@ -5,6 +5,9 @@
*/
#ifndef QPMS_ERROR_H
#define QPMS_ERROR_H
#ifdef __cplusplus
extern "C" {
#endif
#include "optim.h"
@ -266,4 +269,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
}\
}
#ifdef __cplusplus
}
#endif
#endif

View File

@ -3,6 +3,10 @@
*/
#ifndef QPMS_SPECFUNC_H
#define QPMS_SPECFUNC_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <gsl/gsl_sf_legendre.h>
@ -11,25 +15,25 @@
******************************************************************************/
// TODO unify types
qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array);
qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, _Complex double x, _Complex double *result_array);
typedef struct {
qpms_l_t lMax;
double *akn; // coefficients as in DLMF 10.49.1
//complex double *bkn; // coefficients of the derivatives
//_Complex double *bkn; // coefficients of the derivatives
} qpms_sbessel_calculator_t;
qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(void);
void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c);
qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax,
double x, complex double *result_array);
double x, _Complex double *result_array);
complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x);
_Complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, _Complex double x);
qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax,
complex double x, complex double *result_array);
_Complex double x, _Complex double *result_array);
/******************************************************************************
@ -107,4 +111,7 @@ double qpms_legendre0(int m, int n);
// Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2)
double qpms_legendred0(int m, int n);
#ifdef __cplusplus
}
#endif
#endif // QPMS_SPECFUNC_H

View File

@ -3,10 +3,13 @@
*/
#ifndef QPMS_TYPES_H
#define QPMS_TYPES_H
#include <complex.h>
#include <stdbool.h>
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
#include <stdint.h>
#include <stdbool.h>
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L)
@ -99,7 +102,9 @@ typedef enum {
* These bit flags are used by the functions declared in normalisation.h
* that return the appropriate convention-dependent factors.
*
* See @ref vswf_conventions for comparison of the various conventions used.
* \see @ref vswf_conventions for comparison of the various conventions used.
* \see @ref swf_conventions_qpms for description how the conventions are used internally and in the QPMS API.
*
*/
typedef enum {
QPMS_NORMALISATION_UNDEF = 0, ///< Convention undefined. This should not happen.
@ -189,7 +194,7 @@ typedef struct cart3_t {
/// 3D complex (actually 6D) coordinates. See also vectors.h.
typedef struct ccart3_t {
complex double x, y, z;
_Complex double x, y, z;
} ccart3_t;
/// 3D complex vector pair (represents the E, H fields).
@ -210,13 +215,13 @@ typedef struct sph_t {
/// Spherical coordinates with complex radial component. See also vectors.h.
typedef struct csph_t { // Do I really need this???
complex double r;
_Complex double r;
double theta, phi;
} csph_t;
/// 3D complex vector components in local spherical basis. See also vectors.h.
typedef struct csphvec_t {
complex double rc, thetac, phic;
_Complex double rc, thetac, phic;
} csphvec_t;
/// 2D polar coordinates. See also vectors.h.
@ -259,7 +264,7 @@ typedef enum {
* See quaternions.h for "methods".
*/
typedef struct qpms_quat_t {
complex double a, b;
_Complex double a, b;
} qpms_quat_t;
/// Quaternion type as four doubles.
@ -344,7 +349,7 @@ typedef struct qpms_tmatrix_t {
* there is also one with order \a m.
*/
const qpms_vswf_set_spec_t *spec;
complex double *m; ///< Matrix elements in row-major order.
_Complex double *m; ///< Matrix elements in row-major order.
bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
} qpms_tmatrix_t;
@ -420,11 +425,22 @@ typedef struct qpms_abstract_tmatrix_t {
/// A type holding electric permittivity and magnetic permeability of a material.
typedef struct qpms_epsmu_t {
complex double eps; ///< Relative permittivity.
complex double mu; ///< Relative permeability.
_Complex double eps; ///< Relative permittivity.
_Complex double mu; ///< Relative permeability.
} qpms_epsmu_t;
struct qpms_tolerance_spec_t; // See tolerances.h
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#ifdef __cplusplus
}
#endif
#endif // QPMS_TYPES

View File

@ -14,6 +14,10 @@
*/
#ifndef QPMSBLAS_H
#define QPMSBLAS_H
#ifdef __cplusplus
extern "C" {
#endif
#define QPMS_BLAS_INDEX_T long long int
#ifndef CBLAS_H
@ -31,4 +35,7 @@ void qpms_zgemm(CBLAS_LAYOUT Order, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE Tran
const _Complex double *B, const QPMS_BLAS_INDEX_T ldb,
const _Complex double *beta, _Complex double *C, const QPMS_BLAS_INDEX_T ldc);
#ifdef __cplusplus
}
#endif
#endif //QPMSBLAS_H

View File

@ -3,6 +3,9 @@
*/
#ifndef QPMS_WIGNER_H
#define QPMS_WIGNER_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include "vectors.h"
@ -93,7 +96,7 @@ static inline double qpms_quat_norm(const qpms_quat_t q) {
}
/// Test approximate equality of quaternions.
static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
static inline _Bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
}
@ -119,7 +122,7 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
}
/// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
static inline _Bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
return qpms_quat_norm(qpms_quat_sub(
qpms_quat_standardise(p, atol),
qpms_quat_standardise(q, atol))) <= atol;
@ -202,17 +205,17 @@ static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) {
* The D matrix are calculated using formulae (3), (4), (6), (7) from
* http://moble.github.io/spherical_functions/WignerDMatrices.html
*/
complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
_Complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
qpms_m_t mp, qpms_m_t m);
/// A VSWF representation element of the O(3) group.
/**
* TODO more doc.
*/
complex double qpms_vswf_irot_elem_from_irot3(
_Complex double qpms_vswf_irot_elem_from_irot3(
const qpms_irot3_t q, ///< The O(3) element in the quaternion representation.
qpms_l_t l, qpms_m_t mp, qpms_m_t m,
bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
_Bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
);
@ -251,7 +254,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
}
/// Test approximate equality of irot3.
static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
static inline _Bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
}
@ -303,4 +306,7 @@ static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
return qpms_irot3_zrot_angle(2 * M_PI * k / N);
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_WIGNER_H

View File

@ -24,6 +24,7 @@
#include "tolerances.h"
#include "beyn.h"
#include "tiny_inlines.h"
#include "lattices.h"
#ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
#include "qpmsblas.h"
@ -2049,6 +2050,61 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
cvf, where);
}
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
ccart3_t *target,
const qpms_scatsys_t *ss,
const qpms_ss_pi_t pi,
const qpms_bessel_t btyp,
const complex double k,
const cart3_t where
) {
qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
QPMS_UNTESTED;
const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates
QPMS_CRASHING_MALLOC(vswfs_sph, bspec->n * sizeof(*vswfs_sph));
const cart3_t particle_pos = ss->p[pi].pos;
const csph_t kr = sph_cscale(k, cart2sph(
cart3_substract(where, particle_pos)));
QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp));
for(size_t i = 0; i < bspec->n; ++i)
target[i] = csphvec2ccart_csph(vswfs_sph[i], kr);
free(vswfs_sph);
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, const qpms_ss_pi_t pi,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsys_scattered_field_basis_pi(target, ssw->ss, pi, btyp,
ssw->wavenumber, where);
}
qpms_errno_t qpms_scatsys_scattered_field_basis(
ccart3_t *target,
const qpms_scatsys_t *ss,
const qpms_bessel_t btyp,
const complex double k,
const cart3_t where
) {
qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
QPMS_UNTESTED;
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi)
QPMS_ENSURE_SUCCESS(qpms_scatsys_scattered_field_basis_pi(
target + ss->fecv_pstarts[pi], ss, pi, btyp, k, where));
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsysw_scattered_field_basis(
ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsys_scattered_field_basis(target, ssw->ss, btyp,
ssw->wavenumber, where);
}
const qpms_uvswfi_t ELDIPILIST[] = {6, 10, 14};
#define DIPSPECN 3 // We have three basis vectors
// Evaluates the regular electric dipole waves in the origin. The returned
@ -2058,11 +2114,15 @@ static inline const qpms_vswf_set_spec_t qpms_fill_regdipoles_0(
// bspec containing only electric dipoles
const qpms_vswf_set_spec_t dipspec = {
.n = DIPSPECN,
#if 0
.ilist = (qpms_uvswfi_t[]){
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
},
#else
.ilist = ELDIPILIST,
#endif
.lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
.capacity=0,
.norm = normalisation,
@ -2127,10 +2187,11 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
// For periodic lattices, we use directly the "alternative" implementation,
// using translation operator and regular dipole waves at zero
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp,
const complex double *cvf,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2162,13 +2223,13 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2182,11 +2243,17 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
return res;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) {
return qpms_scatsyswk_scattered_E_e(sswk, btyp, cvf, where, QPMS_EWALD_FULL);
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target,
const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2220,13 +2287,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2242,6 +2309,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsyswk_scattered_field_basis_e(
target, sswk, btyp, where, QPMS_EWALD_FULL);
}
#if 0
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, const complex double *cvr, cart3_t where) {

View File

@ -10,6 +10,10 @@
*/
#ifndef QPMS_SCATSYSTEM_H
#define QPMS_SCATSYSTEM_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include "vswf.h"
#include "tmatrices.h"
@ -107,7 +111,7 @@ typedef struct qpms_ss_orbit_type_t {
*
* TODO doc.
*/
complex double *irbases;
_Complex double *irbases;
/// TODO doc.
size_t instance_count;
/// Cumulative sum of the preceding ot->siza * ot->instance_count;
@ -243,9 +247,9 @@ typedef struct qpms_scatsys_at_omega_t {
* i.e in the order corresponding to \a ss->tm.
*/
qpms_tmatrix_t **tm;
complex double omega; ///< Angular frequency
_Complex double omega; ///< Angular frequency
qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
complex double wavenumber; ///< Background medium wave number
_Complex double wavenumber; ///< Background medium wave number
} qpms_scatsys_at_omega_t;
@ -291,7 +295,7 @@ typedef struct qpms_scatsys_at_omega_t {
*
*/
qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym,
complex double omega, const struct qpms_tolerance_spec_t *tol);
_Complex double omega, const struct qpms_tolerance_spec_t *tol);
/// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
void qpms_scatsys_free(qpms_scatsys_t *s);
@ -303,50 +307,50 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
/// Evaluates scattering system T-matrices at a given frequency.
/** Free the result using qpms_scatsys_at_omega_free() when done. */
qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
complex double omega);
_Complex double omega);
/// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis.
/** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
*
* TODO doc about shape etc.
*/
complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
_Complex double *qpms_scatsys_irrep_transform_matrix(_Complex double *target_U,
const qpms_scatsys_t *ss, qpms_iri_t iri);
/// Projects a "big" matrix onto an irrep (slow reference implementation).
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_matrix_stupid(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis (slow reference implementation).
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_matrix_stupid(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Projects a "big" matrix onto an irrep.
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_matrix(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis.
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_matrix(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Projects a "big" vector onto an irrep.
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_vector(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" vector into the full basis.
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_vector(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Global translation matrix.
@ -354,31 +358,31 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
* The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves).
* This may change in the future.
*/
complex double *qpms_scatsys_build_translation_matrix_full(
_Complex double *qpms_scatsys_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double k ///< Wave number to use in the translation matrix.
_Complex double k ///< Wave number to use in the translation matrix.
);
/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_full(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw
);
/// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
/** Might be useful for evaluation of cross sections and testing.
*/
complex double *qpms_scatsys_build_translation_matrix_e_full(
_Complex double *qpms_scatsys_build_translation_matrix_e_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double k, ///< Wave number to use in the translation matrix.
_Complex double k, ///< Wave number to use in the translation matrix.
qpms_bessel_t J
);
@ -387,12 +391,12 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
* The diagonal (particle self-) blocks are currently filled with zeros.
* This may change in the future.
*/
complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
_Complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
qpms_iri_t iri,
complex double k, ///< Wave number to use in the translation matrix.
_Complex double k, ///< Wave number to use in the translation matrix.
qpms_bessel_t J
);
@ -400,9 +404,9 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -410,9 +414,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -420,9 +424,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -435,7 +439,7 @@ typedef struct qpms_ss_LU {
bool full; ///< true if full matrix; false if irrep-packed.
qpms_iri_t iri; ///< Irrep index if `full == false`.
/// LU decomposition array.
complex double *a;
_Complex double *a;
/// Pivot index array, size at least max(1,min(m, n)).
int *ipiv;
} qpms_ss_LU;
@ -443,14 +447,14 @@ void qpms_ss_LU_free(qpms_ss_LU);
/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw
);
/// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri
@ -458,7 +462,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
/// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
_Complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems.
const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
@ -466,16 +470,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
/// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
_Complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri
);
/// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$.
complex double *qpms_scatsys_scatter_solve(
complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
const complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
_Complex double *qpms_scatsys_scatter_solve(
_Complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
const _Complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data.
);
@ -495,33 +499,33 @@ typedef struct qpms_scatsys_at_omega_k_t {
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
_Complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_k_t *sswk
);
/// Global translation matrix.
complex double *qpms_scatsys_periodic_build_translation_matrix_full(
_Complex double *qpms_scatsys_periodic_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double wavenumber, ///< Wave number to use in the translation matrix.
_Complex double wavenumber, ///< Wave number to use in the translation matrix.
const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates.
double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
);
/// Global translation matrix.
complex double *qpms_scatsyswk_build_translation_matrix_full(
_Complex double *qpms_scatsyswk_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_k_t *sswk
);
/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only.
qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_k_t *sswk
);
@ -539,7 +543,7 @@ struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
const qpms_scatsys_t *ss,
/// Wavevector in cartesian coordinates (must lie in the lattice plane).
const double k[3],
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -561,12 +565,12 @@ struct qpms_finite_group_t;
/// Constructs a "full matrix action" of a point group element for an orbit type.
/** TODO detailed doc */
complex double *qpms_orbit_action_matrix(
_Complex double *qpms_orbit_action_matrix(
/// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -579,12 +583,12 @@ complex double *qpms_orbit_action_matrix(
/// Constructs a dense matrix representation of a irrep projector for an orbit type.
/** TODO detailed doc */
complex double *qpms_orbit_irrep_projector_matrix(
_Complex double *qpms_orbit_irrep_projector_matrix(
/// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -596,14 +600,14 @@ complex double *qpms_orbit_irrep_projector_matrix(
const qpms_iri_t iri);
/// TODO DOC!!!!!
complex double *qpms_orbit_irrep_basis(
_Complex double *qpms_orbit_irrep_basis(
/// Here theh size of theh basis shall be saved,
size_t *basis_size,
/// Target array. If NULL, a new one is allocated.
/** The size of the array is basis_size * (orbit->size * bspec->n)
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -618,10 +622,10 @@ complex double *qpms_orbit_irrep_basis(
/// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc!
* \returns target_full if target_full was not NULL, otherwise the newly allocated array. */
complex double *qpms_scatsys_incident_field_vector_full(
_Complex double *qpms_scatsys_incident_field_vector_full(
/// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */
complex double *target_full,
_Complex double *target_full,
const qpms_scatsys_t *ss,
qpms_incfield_t field_at_point,
const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
@ -629,9 +633,9 @@ complex double *qpms_scatsys_incident_field_vector_full(
);
/// Applies T-matrices onto an incident field vector in the full basis.
complex double *qpms_scatsysw_apply_Tmatrices_full(
complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
_Complex double *qpms_scatsysw_apply_Tmatrices_full(
_Complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
const _Complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
const qpms_scatsys_at_omega_t *ssw
);
@ -650,7 +654,7 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
const qpms_scatsys_t *ss,
/// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
qpms_iri_t iri,
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -673,7 +677,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
const qpms_scatsys_t *ss,
double eta, ///< Ewald sum parameter
const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system.
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -687,10 +691,10 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
#if 0
/// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc! */
complex double *qpms_scatsys_incident_field_vector_irrep_packed(
_Complex double *qpms_scatsys_incident_field_vector_irrep_packed(
/// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */
complex double *target_full,
_Complex double *target_full,
const qpms_scatsys_t *ss,
const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym.
qpms_incfield_t field_at_point,
@ -706,7 +710,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsysw_scattered_E()
* \see qpms_scatsysw_scattered_E(), qpms_scatsys_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_E() for periodic systems.
*
@ -715,8 +719,8 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
ccart3_t qpms_scatsys_scattered_E(
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
complex double wavenumber, ///< Wavenumber of the background medium.
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
_Complex double wavenumber, ///< Wavenumber of the background medium.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -727,14 +731,96 @@ ccart3_t qpms_scatsys_scattered_E(
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsys_scattered_E()
* \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_E() for periodic systems.
*/
ccart3_t qpms_scatsysw_scattered_E(
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field at a given point.
/**
* This function evaluates all the included VSWFs from the particles in the system, evaluated
* at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
* would be equivalent to the result of qpms_scatsys_scattered_E().
*
* \see qpms_scatsysw_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
*
*/
qpms_errno_t qpms_scatsys_scattered_field_basis(
ccart3_t *target, ///< Target array of length \a ss->fecv_size
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
_Complex double wavenumber, ///< Wavenumber of the background medium.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field of a given particle at a given point.
/**
* This function evaluates all the included VSWFs from a particle in the system, evaluated
* at a given point.
*
* Finite systems only.
*
* \see qpms_scatsys_scattered_field_basis() for the whole-array equivalent.
*
* \see qpms_scatsysw_scattered_field_basis_pi()
*/
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
const qpms_scatsys_t *ss,
qpms_ss_pi_t pi, ///< Particle index
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
_Complex double wavenumber, ///< Wavenumber of the background medium
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field at a given point.
/**
* This function evaluates all the included VSWFs from the particles in the system, evaluated
* at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
* would be equivalent to the result of qpms_scatsysw_scattered_E().
*
* Note that this might require a relatively large amount of memory. Depending
* on the application, qpms_scatsysw_scattered_field_basis_pi() might be a better
* alternative.
*
* \see qpms_scatsysw_scattered_field_basis_pi() for the single-particle equivalent.
*
* \see qpms_scatsys_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
*
*/
qpms_errno_t qpms_scatsysw_scattered_field_basis(
ccart3_t *target, ///< Target array of length \a ss->fecv_size
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field of a given particle at a given point.
/**
* This function evaluates all the included VSWFs from a particle in the system, evaluated
* at a given point.
*
* Finite systems only.
*
* \see qpms_scatsysw_scattered_field_basis() for the whole-array equivalent.
*
* \see qpms_scatsys_scattered_field_basis_pi()
*/
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
const qpms_scatsys_at_omega_t *ssw,
qpms_ss_pi_t pi, ///< Particle index
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -750,8 +836,8 @@ ccart3_t qpms_scatsysw_scattered_E(
ccart3_t qpms_scatsys_scattered_E__alt(
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
complex double wavenumber, ///< Wavenumber of the background medium.
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
_Complex double wavenumber, ///< Wavenumber of the background medium.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -767,7 +853,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(
ccart3_t qpms_scatsysw_scattered_E__alt(
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -786,10 +872,18 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
ccart3_t qpms_scatsyswk_scattered_E(
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
ccart3_t qpms_scatsyswk_scattered_E_e(
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
qpms_ewald_part parts
);
/// Evaluates "scattered" field basis functions in a periodic system.
/**
*
@ -802,10 +896,18 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis is evaluated.
);
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
);
/// Adjusted Ewadl parameter to avoid high-frequency breakdown.
// TODO DOC
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, _Complex double wavenumber, const double wavevector[3]);
#if 0
/** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
@ -815,10 +917,13 @@ double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber,
*/
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, ///< Irreducible representation
const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
const _Complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
cart3_t where, ///< Evaluation point.
complex double k ///< Wave number.
_Complex double k ///< Wave number.
);
#endif
#ifdef __cplusplus
}
#endif
#endif //QPMS_SCATSYSTEM_H

48
qpms/scatsystem_dbg.c Normal file
View File

@ -0,0 +1,48 @@
#include "scatsystem_dbg.h"
#include "translations_dbg.h"
#include "indexing.h"
// analogue of qpms_scatsyswk_scattered_field_basis_e()
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
complex double *target, ///< Target array of size (2 * sswk->ssw->ss->c->lMax + 1) * sswk->ssw->ss->p_count
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t where, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
QPMS_NOT_IMPLEMENTED("Only scattered field with first kind Hankel functions currently implemented.");
const qpms_scatsys_t *ss = sswk->ssw->ss;
if (ss->lattice_dimension != 2)
QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented");
//ccart3_t res = {0,0,0};
//ccart3_t res_kc = {0,0,0}; // kahan sum compensation
const size_t particle_nelem = qpms_lMax2nelem_sc(2*ss->c->lMax + 1);
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
const cart3_t particle_pos = ss->p[pi].pos;
const cart3_t origin_cart = {0, 0, 0};
QPMS_ASSERT(sswk->k[2] == 0); // At least not implemented now
QPMS_ASSERT(ss->per.lattice_basis[0].z == 0);
QPMS_ASSERT(ss->per.lattice_basis[1].z == 0);
// Same choices as in qpms_ss_ppair_W32xy; TODO make it more dynamic
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_test_sswf_e32(ss->c,
target + pi * particle_nelem, NULL,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK, parts));
}
return QPMS_SUCCESS;
}

10
qpms/scatsystem_dbg.h Normal file
View File

@ -0,0 +1,10 @@
#include "scatsystem.h"
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
_Complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
);

View File

@ -21,40 +21,44 @@
*/
#ifndef SYMMETRIES_H
#define SYMMETRIES_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <cblas.h>
/// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring).
complex double *qpms_zflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of the y coordinate sign flip operation (xz-plane mirroring).
complex double *qpms_yflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_yflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of the x coordinate sign flip operation (yz-plane mirroring).
complex double *qpms_xflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_xflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of a rotation around the z-axis
complex double *qpms_zrot_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zrot_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
double phi ///< Rotation angle
);
/// Dense matrix representation of a "rational" rotation around the z-axis
/** Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
*/
complex double *qpms_zrot_rational_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zrot_rational_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
int N,
int w
);
/// Dense matrix (uvswi-indexed) representation of any O(3) transformation.
complex double *qpms_irot3_uvswfi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_irot3_uvswfi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
const qpms_irot3_t transf);
@ -74,5 +78,9 @@ size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol);
* Works on real and imaginary parts separately.
* TODO doc.
*/
size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol);
size_t qpms_czero_roundoff_clean(_Complex double *arr, size_t nmemb, double atol);
#ifdef __cplusplus
}
#endif
#endif // SYMMETRIES_H

View File

@ -4,6 +4,7 @@
#ifndef TINY_INLINES_H
#define TINY_INLINES_H
#include <stdlib.h>
#include <complex.h>
static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; }
@ -19,8 +20,8 @@ static inline int min1pow_m_neg(int m) {
#if 0
#ifdef __GSL_SF_LEGENDRE_H__
static inline complex double
spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, complex double exp_imf) {
static inline _Complex double
spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, _Complex double exp_imf) {
return;
}
@ -28,9 +29,9 @@ spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m,
#endif
// this has shitty precision:
// static inline complex double ipow(int x) { return cpow(I, x); }
// static inline _Complex double ipow(int x) { return cpow(I, x); }
static inline complex double ipow(int x) {
static inline _Complex double ipow(int x) {
x = ((x % 4) + 4) % 4;
switch(x) {
case 0:

View File

@ -882,7 +882,9 @@ qpms_errno_t qpms_tmatrix_axialsym_fill(
// R now contains -T^T.
for(size_t i1 = 0; i1 < bspec->n; ++i1)
for(size_t i2 = 0; i2 < bspec->n; ++i2) {
if (reindex[i1] == ~(size_t) 0 && reindex[i2] == ~(size_t) 0) QPMS_WTF;
if (reindex[i1] == QPMS_VSWF_SET_REINDEX_UNMAPPED &&
reindex[i2] == QPMS_VSWF_SET_REINDEX_UNMAPPED)
QPMS_WTF;
const size_t it = i1 * bspec->n + i2;
const size_t iQR = reindex[i1] + reindex[i2] * bspecQR->n;
t->m[it] = -R[iQR];
@ -988,6 +990,12 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
return QPMS_SUCCESS;
}
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
complex double omega_ignored, const void *params_ignored) {
memset(t->m, 0, SQ(t->spec->n)*sizeof(*(t->m)));
return QPMS_SUCCESS;
}
void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
switch(f->typ) {
case QPMS_TMATRIX_OPERATION_NOOP:

View File

@ -3,6 +3,10 @@
*/
#ifndef TMATRICES_H
#define TMATRICES_H
#ifdef __cplusplus
extern "C" {
#endif
// #include "qpms_types.h" // included via materials.h
// #include <gsl/gsl_spline.h> // included via materials.h
#include "materials.h"
@ -14,7 +18,7 @@ struct qpms_finite_group_t;
typedef struct qpms_finite_group_t qpms_finite_group_t;
/// Returns a pointer to the beginning of the T-matrix row \a rowno.
static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
static inline _Complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
return t->m + (t->spec->n * rowno);
}
@ -42,7 +46,7 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t);
* This function actually checks for identical vswf specs.
* TODO define constants with "default" atol, rtol for this function.
*/
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
_Bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
const double rtol, const double atol);
/// Creates a T-matrix from another matrix and a symmetry operation.
@ -53,7 +57,7 @@ bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
*/
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
const qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
@ -64,7 +68,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
*/
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Symmetrizes a T-matrix with an involution symmetry operation.
@ -75,7 +79,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
*/
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
const qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
@ -105,7 +109,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
*/
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
@ -149,7 +153,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
double **freqs_su, ///< Frequencies in SCUFF units (optional).
/// The resulting T-matrices (optional).
qpms_tmatrix_t **tmatrices_array,
complex double **tmdata ///< The T-matrices raw contents
_Complex double **tmdata ///< The T-matrices raw contents
);
/// Tells whether qpms_load_scuff_tmatrix should crash if fopen() fails.
@ -161,7 +165,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
* This is desirable e.g. when used in Python (so that proper exception can
* be thrown).
*/
extern bool qpms_load_scuff_tmatrix_crash_on_failure;
extern _Bool qpms_load_scuff_tmatrix_crash_on_failure;
/// Loads a scuff-tmatrix generated file.
/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a
@ -189,7 +193,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* is accessed via
* (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
*/
complex double ** tmdata
_Complex double ** tmdata
);
/// In-place application of point group elements on raw T-matrix data.
@ -200,7 +204,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* TODO more doc.
*/
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
complex double *tmdata, const size_t tmcount,
_Complex double *tmdata, const size_t tmcount,
const qpms_vswf_set_spec_t *bspec,
size_t n_symops,
const qpms_irot3_t *symops
@ -213,7 +217,7 @@ qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
* TODO more doc.
*/
qpms_errno_t qpms_symmetrise_tmdata_finite_group(
complex double *tmdata, const size_t tmcount,
_Complex double *tmdata, const size_t tmcount,
const qpms_vswf_set_spec_t *bspec,
const qpms_finite_group_t *pointgroup
);
@ -242,9 +246,9 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
);
/// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$.
complex double *qpms_apply_tmatrix(
complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
const complex double *a, ///< Incident field coefficient array of size T->spec->n.
_Complex double *qpms_apply_tmatrix(
_Complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
const _Complex double *a, ///< Incident field coefficient array of size T->spec->n.
const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
);
@ -254,10 +258,11 @@ complex double *qpms_apply_tmatrix(
* qpms_tmatrix_generator_interpolator()
* qpms_tmatrix_generator_sphere()
* qpms_tmatrix_generator_constant()
* qpms_tmatrix_generator_zero()
*/
typedef struct qpms_tmatrix_generator_t {
qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *params ///< Implementation dependent parameters.
);
const void *params; ///< Parameter pointer passed to the function.
@ -267,7 +272,7 @@ typedef struct qpms_tmatrix_generator_t {
qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
const qpms_vswf_set_spec_t *bspec,
qpms_tmatrix_generator_t gen,
complex double omega);
_Complex double omega);
/// Implementation of qpms_matrix_generator_t that just copies a constant matrix.
@ -275,11 +280,20 @@ qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
* the same base spec.
*/
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
complex double omega,
_Complex double omega_ignored,
/// Source T-matrix, real type is (const qpms_tmatrix_t*).
const void *tmatrix_orig
);
/// Dummy implementation of qpms_tmatrix_generator_t, yielding a zero matrix.
/**
* This effectively represents a part of the background medium without
* any scatterer present.
*/
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
_Complex double omega_ignored,
const void *params_ignored);
/* Fuck this, include the whole <gsl/gsl_spline.h>
typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
typedef struct gsl_interp_type gsl_interp_type;
@ -295,7 +309,7 @@ extern const gsl_interp_type * gsl_interp_steffen;
// struct gsl_interp_accel; // use if lookup proves to be too slow
typedef struct qpms_tmatrix_interpolator_t {
const qpms_vswf_set_spec_t *bspec;
//bool owns_bspec;
//_Bool owns_bspec;
gsl_spline **splines_real; ///< There will be a spline object for each nonzero element
gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
// gsl_interp_accel **accel_real;
@ -317,7 +331,7 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
const gsl_interp_type *iptype
//, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
//, _Bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
/*, ...? */);
@ -325,7 +339,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
/** As in qpms_tmatrix_interpolator_eval(), the imaginary part of frequency is discarded!
*/
qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *interpolator ///< Parameter of type qpms_tmatrix_interpolator_t *.
);
@ -342,14 +356,14 @@ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matri
*
* TODO better doc.
*/
complex double *qpms_mie_coefficients_reflection(
complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
_Complex double *qpms_mie_coefficients_reflection(
_Complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e, ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e, ///< Relative permeability of the surrounding medium.
qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR.
qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS.
);
@ -357,10 +371,10 @@ complex double *qpms_mie_coefficients_reflection(
/// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory.
qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e ///< Relative permeability of the surrounding medium.
);
/// Parameter structure for qpms_tmatrix_generator_sphere().
@ -372,7 +386,7 @@ typedef struct qpms_tmatrix_generator_sphere_param_t {
/// T-matrix generator for spherical particles using Lorentz-Mie solution.
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
complex double omega,
_Complex double omega,
const void *params ///< Of type qpms_tmatrix_generator_sphere_param_t.
);
@ -380,10 +394,10 @@ qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e ///< Relative permeability of the surrounding medium.
) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
@ -395,8 +409,8 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
double a, ///< Radius of the sphere.
double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium.
_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
_Complex double epsilon_bg ///< Relative permittivity of the background medium.
);
/// Convenience function to calculate T-matrix of a non-magnetic spherical particle using the permittivity values.
@ -404,8 +418,8 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere.
double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium.
_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
_Complex double epsilon_bg ///< Relative permittivity of the background medium.
) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
@ -460,7 +474,7 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta,
*/
qpms_errno_t qpms_tmatrix_axialsym_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
qpms_epsmu_t outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -474,7 +488,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill(
/// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry.
static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
const qpms_vswf_set_spec_t *bspec,
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
qpms_epsmu_t outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -501,22 +515,22 @@ typedef struct qpms_tmatrix_generator_axialsym_param_t {
/// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t
qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *params ///< Parameters of type qpms_tmatrix_generator_axialsym_param_t.
);
/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful for debugging).
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target,
complex double omega,
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(_Complex double *target,
_Complex double omega,
const qpms_tmatrix_generator_axialsym_param_t *param,
qpms_normalisation_t norm,
qpms_bessel_t J
);
/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging).
qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target,
complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(_Complex double *target,
_Complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm,
qpms_bessel_t J ///< Use QPMS_BESSEL_REGULAR to calculate \f$ R^T\f$ or QPMS_HANKEL_PLUS to calculate \f$ Q^T\f$.
);
@ -537,7 +551,7 @@ typedef struct qpms_tmatrix_function_t {
/// Convenience function to create a new T-matrix from qpms_tmatrix_function_t.
// FIXME the name is not very intuitive.
static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) {
static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, _Complex double omega) {
return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
}
@ -558,9 +572,9 @@ typedef enum {
struct qpms_tmatrix_operation_lrmatrix {
/// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */
complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
bool owns_m; ///< Whether \a m is owned by this;
_Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
_Bool owns_m; ///< Whether \a m is owned by this;
};
struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
@ -591,9 +605,9 @@ struct qpms_tmatrix_operation_compose_chain {
struct qpms_tmatrix_operation_scmulz {
/// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */
complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
bool owns_m; ///< Whether \a m is owned by this.
_Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
_Bool owns_m; ///< Whether \a m is owned by this.
};
/// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t.
@ -602,7 +616,7 @@ struct qpms_tmatrix_operation_scmulz {
struct qpms_tmatrix_operation_irot3arr {
size_t n; ///< Number of rotoreflections;
qpms_irot3_t *ops; ///< Rotoreflection array of size \a n.
bool owns_ops; ///< Whether \a ops array is owned by this.
_Bool owns_ops; ///< Whether \a ops array is owned by this.
};
/// A generic T-matrix transformation operator.
@ -689,4 +703,7 @@ typedef struct qpms_abstract_particle_t {
typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
#endif // 0
#ifdef __cplusplus
}
#endif
#endif //TMATRICES_H

View File

@ -1,6 +1,9 @@
/*! \file tolerances.h */
#ifndef QPMS_TOLERANCES_H
#define QPMS_TOLERANCES_H
#ifdef __cplusplus
extern "C" {
#endif
// TODO DOC
@ -12,4 +15,7 @@ typedef struct qpms_tolerance_spec_t {
/// A rather arbitrary default tolerance.
static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-9, .rtol = 1e-8};
#ifdef __cplusplus
}
#endif
#endif // QPMS_TOLERANCES_H

View File

@ -1,5 +1,6 @@
#include <math.h>
#include "qpms_types.h"
#include <complex.h>
#include "qpms_specfunc.h"
#include "gaunt.h"
#include "translations.h"
@ -15,6 +16,10 @@
#include "normalisation.h"
#include "translations_inlines.h"
#if defined LATTICESUMS31 || defined LATTICESUMS32
#include "lattices.h"
#endif
/*
* Define macros with additional factors that "should not be there" according
* to the "original" formulae but are needed to work with my vswfs.
@ -825,7 +830,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = Aerr || Berr;
const bool do_sigma0 = (particle_shift == 0)//DIFF21((particle_shift.x == 0) && (particle_shift.y == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@ -846,15 +851,15 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
QPMS_ENSURE_SUCCESS(ewald31z_sigma_short_points_and_shift(sigmas_short, serr_short, //DIFF21
c->e3c, eta, k, nRpoints, Rpoints, beta, particle_shift));
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = serr_short[y] + serr_long[y];
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0) {
QPMS_ENSURE_SUCCESS(ewald31z_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
@ -872,7 +877,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
// TODO skip if ... (N.B. skip will be different for 31z and 32)
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p = n + nu - 2*q;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p);
const complex double multiplier = c->A_multipliers[i][q];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Asum, &Asumc, multiplier * sigma);
@ -887,7 +892,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p_ = n + nu - 2*q + 1;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_);
const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Bsum, &Bsumc, multiplier * sigma);
@ -937,7 +942,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = Aerr || Berr;
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@ -1006,17 +1011,17 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
PGen_destroy(&Rgen);
}
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? sigmas_long[y] : 0);
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? serr_long[y] : 0);
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
@ -1034,7 +1039,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
// TODO skip if ...
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p = n + nu - 2*q;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p);
const complex double multiplier = c->A_multipliers[i][q];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Asum, &Asumc, multiplier * sigma);
@ -1049,7 +1054,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p_ = n + nu - 2*q + 1;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_);
const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Bsum, &Bsumc, multiplier * sigma);

View File

@ -26,9 +26,12 @@
*/
#ifndef QPMS_TRANSLATIONS_H
#define QPMS_TRANSLATIONS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "vectors.h"
#include "qpms_types.h"
#include <complex.h>
#include <stdbool.h>
#include <stddef.h>
@ -37,10 +40,10 @@
#endif
complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
_Complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
_Complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
/// Structure holding the constant factors in normalisation operators.
@ -65,8 +68,8 @@ typedef struct qpms_trans_calculator {
qpms_normalisation_t normalisation;
qpms_l_t lMax;
qpms_y_t nelem;
complex double **A_multipliers;
complex double **B_multipliers;
_Complex double **A_multipliers;
_Complex double **B_multipliers;
#if 0
// Normalised values of the Legendre functions and derivatives
// for θ == π/2, i.e. for the 2D case.
@ -90,44 +93,44 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation
/// Destructor for qpms_trans_calculator_t.
void qpms_trans_calculator_free(qpms_trans_calculator *);
complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
_Complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
_Complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
size_t deststride, size_t srcstride,
csph_t kdlj, bool r_ge_d, qpms_bessel_t J);
// TODO update the types later
complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double *Adest, _Complex double *Bdest,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
size_t deststride, size_t srcstride,
complex double kdlj_r, double kdlj_theta, double kdlj_phi,
_Complex double kdlj_r, double kdlj_theta, double kdlj_phi,
int r_ge_d, int J);
// Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c,
complex double *target,
_Complex double *target,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
@ -138,12 +141,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
/// and with automatic \a r_ge_d = `false`.
qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
const qpms_trans_calculator *c,
complex double *target,
_Complex double *target,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
complex double k, cart3_t destpos, cart3_t srcpos,
_Complex double k, cart3_t destpos, cart3_t srcpos,
qpms_bessel_t J
/// Workspace has to be at least 2 * c->neleme**2 long
);
@ -156,10 +159,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -167,10 +170,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
);
int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -180,12 +183,12 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
// Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
complex double *target, double *err,
_Complex double *target, double *err,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -193,12 +196,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
);
qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
complex double *target, double *err,
_Complex double *target, double *err,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -212,8 +215,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
// e31z means that the particles are positioned along the z-axis;
// their positions and K-values are then denoted by a single z-coordinate
int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const double k,
const double unitcell_area, // just lattice period
@ -244,4 +247,7 @@ int qpms_cython_trans_calculator_get_AB_arrays_loop(
#endif //QPMS_COMPILE_PYTHON_EXTENSIONS
#ifdef __cplusplus
}
#endif
#endif // QPMS_TRANSLATIONS_H

149
qpms/translations_dbg.c Normal file
View File

@ -0,0 +1,149 @@
#include "qpms_types.h"
#include "qpms_specfunc.h"
#include "translations_dbg.h"
#include "indexing.h"
#include "tiny_inlines.h"
#include "qpms_error.h"
#include "normalisation.h"
#include "lattices.h"
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
complex double *dest, const csph_t kdlj, const qpms_bessel_t J) {
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l)
for (qpms_m_t m = 0; m <= l; ++m) {
const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
dest[y_sc] = NAN+I*NAN;
}
// TODO warn? different return value?
return 0;
}
double legendre_buf[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)];
complex double bessel_buf[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
double costheta = cos(kdlj.theta);
QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
costheta,-1,legendre_buf));
QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf));
for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l)
for (qpms_m_t m = -l; m <= l; ++m) {
const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
dest[y_sc] = bessel_buf[l]
* legendre_buf[gsl_sf_legendre_array_index(l, abs(m))]
* cexp(I * m * kdlj.phi);
}
return 0;
}
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J) {
csph_t kdlj = cart2csph(r);
kdlj.r *= k;
return qpms_trans_calculator_test_sswf(c, dest, kdlj, J);
}
#ifdef LATTICESUMS32
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
complex double * const sigmas_total, double * serr_total,
const double eta, const complex double k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = !!serr_total;
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
complex double *sigmas_short = malloc(sizeof(complex double)*nelem2_sc);
complex double *sigmas_long = malloc(sizeof(complex double)*nelem2_sc);
//complex double *sigmas_total = malloc(sizeof(complex double)*nelem2_sc);
double *serr_short, *serr_long;//, *serr_total;
if(doerr) {
serr_short = malloc(sizeof(double)*nelem2_sc);
serr_long = malloc(sizeof(double)*nelem2_sc);
//serr_total = malloc(sizeof(double)*nelem2_sc);
} else serr_short = serr_long = serr_total = NULL;
const double unitcell_area = l2d_unitcell_area(b1, b2);
cart2_t rb1, rb2; // reciprocal basis
QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2));
if (parts & QPMS_EWALD_LONG_RANGE) {
PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL,
#ifdef GEN_KSHIFTEDPOINTS
beta,
#else
CART2_ZERO,
#endif
0, true, maxK, false);
QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k,
unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
#ifdef GEN_KSHIFTEDPOINTS
true,
#else
false,
#endif
cart22cart3xy(beta), particle_shift));
if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Kgen);
}
if (parts & QPMS_EWALD_SHORT_RANGE) {
PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL,
#ifdef GEN_RSHIFTEDPOINTS
cart2_scale(-1 /*CHECKSIGN*/, cart3xy2cart2(particle_shift)),
#else
CART2_ZERO,
#endif
0, !do_sigma0, maxR, false);
#ifdef GEN_RSHIFTEDPOINTS // rather ugly hacks, LPTODO cleanup
if (particle_shift.z != 0) {
const cart3_t zshift = {0, 0, -particle_shift.z /*CHECKSIGN*/};
Rgen = Pgen_shifted_new(Rgen, zshift);
}
#endif
QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
particle_shift.z ? LAT_2D_IN_3D : LAT_2D_IN_3D_XYONLY, &Rgen,
#ifdef GEN_RSHIFTEDPOINTS
true,
#else
false,
#endif
cart22cart3xy(beta), particle_shift));
if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Rgen);
}
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? sigmas_long[y] : 0);
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? serr_long[y] : 0);
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
free(sigmas_short);
free(sigmas_long);
//free(sigmas_total);
if(doerr) {
free(serr_short);
free(serr_long);
//free(serr_total);
}
return 0;
}
#endif // LATTICESUMS32

52
qpms/translations_dbg.h Normal file
View File

@ -0,0 +1,52 @@
/*! \file translations_dbg.h
* \brief Auxiliary functions for debugging and testing (complementary to translations.h)
*
* The following functions evaluate the SSWFs the same way as
* certain functions in translations.h, but they leave them
* as they are, without rearranging them into the translation
* operators.
*/
#ifndef TRANSLATIONS_DBG_H
#include "translations.h"
/// Evaluates scalar spherical wavefunctions at a point (spherical coordinates)
/**
* \f[
* \phi_l^m (kr, \theta, \varphi) = h_l(kr) \rawLeg{l}{m}(\cos \theta) e^{im\varphi}
* \f]
*
* The wavefunctions are evaluated up to the degree `2 * c->lMax + 1`.
*
* Currently, THIS ALWAYS USES NON-NORMALISED LEGENDRE FUNCTIONS (TODO FIXME?).
*/
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
_Complex double *dest, ///< Result array, length `qpms_lMax2nelem_sc(c->lMax)`
const csph_t kdlj, ///< Evaluation point (incl. wavenumber)
const qpms_bessel_t J);
/// Evaluates scalar spherical wavefunctions at a point (cartesian coordinates)
/**
* Currently, THIS ALWAYS USES NON-NORMALISED LEGENDRE FUNCTIONS (TODO FIXME?).
*
* \see qpms_trans_calculator_test_sswf()
*/
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
_Complex double *dest, ///< Result array, length `qpms_lMax2nelem_sc(c->lMax)`
const _Complex double k, ///< Wavenumber
const cart3_t r, ///< Evaluation point (excl. wavenumber)
const qpms_bessel_t J);
#ifdef LATTICESUMS32
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
_Complex double * const sigmas_total, double * serr_total,
const double eta, const _Complex double k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts);
#endif // LATTICESUMS32
#endif // TRANSLATIONS_DBG_H

View File

@ -6,12 +6,12 @@
/// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix.
// TODO DOC
static inline void qpms_trans_array_from_AB(
complex double *t,
_Complex double *t,
const qpms_vswf_set_spec_t *const t_destspec,
const size_t t_deststride,
const qpms_vswf_set_spec_t *const t_srcspec,
const size_t t_srcstride,
const complex double *const A, const complex double *const B,
const _Complex double *const A, const _Complex double *const B,
/// A and B matrices' lMax.
/** This also determines their size and stride: they are assumed to
* be square matrices of size `nelem * nelem` where

View File

@ -1,4 +1,5 @@
#include "groups.h"
#include <complex.h>
/// Trivial group, with one (reduntant) generator.
/**

View File

@ -4,6 +4,7 @@
#ifndef VECTORS_H
#define VECTORS_H
#include <math.h>
#include <complex.h>
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
#endif
@ -198,12 +199,12 @@ static inline double cart3_dist(const cart3_t a, const cart3_t b) {
return cart3norm(cart3_substract(a,b));
}
static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
static inline _Bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
}
/// Complex 3D vector scaling.
static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) {
static inline ccart3_t ccart3_scale(const _Complex double c, const ccart3_t v) {
ccart3_t res = {c * v.x, c * v.y, c * v.z};
return res;
}
@ -221,7 +222,7 @@ static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
}
/// Complex 3D cartesian vector "dot product" without conjugation.
static inline complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
static inline _Complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
@ -244,13 +245,13 @@ static inline csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
}
/// Complex 3D vector (geographic coordinates) scaling.
static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
static inline csphvec_t csphvec_scale(_Complex double c, const csphvec_t v) {
csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
return res;
}
/// Complex 3D vector (geographic coordinates) "dot product" without conjugation.
static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
static inline _Complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
//N.B. no complex conjugation done here
return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
}
@ -262,7 +263,7 @@ static inline sph_t sph_scale(double c, const sph_t s) {
}
/// "Complex spherical" coordinate system scaling.
static inline csph_t sph_cscale(complex double c, const sph_t s) {
static inline csph_t sph_cscale(_Complex double c, const sph_t s) {
csph_t res = {c * s.r, s.theta, s.phi};
return res;
}
@ -358,8 +359,9 @@ static inline csph_t ccart2csph(const ccart3_t cart) {
/// Real 3D cartesian to spherical (complex r) coordinates conversion. See @ref coord_conversions.
static inline csph_t cart2csph(const cart3_t cart) {
csph_t sph;
sph.r = cart3norm(cart);
sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
const double r = cart3norm(cart);
sph.r = r;
sph.theta = sph.r ? acos(cart.z / r) : M_PI_2;
sph.phi = atan2(cart.y, cart.x);
return sph;
}
@ -578,7 +580,6 @@ static inline double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) {
" to 1D not allowed.");
}
/// Coordinate conversion of point arrays (something to something).
/** The dest and src arrays must not overlap */
static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
@ -589,31 +590,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
sph_t *d = (sph_t *) dest;
switch (tsrc & QPMS_COORDS_BITRANGE) {
case QPMS_COORDS_SPH: {
const sph_t *s = src;
const sph_t *s = (const sph_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_CART3: {
const cart3_t *s = src;
const cart3_t *s = (const cart3_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart2sph(s[i]);
return;
} break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2sph_equator(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart22sph(s[i]);
return;
} break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart12sph_zaxis(s[i]);
return;
@ -627,31 +628,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
cart3_t *d = (cart3_t *) dest;
switch (tsrc & QPMS_COORDS_BITRANGE) {
case QPMS_COORDS_SPH: {
const sph_t *s = src;
const sph_t *s = (const sph_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = sph2cart(s[i]);
return;
} break;
case QPMS_COORDS_CART3: {
const cart3_t *s = src;
const cart3_t *s = (const cart3_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2cart3_equator(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart22cart3xy(s[i]);
return;
} break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *)src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart12cart3z(s[i]);
return;
@ -669,13 +670,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed");
break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart2pol(s[i]);
return;
@ -696,13 +697,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed");
break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2cart(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
@ -727,7 +728,7 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 1D coordinates not allowed");
break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
@ -899,6 +900,6 @@ static inline cart2_t cart2_from_double_array(const double a[]) {
typedef double matrix3d[3][3];
typedef double matrix2d[2][2];
typedef complex double cmatrix3d[3][3];
typedef complex double cmatrix2d[2][2];
typedef _Complex double cmatrix3d[3][3];
typedef _Complex double cmatrix2d[2][2];
#endif //VECTORS_H

View File

@ -9,6 +9,7 @@
#include <string.h>
#include "qpms_error.h"
#include "normalisation.h"
#include <stdbool.h>
qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
@ -52,7 +53,7 @@ qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswf
return QPMS_SUCCESS;
}
bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
_Bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
const qpms_vswf_set_spec_t *b) {
if (a == b) return true;
if (a->norm != b->norm) return false;
@ -131,9 +132,9 @@ size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf
while(big_pairs[bi].ui < small_pairs[si].ui)
++bi;
if(big_pairs[bi].ui == small_pairs[si].ui)
r[small_pairs[si].i_orig] = big_pairs[si].i_orig;
r[small_pairs[si].i_orig] = big_pairs[bi].i_orig;
else
r[small_pairs[si].i_orig] = ~(size_t)0;
r[small_pairs[si].i_orig] = QPMS_VSWF_SET_REINDEX_UNMAPPED;
}
free(small_pairs);
@ -307,7 +308,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
qpms_y_t nelem = qpms_lMax2nelem(lMax);
csphvec_t *a;
QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t))
csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem;
csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a1 + 2 * nelem;
QPMS_ENSURE_SUCCESS(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm));
const csphvec_t *p1 = a1;
const csphvec_t *p2 = a2;
@ -478,7 +479,7 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex
}
qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, const void *args, bool add) {
const cart3_t evalpoint, const void *args, _Bool add) {
QPMS_UNTESTED;
const qpms_incfield_planewave_params_t *p = args;

View File

@ -7,6 +7,10 @@
*/
#ifndef QPMS_VSWF_H
#define QPMS_VSWF_H
#ifdef __cplusplus
extern "C" {
#endif
#include <unistd.h> // ssize_t
#include "qpms_types.h"
#include <gsl/gsl_sf_legendre.h>
@ -16,7 +20,7 @@
/// Calculates the (regular VSWF) expansion coefficients of an external incident field.
typedef qpms_errno_t (*qpms_incfield_t)(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
complex double *target,
_Complex double *target,
const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
const void *args, ///< Pointer to additional function-specific arguments.
@ -61,12 +65,14 @@ static inline ssize_t qpms_vswf_set_spec_find_uvswfi(const qpms_vswf_set_spec_t
* It's not lossless if the two bspecs contain different combinations of waves.
*
* Preferably, big->ilist contains everything small->ilist does.
* If small->ilist[i] is not found in big->ilist, r[i] will be set to ~(size_t)0.
* If small->ilist[i] is not found in big->ilist, r[i] will be set to
* \ref QPMS_VSWF_SET_REINDEX_UNMAPPED = ~(size_t)0.
*
* Discard with free() after use.
*/
size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf_set_spec_t *big);
static const size_t QPMS_VSWF_SET_REINDEX_UNMAPPED = ~(size_t)0;
/// Evaluates a set of VSWF basis functions at a given point.
/** The list of basis wave indices is specified in \a setspec;
@ -82,7 +88,7 @@ qpms_errno_t qpms_uvswf_fill(
/** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist.
*/
csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
const complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
const _Complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
csph_t kr, ///< Evaluation point.
qpms_bessel_t btyp);
@ -118,7 +124,7 @@ typedef struct qpms_incfield_planewave_params_t {
*/
qpms_errno_t qpms_incfield_planewave(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
complex double *target,
_Complex double *target,
const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)).
@ -227,23 +233,26 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir, ccart3_t amplitude,
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
qpms_l_t lMax, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
qpms_l_t lMax, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf(sph_t where,
complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf_csph(csph_t where,
complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI
qpms_bessel_t btyp, qpms_normalisation_t norm);
void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);
#ifdef __cplusplus
}
#endif
#endif // QPMS_VSWF_H

View File

@ -52,7 +52,7 @@ complex double qpms_wignerD_elem(const qpms_quat_t R,
}
complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q,
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) {
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, _Bool pseudo) {
#ifndef NDEBUG
qpms_irot3_checkdet(q);
#endif

View File

@ -111,8 +111,8 @@ cywaves = Extension('qpms.cywaves',
libraries = common_libs,
include_dirs = [numpy_includes],
)
cytranslations = Extension('qpms.cytranslations',
sources = ['qpms/cytranslations.pyx',
cytranslations = Extension('qpms.cytranslations_legacy',
sources = ['qpms/cytranslations_legacy.pyx',
'qpms/translations_python.c',
],
extra_compile_args=['-std=c99',
@ -150,7 +150,7 @@ qpms_c = Extension('qpms.qpms_c',
sources = ['qpms/qpms_c.pyx',],
libraries = common_libs,
include_dirs=['amos', 'qpms', numpy_includes],
extra_link_args=['-fopenmp'],
extra_link_args=['-fopenmp'], # needed?
extra_compile_args=['-fopenmp'],
)

View File

@ -0,0 +1,7 @@
find_package(Catch2 REQUIRED)
add_executable(catchtest all_includes.C t_vswf.C catch_aux.C)
# TODO ensure that the tests are linked against the "working tree" version of qpms,
# not the installed one
target_link_libraries(catchtest PRIVATE Catch2::Catch2WithMain qpms lapacke)

View File

@ -0,0 +1,30 @@
// Just to test that all headers are C++ compatible
#include <catch2/catch.hpp>
#include <vswf.h>
#include <beyn.h>
#include <gaunt.h>
#include <groups.h>
#include <kahansum.h>
#include <materials.h>
#include <optim.h>
#include <oshacks.h>
#include <parsing.h>
#include <qpms_error.h>
#include <qpms_specfunc.h>
#include <qpms_types.h>
#include <scatsystem.h>
#include <symmetries.h>
#include <tiny_inlines.h>
#include <tmatrices.h>
#include <tolerances.h>
#include <translations.h>
#include <vectors.h>
#include <pointgroups.h>
#include <lattices_types.h>
#include <indexing.h>
#include <ewald.h>
// C++ type conversion issues:
#include <quaternions.h>
//#include <lattices.h>
//#include <normalisation.h>
//#include <qpmsblas.h>

24
tests/catch/catch_aux.C Normal file
View File

@ -0,0 +1,24 @@
#include "catch_aux.h"
#include "complex.h"
namespace qpmstest{
/// Creates a new vector<double> from complex number array
std::vector<double> pointer2dvec(const _Complex double *arr, size_t siz){
std::vector<double> vec(2*siz);
for(size_t i = 0; i < siz; ++i) {
vec[2*i] = creal(arr[i]);
vec[2*i+1] = cimag(arr[i]);
}
return vec;
}
std::vector<double> pointer2dvec(const double *arr, size_t siz){
std::vector<double> vec(siz);
for(size_t i = 0; i < siz; ++i)
vec[i] = arr[i];
return vec;
}
std::vector<double> pointer2dvec(const csphvec_t *arr, size_t siz) {
return pointer2dvec(&arr->rc, 3*siz);
}
}

44
tests/catch/t_vswf.C Normal file
View File

@ -0,0 +1,44 @@
#include <catch2/catch.hpp>
#include <qpms_error.h>
#include <vswf.h>
#include "catch_aux.h"
using namespace qpmstest;
constexpr qpms_l_t lMax = 6;
static std::vector<double> vswfvals_dvector(double r, double theta, double phi,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
{
std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2))
* 3 /* 3D vector components */ * 2 /*for real and imag*/);
qpms_vswf_fill(
nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
lMax, sph_t{r, theta, phi},
btyp, norm);
return vals;
}
static std::vector<double> vswfvals_dvector_alternative(double r, double theta, double phi,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
{
std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2))
* 3 /* 3D vector components */ * 2 /*for real and imag*/);
qpms_vswf_fill_alternative(
nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
lMax, sph_t{r, theta, phi},
btyp, norm);
return vals;
}
TEST_CASE("VSWF alternative implementation") {
using Catch::Matchers::Approx;
auto norms = GENERATE(QPMS_NORMALISATION_DEFAULT, QPMS_NORMALISATION_NORM_POWER,
QPMS_NORMALISATION_NORM_SPHARM, QPMS_NORMALISATION_CONVENTION_SCUFF);
auto thetas = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
auto phis = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
auto rs_positive = GENERATE(0.0001, 0.001, 0.01, 0.1, 1., 10., 100.);
CHECK_THAT(
vswfvals_dvector(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms),
Approx(vswfvals_dvector_alternative(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms))
) ;
}

65
tests/test_ewald_consistency.py Executable file
View File

@ -0,0 +1,65 @@
#!/usr/bin/env python3
'''
Test 2D Ewald summation consistency for basis VSWFs, using varying values of eta.
'''
import sys
# Test parameters, rather arbitrary
npoints = 10
etafactors = [0.9, 1.1]
zpointmax = 0. if len(sys.argv) == 1 else float(sys.argv[1])
lattice_basis=[(580e-9,0),(0,580e-9)]
positions = [(0,0)] # particle positions
wavevector = [0.,0.,0.]
bgparam = 1.52
omega_eh = 2.1
lMax = 2
import numpy as np
from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, EwaldPart, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar, c
from qpms.symmetries import point_group_info
eh = eV/hbar
np.random.seed(666)
omega = omega_eh * eh
bspec = BaseSpec(lMax=lMax)
medium = EpsMuGenerator(bgparam)
emg_dummy = EpsMuGenerator(10)
tmg_dummy = TMatrixGenerator.sphere(medium, emg_dummy, 50e-9)
particles = [Particle(pos, tmg_dummy, bspec) for pos in positions]
ss, ssw = ScatteringSystem.create(particles, medium, omega, latticebasis=lattice_basis)
wavevector = np.array(wavevector)
sswk = ssw._sswk(wavevector)
eta_orig = sswk.eta
print("default eta: %g" % (eta_orig,))
print("lattice basis:\n", ss.lattice_basis)
points = (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[0] + (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[1]
points += (np.random.rand(npoints) - .5)[:,None] * np.array([0.,0.,zpointmax])
fields_reference = sswk.scattered_field_basis(points)
fields = np.empty((len(etafactors),) + fields_reference.shape, dtype=complex)
fails = 0
for n, etafactor in enumerate(etafactors):
sswk.eta = eta_orig * etafactor
fields[n] = sswk.scattered_field_basis(points)
print(fields[n])
print(fields[n]-fields_reference)
if np.sum(1-np.isclose(fields[n], fields_reference)) > 0:
fails += 1
print(fails)
sys.exit(int(fails))