Mention submodule in README
[qpms.git] / qpms / scatsystem.h
blob7325a6e5e5039a85220a5ed64373812ff301883e
1 /*! \file scatsystem.h
2 * \brief Modern interface for finite lattice calculations, including symmetries.
4 * N.B. Only "reasonably normalised" waves are supported now in most of the
5 * functions defined here, i.e. those that can be rotated by the usual
6 * Wigner matrices, i.e. the "power" or "spharm" -normalised ones.
7 *
8 * TODO FIXME check whether Condon-Shortley phase can have some nasty influence
9 * here; I fear that yes.
11 #ifndef QPMS_SCATSYSTEM_H
12 #define QPMS_SCATSYSTEM_H
13 #include "qpms_types.h"
14 #include "vswf.h"
15 #include "tmatrices.h"
16 #include <stdbool.h>
18 /// Overrides the number of threads spawned by the paralellized functions.
19 /** TODO MORE DOC which are those? */
20 void qpms_scatsystem_set_nthreads(long n);
22 /// A particle, defined by its T-matrix and position.
23 /** This is rather only an auxillary intermediate structure to ultimately
24 * build an qpms_scatsys_t instance */
25 typedef struct qpms_particle_t {
26 // Does it make sense to ever use other than cartesian coords for this?
27 cart3_t pos; ///< Particle position in cartesian coordinates.
28 const qpms_tmatrix_function_t *tmg; ///< T-matrix function; not owned by qpms_particle_t.
29 qpms_tmatrix_operation_t op; ///< T-matrix transformation operation w.r.t. \a tmg.
30 } qpms_particle_t;
32 struct qpms_finite_group_t;
33 typedef struct qpms_finite_group_t qpms_finite_group_t;
35 /// A particle, defined by its T-matrix INDEX and position, to be used in qpms_scatsys_t.
36 typedef struct qpms_particle_tid_t {
37 // Does it make sense to ever use other than cartesian coords for this?
38 cart3_t pos; ///< Particle position in cartesian coordinates.
39 qpms_ss_tmi_t tmatrix_id; ///< T-matrix index
40 } qpms_particle_tid_t;
43 typedef qpms_gmi_t qpms_ss_orbit_pi_t; ///< Auxilliary type used in qpms_ss_orbit_type_t for labeling particles inside orbits.
44 typedef qpms_ss_tmi_t qpms_ss_oti_t; ///< Auxilliary type used for labeling orbit types.
46 /// Structure describing a particle's "orbit type" under symmetry group actions in a system.
47 /**
48 * Each particle has its orbit with respect to a symmetry group of the system in which the particle lies,
49 * i.e. a set of particles onto which the group operations map the original particle.
51 * (TODO DOC improve the following explanation:)
52 * Typically, there will be only a small number of different (T-matrix, particle
53 * <a href="https://en.wikipedia.org/wiki/Group_action#Fixed_points_and_stabilizer_subgroups">stabiliser</a>)
54 * pairs in the system. We can group the particles accordingly, into the same "orbit types"
55 * that will allow to do certain operations only once for each "orbit type", saving memory and (perhaps) time.
57 * Each particle will then have assigned:
58 * 1. an orbit type,
59 * 2. an ID inside that orbit.
62 * TODO DOC how the process of assigning the particle IDs actually work, orbit type (non-)uniqueness.
65 * Memory is managed by qpms_scatspec_t; qpms_ss_orbit_type_t does not own anything.
68 typedef struct qpms_ss_orbit_type_t {
69 qpms_ss_orbit_pi_t size; ///< Size of the orbit (a divisor of the group order), i.e. number of particles on the orbit.
70 size_t bspecn; ///< Individual particle's coefficient vector length. The same as ss->tm[this.tmatrices[0]]->spec->n.
71 /// Action of the group elements onto the elements in this orbit.
72 /** Its size is sym->order * this.size
73 * and its values lie between 0 and \a this.size − 1.
75 * Action of the group element g onto the pi-th particle
76 * is given by action[g + pi*sym->order].
79 qpms_ss_orbit_pi_t *action;
80 /// T-matrix IDs of the particles on this orbit (type).
81 /**
82 * We save all the tmi's of the particles on the orbit here to make the number of array lookups
83 * and pointer dereferences constant.
85 * The size of this array is \a size.
87 qpms_ss_tmi_t *tmatrices;
88 /// Sizes of the per-orbit irrep bases.
89 /**
90 * The order of the irreps corresponds to the order in \a ss->sym->irreps.
91 * The size of this array is (obviously) \a ss->sym->nirreps.
93 * TODO different type?
94 * TODO doc.
96 size_t *irbase_sizes;
97 //The following are pretty redundant, TODO reduce them at some point.
98 /// Cumulative sums of irbase_sizes.
99 size_t *irbase_cumsizes;
100 /// TODO doc.
101 size_t *irbase_offsets;
102 /// Per-orbit irreducible representation orthonormal bases.
103 /** This also defines the unitary operator that transforms the orbital excitation coefficients
104 * in the symmetry-adapted basis.
106 * The size is (\a this->size * \a this->tmatrices[0].spec->n)**2.
108 * TODO doc.
110 complex double *irbases;
111 /// TODO doc.
112 size_t instance_count;
113 /// Cumulative sum of the preceding ot->siza * ot->instance_count;
114 qpms_ss_pi_t p_offset;
115 } qpms_ss_orbit_type_t;
118 typedef ptrdiff_t qpms_ss_osn_t; ///< "serial number" of av orbit in a given type.
120 /// Auxillary type used in qpms_scatsys_t that identifies the particle's orbit and its id inside that orbit.
121 typedef struct qpms_ss_particle_orbitinfo {
122 qpms_ss_oti_t t; ///< Orbit type.
123 #define QPMS_SS_P_ORBITINFO_UNDEF (-1) ///< This labels that the particle has not yet been assigned to an orbit.
124 qpms_ss_osn_t osn; ///< "Serial number" of the orbit in the given type. TODO type and more doc.
125 qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type.
126 } qpms_ss_particle_orbitinfo_t;
128 /// Auxillary type used in qpms_scatsys_t: A recepy to create another T-matrices by symmetry operations.
129 typedef struct qpms_ss_derived_tmatrix_t {
130 qpms_ss_tmgi_t tmgi; ///< Index of the corresponding qpms_scatsys_t::tm element.
131 struct qpms_tmatrix_operation_t op; ///< Operation to derive this particular T-matrix.
132 } qpms_ss_derived_tmatrix_t;
134 typedef struct qpms_scatsys_periodic_info_t {
135 /// (Direct) lattice basis of the system (only \a lattice_dimension elements are used)
136 /** This is mandatory for \a lattice_dimension != 0 */
137 cart3_t lattice_basis[3];
139 /// Reciprocal lattice basis.
140 /**(TODO specify whether it includes 2π or not) */
141 cart3_t reciprocal_basis[3];
143 /// Unitcell volume (irrelevant for non-periodic systems).
144 /** The dimensionality of the volume corresponds to lattice_dimension, so
145 * for lattice_dimension == 1, this will actually be lenght and for
146 * lattice_dimension == 2, a 2D area.
148 double unitcell_volume;
150 /// Default Ewald parameter \f$ \eta \f$.
151 /** Normally, this just gets copied into qpms_scatsys_at_omega_t,
152 * which is then used in the Ewald sums.
153 * However, for higher frequencies it must be adjusted to avoid
154 * numerical instability.
156 double eta;
157 } qpms_scatsys_periodic_info_t;
160 struct qpms_trans_calculator;
161 struct qpms_epsmu_generator_t;
163 /// Common "class" for system of scatterers, both periodic and non-periodic.
165 * Infinite periodic structures (those with \a lattice_dimension > 0)
166 * have the \a per filled.
167 * These are ignored for finite systems (lattice_dimension == 0).
169 typedef struct qpms_scatsys_t {
170 /// Number of dimensions in which the system is periodic from the range 0–3.
171 int lattice_dimension;
172 struct qpms_epsmu_generator_t medium; ///< Optical properties of the background medium.
174 /// (Template) T-matrix functions in the system.
175 /** The qpms_abstract_tmatrix_t objects (onto which this array member point)
176 * are NOT owned by this and must be kept alive for the whole lifetime
177 * of all qpms_scatsys_t objects that are built upon them.
179 qpms_tmatrix_function_t *tmg;
180 qpms_ss_tmgi_t tmg_count; ///< Number of all different original T-matrix generators in the system.
182 /// All the different T-matrix functions in the system, including those derived from \a tmg elements by symmetries.
183 qpms_ss_derived_tmatrix_t *tm;
184 qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices in the system (length of tm[]).
185 qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[].
187 qpms_particle_tid_t *p; ///< Particles.
188 qpms_ss_pi_t p_count; ///< Size of particles array.
189 qpms_ss_pi_t p_capacity; ///< Capacity of p[].
191 //TODO the index types do not need to be so big.
192 const struct qpms_finite_group_t *sym; ///< Symmetry group of the array
193 qpms_ss_pi_t *p_sym_map; ///< Which particles map onto which by the symmetry ops.
194 ///< p_sym_map[idi + pi * sym->order] gives the index of pi-th particle under the idi'th sym op.
195 qpms_ss_tmi_t *tm_sym_map; ///< Which t-matrices map onto which by the symmetry ops. Lookup by tm_sum_map[idi + tmi * sym->order].
197 qpms_ss_oti_t orbit_type_count;
198 qpms_ss_orbit_type_t *orbit_types; ///< (Array length is \a orbit_type_count.)
200 qpms_ss_particle_orbitinfo_t *p_orbitinfo; ///< Orbit type identification of each particle. (Array length is \a p_count.)
202 size_t fecv_size; ///< Number of elements of a full excitation coefficient vector size.
203 size_t *saecv_sizes; ///< Number of elements of symmetry-adjusted coefficient vector sizes (order as in sym->irreps).
205 size_t *fecv_pstarts; ///< Indices of where pi'th particle's excitation coeffs start in a full excitation coefficient vector.
206 size_t *saecv_ot_offsets; ///< TODO DOC. In the packed vector, saecv_ot_offsets[iri * orbit_type_count + oti] indicates start of ot
207 /**< TODO maybe move it to qpms_ss_orbit_type_t, ffs. */
208 //size_t **saecv_pstarts; ///< NI. Indices of where pi'th particle's excitation coeff start in a symmetry-adjusted e.c.v.
209 ///**< First index is irrep index as in sym->irreps, second index is particle index. */
211 // TODO shifted origin of the symmetry group etc.
212 // TODO some indices for fast operations here.
213 // private
214 size_t max_bspecn; ///< Maximum tm[...]->spec->n. Mainly for workspace allocation.
216 /// Particles grouped by orbit (in the order corresponding to the packed memory layout).
217 qpms_ss_pi_t *p_by_orbit;
219 // We keep the p_orbitinfo arrays in this chunk in order to avoid memory fragmentation
220 char *otspace;
221 char *otspace_end;
222 double lenscale; // radius of the array, used as a relative tolerance measure
223 struct qpms_trans_calculator *c;
225 /// Periodic lattice metadata.
226 qpms_scatsys_periodic_info_t per;
227 } qpms_scatsys_t;
229 /// Retrieve the bspec of \a tmi'th element of \a ss->tm.
230 static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) {
231 return ss->tmg[ss->tm[tmi].tmgi].spec;
234 /// Retrieve the bspec of \a pi'th particle in \a ss->p.
235 static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) {
236 return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec;
239 typedef struct qpms_scatsys_at_omega_t {
240 const qpms_scatsys_t *ss; ///< Parent scattering system.
241 /// T-matrices from \a ss, evaluated at \a omega.
242 /** The T-matrices are in the same order as in \a ss,
243 * i.e in the order corresponding to \a ss->tm.
245 qpms_tmatrix_t **tm;
246 complex double omega; ///< Angular frequency
247 qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
248 complex double wavenumber; ///< Background medium wave number
249 } qpms_scatsys_at_omega_t;
252 /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed.
253 /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances,
254 * so keep them alive until scatsys is destroyed.
256 * The following fields must be filled in the "proto- scattering system" \a orig:
257 * * orig->lattice_dimension
258 * * orig->medium – The pointers are copied to the new qpms_scatsys_t instance;
259 * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting
260 * qpms_scatsys_t instances are properly destroyed.
261 * * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance;
262 * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting
263 * qpms_scatsys_t instances are properly destroyed. The pointers from orig->tmg, however, are copied.
264 * * orig->tmg_count
265 * * orig->tm – Must be filled, although the operations will typically be identities
266 * (QPMS_TMATRIX_OPERATION_NOOP). N.B. these NOOPs might be replaced with some symmetrisation operation
267 * in the resulting "full" qpms_scatsys_t instance.
268 * * orig->tm_count
269 * * orig->p
270 * * orig->p_count
272 * For periodic systems, the corresponding number of orig->per->lattice_basis[] elements
273 * must be filled as well.
275 * For periodic systems, only trivial group is currently supported. Non-trivial
276 * groups will cause undefined behaviour.
278 * The resulting qpms_scatsys_t is obtained by actually evaluating the T-matrices
279 * at the given frequency \a omega and where applicable, these are compared
280 * by their values with given tolerances. The T-matrix generators are expected
281 * to preserve the point group symmetries for all frequencies.
283 * This particular function will likely change in the future.
285 * \returns An instance \a sso of qpms_scatsys_omega_t. Note that \a sso->ss
286 * must be saved by the caller before destroying \a sso
287 * (with qpms_scatsys_at_omega_free(), and destroyed only afterwards with
288 * qpms_scatsys_free() when not needed anymore.
289 * \a sso->ss can be reused for different frequency by a
290 * qpms_scatsys_at_omega() call.
293 qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym,
294 complex double omega, const struct qpms_tolerance_spec_t *tol);
296 /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
297 void qpms_scatsys_free(qpms_scatsys_t *s);
299 /// Destroys a qpms_scatsys_at_omega_t.
300 /** Used on results of qpms_scatsys_apply_symmetry() and qpms_scatsys_at_omega(). */
301 void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
303 /// Evaluates scattering system T-matrices at a given frequency.
304 /** Free the result using qpms_scatsys_at_omega_free() when done. */
305 qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
306 complex double omega);
308 /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis.
309 /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
311 * TODO doc about shape etc.
313 complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
314 const qpms_scatsys_t *ss, qpms_iri_t iri);
316 /// Projects a "big" matrix onto an irrep (slow reference implementation).
317 /** TODO doc */
318 complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
319 const complex double *orig_full, const qpms_scatsys_t *ss,
320 qpms_iri_t iri);
322 /// Transforms a big "packed" matrix into the full basis (slow reference implementation).
323 /** TODO doc */
324 complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
325 const complex double *orig_packed, const qpms_scatsys_t *ss,
326 qpms_iri_t iri, bool add);
328 /// Projects a "big" matrix onto an irrep.
329 /** TODO doc */
330 complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
331 const complex double *orig_full, const qpms_scatsys_t *ss,
332 qpms_iri_t iri);
334 /// Transforms a big "packed" matrix into the full basis.
335 /** TODO doc */
336 complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
337 const complex double *orig_packed, const qpms_scatsys_t *ss,
338 qpms_iri_t iri, bool add);
340 /// Projects a "big" vector onto an irrep.
341 /** TODO doc */
342 complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
343 const complex double *orig_full, const qpms_scatsys_t *ss,
344 qpms_iri_t iri);
346 /// Transforms a big "packed" vector into the full basis.
347 /** TODO doc */
348 complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
349 const complex double *orig_packed, const qpms_scatsys_t *ss,
350 qpms_iri_t iri, bool add);
352 /// Global translation matrix.
354 * The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves).
355 * This may change in the future.
357 complex double *qpms_scatsys_build_translation_matrix_full(
358 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
359 complex double *target,
360 const qpms_scatsys_t *ss,
361 complex double k ///< Wave number to use in the translation matrix.
364 /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
366 * \returns \a target on success, NULL on error.
368 complex double *qpms_scatsysw_build_modeproblem_matrix_full(
369 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
370 complex double *target,
371 const qpms_scatsys_at_omega_t *ssw
374 /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
375 /** Might be useful for evaluation of cross sections and testing.
377 complex double *qpms_scatsys_build_translation_matrix_e_full(
378 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
379 complex double *target,
380 const qpms_scatsys_t *ss,
381 complex double k, ///< Wave number to use in the translation matrix.
382 qpms_bessel_t J
385 /// Global translation matrix with selectable Bessel function, projected on an irrep.
387 * The diagonal (particle self-) blocks are currently filled with zeros.
388 * This may change in the future.
390 complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
391 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
392 complex double *target,
393 const qpms_scatsys_t *ss,
394 qpms_iri_t iri,
395 complex double k, ///< Wave number to use in the translation matrix.
396 qpms_bessel_t J
399 /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form.
401 * \returns \a target on success, NULL on error.
403 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
404 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
405 complex double *target,
406 const qpms_scatsys_at_omega_t *ssw,
407 qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
409 /// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
411 * \returns \a target on success, NULL on error.
413 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
414 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
415 complex double *target,
416 const qpms_scatsys_at_omega_t *ssw,
417 qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
419 /// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
421 * \returns \a target on success, NULL on error.
423 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
424 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
425 complex double *target,
426 const qpms_scatsys_at_omega_t *ssw,
427 qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
430 struct qpms_scatsys_at_omega_k_t; // Defined below.
431 /// LU factorisation (LAPACKE_zgetrf) result holder.
432 typedef struct qpms_ss_LU {
433 const qpms_scatsys_at_omega_t *ssw;
434 const struct qpms_scatsys_at_omega_k_t *sswk; ///< Only for periodic systems, otherwise NULL.
435 bool full; ///< true if full matrix; false if irrep-packed.
436 qpms_iri_t iri; ///< Irrep index if `full == false`.
437 /// LU decomposition array.
438 complex double *a;
439 /// Pivot index array, size at least max(1,min(m, n)).
440 int *ipiv;
441 } qpms_ss_LU;
442 void qpms_ss_LU_free(qpms_ss_LU);
444 /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only.
445 qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
446 complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
447 int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
448 const qpms_scatsys_at_omega_t *ssw
451 /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
452 qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
453 complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
454 int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
455 const qpms_scatsys_at_omega_t *ssw,
456 qpms_iri_t iri
459 /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
460 qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
461 complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
462 int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
463 const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems.
464 const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
467 /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
468 qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
469 complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
470 int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
471 const qpms_scatsys_at_omega_t *ssw,
472 qpms_iri_t iri
475 /// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$.
476 complex double *qpms_scatsys_scatter_solve(
477 complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
478 const complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
479 qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data.
482 // ======================= Periodic system -only related stuff =============================
484 /// Scattering system at a given frequency and k-vector. Used only with periodic systems.
486 * N.B. use as a stack variable now, but this might become heap-allocated in the future (with own con- and destructor)
488 typedef struct qpms_scatsys_at_omega_k_t {
489 const qpms_scatsys_at_omega_t *ssw;
490 double k[3]; ///< The k-vector's cartesian coordinates.
491 double eta; ///< Ewald parameter η.
492 } qpms_scatsys_at_omega_k_t;
494 /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
496 * \returns \a target on success, NULL on error.
498 complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
499 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
500 complex double *target,
501 const qpms_scatsys_at_omega_k_t *sswk
504 /// Global translation matrix.
505 complex double *qpms_scatsys_periodic_build_translation_matrix_full(
506 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
507 complex double *target,
508 const qpms_scatsys_t *ss,
509 complex double wavenumber, ///< Wave number to use in the translation matrix.
510 const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates.
511 double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
514 /// Global translation matrix.
515 complex double *qpms_scatsyswk_build_translation_matrix_full(
516 /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
517 complex double *target,
518 const qpms_scatsys_at_omega_k_t *sswk
522 /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only.
523 qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
524 complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
525 int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
526 const qpms_scatsys_at_omega_k_t *sswk
529 /// Searches for periodic scattering system's eigenmodes using Beyn's algorithm.
531 * Currently, elliptical contour is used.
533 * TODO In the future, this will probably support irrep decomposition as well,
534 * but for the case of periodic / small systems,
535 * the bottleneck is the T-matrix and translation matrix evaluation
536 * rather than the linear algebra.
538 struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
539 const qpms_scatsys_t *ss,
540 /// Wavevector in cartesian coordinates (must lie in the lattice plane).
541 const double k[3],
542 complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
543 double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
544 double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
545 size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
546 double rank_tol, ///< (default: `1e-4`) TODO DOC.
547 size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
548 double res_tol ///< (default: `0.0`) TODO DOC.
551 // ======================= Periodic system -only related stuff end =========================
554 /// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file.
555 qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path);
557 /// NOT IMPLEMENTED Reads a qpms_scatsys_t structure from a file.
558 qpms_scatsys_t *qpms_scatsys_load(char *path);
560 struct qpms_finite_group_t;
562 /// Constructs a "full matrix action" of a point group element for an orbit type.
563 /** TODO detailed doc */
564 complex double *qpms_orbit_action_matrix(
565 /// Target array. If NULL, a new one is allocated.
566 /** The size of the array is (orbit->size * bspec->n)**2
567 * (it makes sense to assume all the T-matrices share their spec).
569 complex double *target,
570 /// The orbit (type).
571 const qpms_ss_orbit_type_t *orbit,
572 /// Base spec of the t-matrices (we don't know it from orbit, as it has
573 /// only T-matrix indices.
574 const qpms_vswf_set_spec_t *bspec,
575 /// The symmetry group used to generate the orbit (must have rep3d filled).
576 const struct qpms_finite_group_t *sym,
577 /// The index of the operation in sym to represent.
578 const qpms_gmi_t g);
580 /// Constructs a dense matrix representation of a irrep projector for an orbit type.
581 /** TODO detailed doc */
582 complex double *qpms_orbit_irrep_projector_matrix(
583 /// Target array. If NULL, a new one is allocated.
584 /** The size of the array is (orbit->size * bspec->n)**2
585 * (it makes sense to assume all the T-matrices share their spec).
587 complex double *target,
588 /// The orbit (type).
589 const qpms_ss_orbit_type_t *orbit,
590 /// Base spec of the t-matrices (we don't know it from orbit, as it has
591 /// only T-matrix indices.
592 const qpms_vswf_set_spec_t *bspec,
593 /// The symmetry group used to generate the orbit (must have rep3d filled).
594 const struct qpms_finite_group_t *sym,
595 /// The index of the irreducible representation of sym.
596 const qpms_iri_t iri);
598 /// TODO DOC!!!!!
599 complex double *qpms_orbit_irrep_basis(
600 /// Here theh size of theh basis shall be saved,
601 size_t *basis_size,
602 /// Target array. If NULL, a new one is allocated.
603 /** The size of the array is basis_size * (orbit->size * bspec->n)
604 * (it makes sense to assume all the T-matrices share their spec).
606 complex double *target,
607 /// The orbit (type).
608 const qpms_ss_orbit_type_t *orbit,
609 /// Base spec of the t-matrices (we don't know it from orbit, as it has
610 /// only T-matrix indices.
611 const qpms_vswf_set_spec_t *bspec,
612 /// The symmetry group used to generate the orbit (must have rep3d filled).
613 const struct qpms_finite_group_t *sym,
614 /// The index of the irreducible representation of sym.
615 const qpms_iri_t iri);
618 /// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points.
619 /** TODO detailed doc!
620 * \returns target_full if target_full was not NULL, otherwise the newly allocated array. */
621 complex double *qpms_scatsys_incident_field_vector_full(
622 /// Target array. If NULL, a new one is allocated.
623 /** The length of the array is ss->fecv_size. */
624 complex double *target_full,
625 const qpms_scatsys_t *ss,
626 qpms_incfield_t field_at_point,
627 const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
628 bool add ///< If true, add to target_full; rewrite target_full if false.
631 /// Applies T-matrices onto an incident field vector in the full basis.
632 complex double *qpms_scatsysw_apply_Tmatrices_full(
633 complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
634 const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
635 const qpms_scatsys_at_omega_t *ssw
638 struct beyn_result_t; // See beyn.h for full definition
640 /// Searches for finite scattering system's eigenmodes using Beyn's algorithm.
642 * Currently, elliptical contour is used.
644 * TODO In the future, this will probably support irrep decomposition as well,
645 * but it does not make much sense for periodic / small systems, as in their
646 * case the bottleneck is the T-matrix and translation matrix evaluation
647 * rather than the linear algebra.
649 struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
650 const qpms_scatsys_t *ss,
651 /// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
652 qpms_iri_t iri,
653 complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
654 double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
655 double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
656 size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
657 double rank_tol, ///< (default: `1e-4`) TODO DOC.
658 size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
659 double res_tol ///< (default: `0.0`) TODO DOC.
662 #if 0
663 /// Searches for scattering system's eigenmodes using Beyn's algorithm.
665 * Currently, elliptical contour is used.
667 * TODO In the future, this will probably support irrep decomposition as well,
668 * but it does not make much sense for periodic / small systems, as in their
669 * case the bottleneck is the T-matrix and translation matrix evaluation
670 * rather than the linear algebra.
672 struct beyn_result_t *qpms_scatsys_find_eigenmodes(
673 const qpms_scatsys_t *ss,
674 double eta, ///< Ewald sum parameter
675 const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system.
676 complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
677 double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
678 double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
679 size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
680 double rank_tol, ///< (default: `1e-4`) TODO DOC.
681 size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
682 double res_tol ///< (default: `0.0`) TODO DOC.
684 #endif
687 #if 0
688 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
689 /** TODO detailed doc! */
690 complex double *qpms_scatsys_incident_field_vector_irrep_packed(
691 /// Target array. If NULL, a new one is allocated.
692 /** The length of the array is ss->fecv_size. */
693 complex double *target_full,
694 const qpms_scatsys_t *ss,
695 const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym.
696 qpms_incfield_t field_at_point,
697 const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
698 bool add ///< If true, add to target_full; rewrite target_full if false.
700 #endif
702 /// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
704 * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any given wavenumber of the
705 * background medium and any given vector of the excitation coefficients \f$ \wckcout \f$.
707 * \return Complex electric field at the point defined by \a evalpoint.
709 * \see qpms_scatsysw_scattered_E()
711 * \see qpms_scatsyswk_scattered_E() for periodic systems.
715 ccart3_t qpms_scatsys_scattered_E(
716 const qpms_scatsys_t *ss,
717 qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
718 complex double wavenumber, ///< Wavenumber of the background medium.
719 const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
720 cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
723 /// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
725 * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any
726 * given vector of the excitation coefficients \f$ \wckcout \f$.
728 * \return Complex electric field at the point defined by \a evalpoint.
730 * \see qpms_scatsys_scattered_E()
732 * \see qpms_scatsyswk_scattered_E() for periodic systems.
734 ccart3_t qpms_scatsysw_scattered_E(
735 const qpms_scatsys_at_omega_t *ssw,
736 qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
737 const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
738 cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
741 /// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
743 * This is an alternative implementation of qpms_scatsys_scattered_E(), and should give the same results
744 * up to rounding errors.
746 * \return Complex electric field at the point defined by \a evalpoint.
748 * \see qpms_scatsys_scattered_E()
750 ccart3_t qpms_scatsys_scattered_E__alt(
751 const qpms_scatsys_t *ss,
752 qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
753 complex double wavenumber, ///< Wavenumber of the background medium.
754 const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
755 cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
758 /// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
760 * This is an alternative implementation of qpms_scatsys_scattered_E(), and should give the same results
761 * up to rounding errors.
763 * \return Complex electric field at the point defined by \a evalpoint.
765 * \see qpms_scatsysw_scattered_E()
767 ccart3_t qpms_scatsysw_scattered_E__alt(
768 const qpms_scatsys_at_omega_t *ssw,
769 qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
770 const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
771 cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
775 /// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients. (Periodic system.)
777 * This function evaluates the field \f$ \vect E (\vect r ) \f$, with any
778 * given vector of the excitation coefficients \f$ \wckcout \f$.
780 * \return Complex electric field at the point defined by \a evalpoint.
782 * \bug Currently implemented only for btyp == QPMS_HANKEL_PLUS.
784 * \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_E() for finite systems.
786 ccart3_t qpms_scatsyswk_scattered_E(
787 const qpms_scatsys_at_omega_k_t *sswk,
788 qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
789 const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
790 cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
793 /// Adjusted Ewadl parameter to avoid high-frequency breakdown.
794 // TODO DOC
795 double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
797 #if 0
798 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
799 * at a given point.
801 * \return Complex electric field at the point defined by \a where.
803 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
804 qpms_iri_t iri, ///< Irreducible representation
805 const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
806 cart3_t where, ///< Evaluation point.
807 complex double k ///< Wave number.
809 #endif
811 #endif //QPMS_SCATSYSTEM_H