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

Merge branch 'catch_test' into ewald_dbg_rebased2

This commit is contained in:
Marek Nečada 2022-04-24 22:26:44 +03:00
parent 736adb6974 16a0a7dc0a
commit cae211e426
44 changed files with 614 additions and 354 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
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)

35
qpms/beyn.h
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

16
qpms/cytmatrices.pyx
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)

1
qpms/ewald.c
View File

@ -7,6 +7,7 @@
#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>

68
qpms/ewald.h
View File

@ -26,14 +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"
#include "lattices_types.h"
#include <complex.h>
/// Use this handler to ignore underflows of incomplete gamma.
@ -51,7 +55,7 @@ typedef struct qpms_ewald3_constants_t {
/// 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:
*
@ -65,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
@ -78,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:
*
@ -92,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
@ -100,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] =
*
@ -108,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.
@ -131,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);
@ -146,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);
@ -172,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.
/**
@ -184,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.
/**
@ -201,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.
@ -218,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.
@ -237,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.
@ -255,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,
@ -285,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,
@ -302,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,
@ -317,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,
@ -336,4 +340,8 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
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

1
qpms/ewaldsf.c
View File

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

7
qpms/gaunt.h
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

8
qpms/groups.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

12
qpms/indexing.h
View File

@ -26,6 +26,9 @@
*/
#ifndef QPMS_INDEXING_H
#define QPMS_INDEXING_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <math.h>
@ -106,7 +109,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
@ -119,7 +122,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;
@ -132,7 +135,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;
@ -147,4 +150,7 @@ static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_INDEXING_H

16
qpms/kahansum.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

25
qpms/lattices.h
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

1
qpms/materials.c
View File

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

42
qpms/materials.h
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

41
qpms/normalisation.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

6
qpms/parsing.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

1
qpms/pointgroups.c
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;

10
qpms/pointgroups.h
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

1
qpms/qpms_cdefs.pxd
View File

@ -493,6 +493,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

6
qpms/qpms_error.h
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

17
qpms/qpms_specfunc.h
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

25
qpms/qpms_types.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)
@ -191,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).
@ -212,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.
@ -261,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.
@ -346,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;
@ -422,8 +425,8 @@ 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
@ -436,4 +439,8 @@ typedef enum {
} qpms_ewald_part;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#ifdef __cplusplus
}
#endif
#endif // QPMS_TYPES

7
qpms/qpmsblas.h
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

18
qpms/quaternions.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

1
qpms/scatsystem.c
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"

167
qpms/scatsystem.h
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,
@ -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.
);
@ -734,7 +738,7 @@ ccart3_t qpms_scatsys_scattered_E(
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.
);
@ -753,7 +757,7 @@ 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.
_Complex double wavenumber, ///< Wavenumber of the background medium.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -773,7 +777,7 @@ qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
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
_Complex double wavenumber, ///< Wavenumber of the background medium
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -832,8 +836,8 @@ qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
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.
);
@ -849,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.
);
@ -868,14 +872,14 @@ 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$.
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
);
@ -903,7 +907,7 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
/// 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)
@ -913,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

2
qpms/scatsystem_dbg.h
View File

@ -1,7 +1,7 @@
#include "scatsystem.h"
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size
_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.

34
qpms/symmetries.h
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

9
qpms/tiny_inlines.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:

6
qpms/tmatrices.c
View File

@ -988,6 +988,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:

125
qpms/tmatrices.h
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

6
qpms/tolerances.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

5
qpms/translations.c
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.

70
qpms/translations.h
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

1
qpms/translations_dbg.c
View File

@ -5,6 +5,7 @@
#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,

8
qpms/translations_dbg.h
View File

@ -8,16 +8,16 @@
#include "translations.h"
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
complex double *dest, const csph_t kdlj, const qpms_bessel_t J);
_Complex double *dest, const csph_t kdlj, const qpms_bessel_t J);
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);
_Complex double *dest, const _Complex double k, const cart3_t r, 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,
_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,

4
qpms/translations_inlines.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

1
qpms/trivialgroup.c
View File

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

53
qpms/vectors.h
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

7
qpms/vswf.c
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;
@ -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;

21
qpms/vswf.h
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.
@ -82,7 +86,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 +122,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 +231,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
#ifdef __cplusplus
}
#endif
#endif // QPMS_VSWF_H

2
qpms/wigner.c
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

7
tests/catch/CMakeLists.txt Normal file
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)

30
tests/catch/all_includes.C Normal file
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>

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))
) ;
}