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

Replace "complex" with "_Complex" in headers 

Needed for compatibility with C++ (in which "complex" is
a template).
This commit is contained in:
Marek Nečada 2022-03-30 09:06:47 +03:00
parent d640176d62
commit 5af1fd7734
26 changed files with 312 additions and 311 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
qpms/beyn.h
View File

@ -5,16 +5,15 @@
#define BEYN_H #define BEYN_H
#include <stddef.h> #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. /// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
/** Pure C array version */ /** 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$. /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
/** Pure C array version */ /** Pure C array version */
typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l, 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); const _Complex double *Vhat, _Complex double z, void *params);
/// Complex plane integration contour structure. /// Complex plane integration contour structure.
typedef struct beyn_contour_t { typedef struct beyn_contour_t {
@ -26,15 +25,15 @@ typedef struct beyn_contour_t {
* It does not have to be a centre in some strictly defined sense, * It does not have to be a centre in some strictly defined sense,
* but it should be "somewhere around" where the contour is. * 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). /// Function testing that a point \a z lies inside the contour (optional).
_Bool (*inside_test)(struct beyn_contour_t *, complex double z); _Bool (*inside_test)(struct beyn_contour_t *, _Complex double z);
complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. _Complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
} beyn_contour_t; } beyn_contour_t;
/// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes. /// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */ /** 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 { typedef enum {
@ -47,11 +46,11 @@ typedef enum {
/// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes. /// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */ /** Free using free(). */
beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints, 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_halfellipse_orientation or);
/// Similar to halfellipse but with rounded corners. /// 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) double rounding, ///< Must be in interval [0, 0.5)
size_t n, beyn_contour_halfellipse_orientation or); size_t n, beyn_contour_halfellipse_orientation or);
@ -60,10 +59,10 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, dou
typedef struct beyn_result_t { typedef struct beyn_result_t {
size_t neig; ///< Number of eigenvalues found. size_t neig; ///< Number of eigenvalues found.
size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec). size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
complex double *eigval; _Complex double *eigval;
complex double *eigval_err; _Complex double *eigval_err;
double *residuals; double *residuals;
complex double *eigvec; // Rows are the eigenvectors _Complex double *eigvec; // Rows are the eigenvectors
double *ranktest_SV; double *ranktest_SV;
} beyn_result_t; } beyn_result_t;

53
qpms/ewald.h
View File

@ -31,7 +31,6 @@
#include <gsl/gsl_sf_legendre.h> #include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_errno.h> #include <gsl/gsl_errno.h>
#include <math.h> // for inlined lilgamma #include <math.h> // for inlined lilgamma
#include <complex.h>
#include "qpms_types.h" #include "qpms_types.h"
#include "lattices.h" #include "lattices.h"
@ -51,7 +50,7 @@ typedef struct qpms_ewald3_constants_t {
/// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`. /// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`.
qpms_l_t *s1_jMaxes; qpms_l_t *s1_jMaxes;
/// The constant factors for the long range part of a 2D Ewald sum. /// 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) /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN: * for m + n EVEN:
* *
@ -65,7 +64,7 @@ typedef struct qpms_ewald3_constants_t {
* *
* s1_constfacs[y(m,n)][j] = 0 * 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) // 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. /// 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 /** If the summation points lie along a different direction, use the formula for
@ -78,7 +77,7 @@ typedef struct qpms_ewald3_constants_t {
// TODO indexing mechanisms // TODO indexing mechanisms
/// The constant factors for the long range part of a 2D Ewald sum. /// 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) /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN: * for m + n EVEN:
* *
@ -92,7 +91,7 @@ typedef struct qpms_ewald3_constants_t {
* *
* S1_constfacs[y(m,n)][j] = 0 * 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. /// 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 /** If the summation points lie along a different direction, use the formula for
* 2D sum with additional factor of * 2D sum with additional factor of
@ -100,7 +99,7 @@ typedef struct qpms_ewald3_constants_t {
*/ */
complex double **s1_constfacs_1Dz; _Complex double **s1_constfacs_1Dz;
/* These are the actual numbers now: /* These are the actual numbers now:
* s1_constfacs_1Dz[n][j] = * s1_constfacs_1Dz[n][j] =
* *
@ -108,7 +107,7 @@ typedef struct qpms_ewald3_constants_t {
* -------------------------- * --------------------------
* j! * 2**(2*j) * (n - 2*j)! * 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 double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
* what the multipliers from translations.c count with. * what the multipliers from translations.c count with.
@ -131,13 +130,13 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
/// Structure for holding complex-valued result of computation and an error estimate. /// Structure for holding complex-valued result of computation and an error estimate.
/** Similar to gsl_sf_result, but with complex val. */ /** Similar to gsl_sf_result, but with complex val. */
typedef struct qpms_csf_result { typedef struct qpms_csf_result {
complex double val; ///< Calculation result. _Complex double val; ///< Calculation result.
double err; ///< Error estimate. double err; ///< Error estimate.
} qpms_csf_result; } qpms_csf_result;
// [1, (A.9)] // [1, (A.9)]
static inline complex double lilgamma(double t) { static inline _Complex double lilgamma(double t) {
t = fabs(t); t = fabs(t);
if (t >= 1) if (t >= 1)
return sqrt(t*t - 1); return sqrt(t*t - 1);
@ -146,8 +145,8 @@ static inline complex double lilgamma(double t) {
} }
// [1, (A.8)], complex version of lilgamma() // [1, (A.8)], complex version of lilgamma()
static inline complex double clilgamma(complex double z) { static inline _Complex double clilgamma(_Complex double z) {
complex double a1 = z - 1, a2 = z + 1; _Complex double a1 = z - 1, a2 = z + 1;
// ensure -pi/2 < arg(z + 1) < 3*pi/2 // ensure -pi/2 < arg(z + 1) < 3*pi/2
if (creal(a2) < 0 && cimag(a2) <= 0) if (creal(a2) < 0 && cimag(a2) <= 0)
a2 = -csqrt(a2); a2 = -csqrt(a2);
@ -172,7 +171,7 @@ static inline complex double clilgamma(complex double z) {
* even if `z1 == z2`, because `-0 == 0` according to IEEE 754. * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`. * 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. /// Incomplete Gamma function as continued fractions.
/** /**
@ -184,7 +183,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. * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`. * 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. /// Incomplete gamma for complex second argument.
/** /**
@ -201,7 +200,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 * Another than principal branch can be selected using non-zero \a m
* argument. * argument.
*/ */
int complex_gamma_inc_e(double a, complex double x, int complex_gamma_inc_e(double a, _Complex double x,
/// Branch index. /// Branch index.
/** If zero, the principal value is calculated. /** If zero, the principal value is calculated.
* Other branches might be chosen using non-zero \a m. * Other branches might be chosen using non-zero \a m.
@ -218,7 +217,7 @@ int complex_gamma_inc_e(double a, complex double x,
/// Exponential integral for complex second argument. /// Exponential integral for complex second argument.
/** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */ /** 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. /// Hypergeometric 2F2, used to calculate some errors.
@ -237,15 +236,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
* unsuitable especially for big values of \a maxn. * 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, void ewald3_2_sigma_long_Delta(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, complex double z); int xbranch, _Complex double z);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /// 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. /** This function always uses Kambe's (corrected) recurrent formula.
* For production, use ewald3_2_sigma_long_Delta() instead. * 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, 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); int xbranch, _Complex double z, _Bool bigimz);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /// 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. /** This function always uses Taylor expansion in \a z.
@ -255,26 +254,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 * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
* of the error. * of the error.
*/ */
void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x, void ewald3_2_sigma_long_Delta_series(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, complex double z); int xbranch, _Complex double z);
// General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space // 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. /// 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). 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(). const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter. 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$. /// 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( 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`. 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(). const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter. double eta, ///< Ewald parameter.
complex double wavenumber, ///< Wavenumber of the background medium. _Complex double wavenumber, ///< Wavenumber of the background medium.
/// Lattice dimensionality. /// Lattice dimensionality.
/** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */ /** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
LatticeDimensionality latdim, LatticeDimensionality latdim,
@ -302,11 +301,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$. /// 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 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`. 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(). const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter. 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). double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality).
/// Lattice dimensionality. /// Lattice dimensionality.
LatticeDimensionality latdim, LatticeDimensionality latdim,

2
qpms/groups.h
View File

@ -35,7 +35,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 /** 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] * m[i * dim * dim + r * dim + c]
*/ */
complex double *m; _Complex double *m;
}; };
/// A point group with its irreducible representations and some metadata. /// A point group with its irreducible representations and some metadata.

10
qpms/kahansum.h
View File

@ -4,8 +4,6 @@
#ifndef KAHANSUM_H #ifndef KAHANSUM_H
#define KAHANSUM_H #define KAHANSUM_H
#include <complex.h>
static inline void kahaninit(double * const sum, double * const compensation) { static inline void kahaninit(double * const sum, double * const compensation) {
*sum = 0; *sum = 0;
*compensation = 0; *compensation = 0;
@ -19,14 +17,14 @@ 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; *sum = 0;
*compensation = 0; *compensation = 0;
} }
static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) { static inline void ckahanadd(_Complex double *sum, _Complex double *compensation, _Complex double input) {
complex double compensated_input = input - *compensation; _Complex double compensated_input = input - *compensation;
complex double nsum = *sum + compensated_input; _Complex double nsum = *sum + compensated_input;
*compensation = (nsum - *sum) - compensated_input; *compensation = (nsum - *sum) - compensated_input;
*sum = nsum; *sum = nsum;
} }

1
qpms/materials.c
View File

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

35
qpms/materials.h
View File

@ -4,6 +4,7 @@
#ifndef QPMS_MATERIALS_H #ifndef QPMS_MATERIALS_H
#define QPMS_MATERIALS_H #define QPMS_MATERIALS_H
#include "qpms_types.h" #include "qpms_types.h"
#include <complex.h>
#include <gsl/gsl_spline.h> #include <gsl/gsl_spline.h>
#ifndef SPEED_OF_LIGHT #ifndef SPEED_OF_LIGHT
@ -20,36 +21,36 @@ typedef struct qpms_epsmu_generator_t {
* qpms_permittivity_interpolator_epsmu_g(), * qpms_permittivity_interpolator_epsmu_g(),
* qpms_lorentzdrude_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; const void *params;
} qpms_epsmu_generator_t; } qpms_epsmu_generator_t;
/// Convenience function for generating material properties at given frequency. /// Convenience function for generating material properties at given frequency.
static inline qpms_epsmu_t qpms_epsmu_generator_eval( 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); return gen.function(omega, gen.params);
} }
/// Constant optical property "generator" for qpms_epsmu_generator_t. /// 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. const void *epsmu ///< Points to the qpms_epsmu_t to be returned.
); );
/// Gets refractive index of a material from its permeability and permittivity. /// Gets refractive index of a material from its permeability and permittivity.
/** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */ /** \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); return csqrt(em.eps * em.mu);
} }
/// Gets wave number \a k from angular frequency and material permeability and permittivity. /// 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] */ /** \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; return qpms_refindex(em)*omega/SPEED_OF_LIGHT;
} }
/// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity. /// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity.
/** \eta_r = \sqrt{\mu_r / \varepsilon_r} \f] */ /** \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); return csqrt(em.mu / em.eps);
} }
@ -67,7 +68,7 @@ typedef struct qpms_ldparams_triple_t {
* \f] * \f]
*/ */
typedef struct qpms_ldparams_t { typedef struct qpms_ldparams_t {
complex double eps_inf; ///< Permittivity at infinity. _Complex double eps_inf; ///< Permittivity at infinity.
double omega_p; ///< Plasma frequency. double omega_p; ///< Plasma frequency.
size_t n; ///< Number of "oscillators". size_t n; ///< Number of "oscillators".
qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters. qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters.
@ -86,14 +87,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. extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten.
/// Lorentz-Drude permittivity. /// 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. /// 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. /// Lorentz-Drude optical properties, with relative permeability set always to one, compatible with qpms_epsmu_generator_t.
qpms_epsmu_t qpms_lorentzdrude_epsmu_g( 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 *. const void *ldparams ///< Lorentz-Drude parameters, in reality const qpms_ldparams_t *.
); );
@ -125,14 +126,14 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
); );
/// Evaluates interpolated material permittivity at a given angular frequency. /// 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); const qpms_permittivity_interpolator_t *interp, double omega_SI);
/// Evaluates interpolated material permittivity at a given angular frequency, qpms_epsmu_generator_t compatible version. /// 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. /** Permeability is always set to one. Imaginary part of omega is discarded.
*/ */
qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g( 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 const void * interpolator ///< Interpolator of type qpms_permittivity_interpolator_t
); );
@ -148,11 +149,11 @@ double qpms_permittivity_interpolator_omega_max(
void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp); void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
/// Relative permittivity from the Drude model. /// Relative permittivity from the Drude model.
static inline complex double qpms_drude_epsilon( static inline _Complex double qpms_drude_epsilon(
complex double eps_inf, ///< Relative permittivity "at infinity". _Complex double eps_inf, ///< Relative permittivity "at infinity".
complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material. _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 gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated. _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)); return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
} }

33
qpms/normalisation.h
View File

@ -9,7 +9,6 @@
#include "qpms_types.h" #include "qpms_types.h"
#include "qpms_error.h" #include "qpms_error.h"
#include <math.h> #include <math.h>
#include <complex.h>
#include "indexing.h" #include "indexing.h"
#include "optim.h" #include "optim.h"
@ -36,8 +35,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 * 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.) * 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) { 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); _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_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -51,8 +50,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 * Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call. * 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) { 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); _Complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
} }
@ -62,8 +61,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 * 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.) * 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) { 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); _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_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -77,14 +76,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 * Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call. * 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) { 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); _Complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; 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. /// 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) return qpms_normalisation_factor_N_noCS(norm, l, m)
/ qpms_normalisation_factor_M_noCS(norm, l, m); / qpms_normalisation_factor_M_noCS(norm, l, m);
} }
@ -95,8 +94,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 * 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.) * 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) { 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); _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_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -109,13 +108,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 * Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call. * 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) { 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); _Complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; 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. /// 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_vswf_type_t t; qpms_m_t m; qpms_l_t l;
qpms_uvswfi2tmn(ui, &t, &m, &l); qpms_uvswfi2tmn(ui, &t, &m, &l);
switch(t) { switch(t) {
@ -156,7 +155,7 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t
* 0 \quad \mbox{if } m>0. \\ * 0 \quad \mbox{if } m>0. \\
* \f] * \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) switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
case 0: case 0:
@ -194,7 +193,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; if(m==0) return 0;
switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {

1
qpms/pointgroups.c
View File

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

4
qpms/pointgroups.h
View File

@ -18,9 +18,9 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
case QPMS_PGS_DN: case QPMS_PGS_DN:
case QPMS_PGS_DND: case QPMS_PGS_DND:
case QPMS_PGS_DNH: case QPMS_PGS_DNH:
return true; return 1;
default: default:
return false; return 0;
} }
} }

10
qpms/qpms_specfunc.h
View File

@ -11,25 +11,25 @@
******************************************************************************/ ******************************************************************************/
// TODO unify types // 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 { typedef struct {
qpms_l_t lMax; qpms_l_t lMax;
double *akn; // coefficients as in DLMF 10.49.1 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_t *qpms_sbessel_calculator_init(void); qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(void);
void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c); 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, 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, 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);
/****************************************************************************** /******************************************************************************

17
qpms/qpms_types.h
View File

@ -3,10 +3,9 @@
*/ */
#ifndef QPMS_TYPES_H #ifndef QPMS_TYPES_H
#define QPMS_TYPES_H #define QPMS_TYPES_H
#include <complex.h>
#include <stdbool.h>
#include <stddef.h> #include <stddef.h>
#include <stdint.h> #include <stdint.h>
#include <stdbool.h>
#ifndef M_PI_2 #ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L) #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L)
@ -191,7 +190,7 @@ typedef struct cart3_t {
/// 3D complex (actually 6D) coordinates. See also vectors.h. /// 3D complex (actually 6D) coordinates. See also vectors.h.
typedef struct ccart3_t { typedef struct ccart3_t {
complex double x, y, z; _Complex double x, y, z;
} ccart3_t; } ccart3_t;
/// 3D complex vector pair (represents the E, H fields). /// 3D complex vector pair (represents the E, H fields).
@ -212,13 +211,13 @@ typedef struct sph_t {
/// Spherical coordinates with complex radial component. See also vectors.h. /// Spherical coordinates with complex radial component. See also vectors.h.
typedef struct csph_t { // Do I really need this??? typedef struct csph_t { // Do I really need this???
complex double r; _Complex double r;
double theta, phi; double theta, phi;
} csph_t; } csph_t;
/// 3D complex vector components in local spherical basis. See also vectors.h. /// 3D complex vector components in local spherical basis. See also vectors.h.
typedef struct csphvec_t { typedef struct csphvec_t {
complex double rc, thetac, phic; _Complex double rc, thetac, phic;
} csphvec_t; } csphvec_t;
/// 2D polar coordinates. See also vectors.h. /// 2D polar coordinates. See also vectors.h.
@ -261,7 +260,7 @@ typedef enum {
* See quaternions.h for "methods". * See quaternions.h for "methods".
*/ */
typedef struct qpms_quat_t { typedef struct qpms_quat_t {
complex double a, b; _Complex double a, b;
} qpms_quat_t; } qpms_quat_t;
/// Quaternion type as four doubles. /// Quaternion type as four doubles.
@ -346,7 +345,7 @@ typedef struct qpms_tmatrix_t {
* there is also one with order \a −m. * there is also one with order \a −m.
*/ */
const qpms_vswf_set_spec_t *spec; 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() bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
} qpms_tmatrix_t; } qpms_tmatrix_t;
@ -422,8 +421,8 @@ typedef struct qpms_abstract_tmatrix_t {
/// A type holding electric permittivity and magnetic permeability of a material. /// A type holding electric permittivity and magnetic permeability of a material.
typedef struct qpms_epsmu_t { typedef struct qpms_epsmu_t {
complex double eps; ///< Relative permittivity. _Complex double eps; ///< Relative permittivity.
complex double mu; ///< Relative permeability. _Complex double mu; ///< Relative permeability.
} qpms_epsmu_t; } qpms_epsmu_t;
struct qpms_tolerance_spec_t; // See tolerances.h struct qpms_tolerance_spec_t; // See tolerances.h

12
qpms/quaternions.h
View File

@ -93,7 +93,7 @@ static inline double qpms_quat_norm(const qpms_quat_t q) {
} }
/// Test approximate equality of quaternions. /// 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; return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
} }
@ -119,7 +119,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$. /// 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( return qpms_quat_norm(qpms_quat_sub(
qpms_quat_standardise(p, atol), qpms_quat_standardise(p, atol),
qpms_quat_standardise(q, atol))) <= atol; qpms_quat_standardise(q, atol))) <= atol;
@ -202,17 +202,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 * The D matrix are calculated using formulae (3), (4), (6), (7) from
* http://moble.github.io/spherical_functions/WignerDMatrices.html * 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); qpms_m_t mp, qpms_m_t m);
/// A VSWF representation element of the O(3) group. /// A VSWF representation element of the O(3) group.
/** /**
* TODO more doc. * 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. 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, 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 +251,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
} }
/// Test approximate equality of irot3. /// 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; return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
} }

160
qpms/scatsystem.h
View File

@ -107,7 +107,7 @@ typedef struct qpms_ss_orbit_type_t {
* *
* TODO doc. * TODO doc.
*/ */
complex double *irbases; _Complex double *irbases;
/// TODO doc. /// TODO doc.
size_t instance_count; size_t instance_count;
/// Cumulative sum of the preceding ot->siza * ot->instance_count; /// Cumulative sum of the preceding ot->siza * ot->instance_count;
@ -243,9 +243,9 @@ typedef struct qpms_scatsys_at_omega_t {
* i.e in the order corresponding to \a ss->tm. * i.e in the order corresponding to \a ss->tm.
*/ */
qpms_tmatrix_t **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 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; } qpms_scatsys_at_omega_t;
@ -291,7 +291,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, 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. /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
void qpms_scatsys_free(qpms_scatsys_t *s); void qpms_scatsys_free(qpms_scatsys_t *s);
@ -303,50 +303,50 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
/// Evaluates scattering system T-matrices at a given frequency. /// Evaluates scattering system T-matrices at a given frequency.
/** Free the result using qpms_scatsys_at_omega_free() when done. */ /** 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, 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. /// 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. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
* *
* TODO doc about shape etc. * 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); const qpms_scatsys_t *ss, qpms_iri_t iri);
/// Projects a "big" matrix onto an irrep (slow reference implementation). /// Projects a "big" matrix onto an irrep (slow reference implementation).
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed, _Complex double *qpms_scatsys_irrep_pack_matrix_stupid(_Complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss, const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri); qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis (slow reference implementation). /// Transforms a big "packed" matrix into the full basis (slow reference implementation).
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full, _Complex double *qpms_scatsys_irrep_unpack_matrix_stupid(_Complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss, const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add); qpms_iri_t iri, bool add);
/// Projects a "big" matrix onto an irrep. /// Projects a "big" matrix onto an irrep.
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, _Complex double *qpms_scatsys_irrep_pack_matrix(_Complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss, const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri); qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis. /// Transforms a big "packed" matrix into the full basis.
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full, _Complex double *qpms_scatsys_irrep_unpack_matrix(_Complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss, const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add); qpms_iri_t iri, bool add);
/// Projects a "big" vector onto an irrep. /// Projects a "big" vector onto an irrep.
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed, _Complex double *qpms_scatsys_irrep_pack_vector(_Complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss, const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri); qpms_iri_t iri);
/// Transforms a big "packed" vector into the full basis. /// Transforms a big "packed" vector into the full basis.
/** TODO doc */ /** TODO doc */
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, _Complex double *qpms_scatsys_irrep_unpack_vector(_Complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss, const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add); qpms_iri_t iri, bool add);
/// Global translation matrix. /// Global translation matrix.
@ -354,31 +354,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). * The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves).
* This may change in the future. * 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. /// 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, 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. /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
/** /**
* \returns \a target on success, NULL on error. * \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. /// 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 const qpms_scatsys_at_omega_t *ssw
); );
/// As qpms_scatsys_build_translation_full() but with choice of Bessel function type. /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
/** Might be useful for evaluation of cross sections and testing. /** 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. /// 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, 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 qpms_bessel_t J
); );
@ -387,12 +387,12 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
* The diagonal (particle self-) blocks are currently filled with zeros. * The diagonal (particle self-) blocks are currently filled with zeros.
* This may change in the future. * 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. /// 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, const qpms_scatsys_t *ss,
qpms_iri_t iri, 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 qpms_bessel_t J
); );
@ -400,9 +400,9 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
/** /**
* \returns \a target on success, NULL on error. * \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. /// 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, const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
); );
@ -410,9 +410,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
/** /**
* \returns \a target on success, NULL on error. * \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. /// 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, const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
); );
@ -420,9 +420,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
/** /**
* \returns \a target on success, NULL on error. * \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. /// 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, const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
); );
@ -435,7 +435,7 @@ typedef struct qpms_ss_LU {
bool full; ///< true if full matrix; false if irrep-packed. bool full; ///< true if full matrix; false if irrep-packed.
qpms_iri_t iri; ///< Irrep index if `full == false`. qpms_iri_t iri; ///< Irrep index if `full == false`.
/// LU decomposition array. /// LU decomposition array.
complex double *a; _Complex double *a;
/// Pivot index array, size at least max(1,min(m, n)). /// Pivot index array, size at least max(1,min(m, n)).
int *ipiv; int *ipiv;
} qpms_ss_LU; } qpms_ss_LU;
@ -443,14 +443,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. /// 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( 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw const qpms_scatsys_at_omega_t *ssw
); );
/// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri qpms_iri_t iri
@ -458,7 +458,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. /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( 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). 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 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. const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
@ -466,16 +466,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. /// 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( 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri 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$. /// 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 *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. _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`). 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. qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data.
); );
@ -495,33 +495,33 @@ typedef struct qpms_scatsys_at_omega_k_t {
/** /**
* \returns \a target on success, NULL on error. * \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. /// 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 const qpms_scatsys_at_omega_k_t *sswk
); );
/// Global translation matrix. /// 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. /// 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, 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. 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. double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
); );
/// Global translation matrix. /// 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. /// 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 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. /// 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( 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). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_k_t *sswk const qpms_scatsys_at_omega_k_t *sswk
); );
@ -539,7 +539,7 @@ struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
/// Wavevector in cartesian coordinates (must lie in the lattice plane). /// Wavevector in cartesian coordinates (must lie in the lattice plane).
const double k[3], 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_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. 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) size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -561,12 +561,12 @@ struct qpms_finite_group_t;
/// Constructs a "full matrix action" of a point group element for an orbit type. /// Constructs a "full matrix action" of a point group element for an orbit type.
/** TODO detailed doc */ /** TODO detailed doc */
complex double *qpms_orbit_action_matrix( _Complex double *qpms_orbit_action_matrix(
/// Target array. If NULL, a new one is allocated. /// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2 /** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec). * (it makes sense to assume all the T-matrices share their spec).
*/ */
complex double *target, _Complex double *target,
/// The orbit (type). /// The orbit (type).
const qpms_ss_orbit_type_t *orbit, const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has /// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -579,12 +579,12 @@ complex double *qpms_orbit_action_matrix(
/// Constructs a dense matrix representation of a irrep projector for an orbit type. /// Constructs a dense matrix representation of a irrep projector for an orbit type.
/** TODO detailed doc */ /** 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. /// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2 /** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec). * (it makes sense to assume all the T-matrices share their spec).
*/ */
complex double *target, _Complex double *target,
/// The orbit (type). /// The orbit (type).
const qpms_ss_orbit_type_t *orbit, const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has /// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -596,14 +596,14 @@ complex double *qpms_orbit_irrep_projector_matrix(
const qpms_iri_t iri); const qpms_iri_t iri);
/// TODO DOC!!!!! /// TODO DOC!!!!!
complex double *qpms_orbit_irrep_basis( _Complex double *qpms_orbit_irrep_basis(
/// Here theh size of theh basis shall be saved, /// Here theh size of theh basis shall be saved,
size_t *basis_size, size_t *basis_size,
/// Target array. If NULL, a new one is allocated. /// Target array. If NULL, a new one is allocated.
/** The size of the array is basis_size * (orbit->size * bspec->n) /** The size of the array is basis_size * (orbit->size * bspec->n)
* (it makes sense to assume all the T-matrices share their spec). * (it makes sense to assume all the T-matrices share their spec).
*/ */
complex double *target, _Complex double *target,
/// The orbit (type). /// The orbit (type).
const qpms_ss_orbit_type_t *orbit, const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has /// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -618,10 +618,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. /// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc! /** TODO detailed doc!
* \returns target_full if target_full was not NULL, otherwise the newly allocated array. */ * \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. /// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */ /** The length of the array is ss->fecv_size. */
complex double *target_full, _Complex double *target_full,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_incfield_t field_at_point, qpms_incfield_t field_at_point,
const void *args, ///< Pointer passed as the last argument to (*field_at_point)() const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
@ -629,9 +629,9 @@ complex double *qpms_scatsys_incident_field_vector_full(
); );
/// Applies T-matrices onto an incident field vector in the full basis. /// Applies T-matrices onto an incident field vector in the full basis.
complex double *qpms_scatsysw_apply_Tmatrices_full( _Complex double *qpms_scatsysw_apply_Tmatrices_full(
complex double *target_full, /// Target vector array. If NULL, a new one is allocated. _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 _Complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
const qpms_scatsys_at_omega_t *ssw const qpms_scatsys_at_omega_t *ssw
); );
@ -650,7 +650,7 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
/// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system. /// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
qpms_iri_t iri, 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_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. 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) size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -673,7 +673,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
double eta, ///< Ewald sum parameter double eta, ///< Ewald sum parameter
const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system. 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_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. 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) size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -687,10 +687,10 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
#if 0 #if 0
/// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points. /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc! */ /** 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. /// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */ /** The length of the array is ss->fecv_size. */
complex double *target_full, _Complex double *target_full,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym. const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym.
qpms_incfield_t field_at_point, qpms_incfield_t field_at_point,
@ -715,8 +715,8 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
ccart3_t qpms_scatsys_scattered_E( ccart3_t qpms_scatsys_scattered_E(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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.
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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -734,7 +734,7 @@ ccart3_t qpms_scatsys_scattered_E(
ccart3_t qpms_scatsysw_scattered_E( ccart3_t qpms_scatsysw_scattered_E(
const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -753,7 +753,7 @@ qpms_errno_t qpms_scatsys_scattered_field_basis(
ccart3_t *target, ///< Target array of length \a ss->fecv_size ccart3_t *target, ///< Target array of length \a ss->fecv_size
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -773,7 +773,7 @@ qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_ss_pi_t pi, ///< Particle index qpms_ss_pi_t pi, ///< Particle index
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -832,8 +832,8 @@ qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
ccart3_t qpms_scatsys_scattered_E__alt( ccart3_t qpms_scatsys_scattered_E__alt(
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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.
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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -849,7 +849,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(
ccart3_t qpms_scatsysw_scattered_E__alt( ccart3_t qpms_scatsysw_scattered_E__alt(
const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). 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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
@ -868,14 +868,14 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
ccart3_t qpms_scatsyswk_scattered_E( ccart3_t qpms_scatsyswk_scattered_E(
const qpms_scatsys_at_omega_k_t *sswk, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported). 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. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
ccart3_t qpms_scatsyswk_scattered_E_e( ccart3_t qpms_scatsyswk_scattered_E_e(
const qpms_scatsys_at_omega_k_t *sswk, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported). 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. cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
qpms_ewald_part parts qpms_ewald_part parts
); );
@ -903,7 +903,7 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
/// Adjusted Ewadl parameter to avoid high-frequency breakdown. /// Adjusted Ewadl parameter to avoid high-frequency breakdown.
// TODO DOC // 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 #if 0
/** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
@ -913,9 +913,9 @@ 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, ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, ///< Irreducible representation 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. cart3_t where, ///< Evaluation point.
complex double k ///< Wave number. _Complex double k ///< Wave number.
); );
#endif #endif

2
qpms/scatsystem_dbg.h
View File

@ -1,7 +1,7 @@
#include "scatsystem.h" #include "scatsystem.h"
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e( 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, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted). 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. cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.

26
qpms/symmetries.h
View File

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

9
qpms/tiny_inlines.h
View File

@ -4,6 +4,7 @@
#ifndef TINY_INLINES_H #ifndef TINY_INLINES_H
#define TINY_INLINES_H #define TINY_INLINES_H
#include <stdlib.h> #include <stdlib.h>
#include <complex.h>
static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; } 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 #if 0
#ifdef __GSL_SF_LEGENDRE_H__ #ifdef __GSL_SF_LEGENDRE_H__
static inline complex double 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) { 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; 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 #endif
// this has shitty precision: // 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; x = ((x % 4) + 4) % 4;
switch(x) { switch(x) {
case 0: case 0:

108
qpms/tmatrices.h
View File

@ -14,7 +14,7 @@ struct qpms_finite_group_t;
typedef struct qpms_finite_group_t 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. /// 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); return t->m + (t->spec->n * rowno);
} }
@ -42,7 +42,7 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t);
* This function actually checks for identical vswf specs. * This function actually checks for identical vswf specs.
* TODO define constants with "default" atol, rtol for this function. * 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); const double rtol, const double atol);
/// Creates a T-matrix from another matrix and a symmetry operation. /// Creates a T-matrix from another matrix and a symmetry operation.
@ -53,7 +53,7 @@ bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
*/ */
qpms_tmatrix_t *qpms_tmatrix_apply_symop( qpms_tmatrix_t *qpms_tmatrix_apply_symop(
const qpms_tmatrix_t *T, ///< the original T-matrix 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. /// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
@ -64,7 +64,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
*/ */
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix 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. /// Symmetrizes a T-matrix with an involution symmetry operation.
@ -75,7 +75,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
*/ */
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution( qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
const qpms_tmatrix_t *T, ///< the original T-matrix 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. /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
@ -105,7 +105,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
*/ */
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace( qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix 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. /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
@ -149,7 +149,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
double **freqs_su, ///< Frequencies in SCUFF units (optional). double **freqs_su, ///< Frequencies in SCUFF units (optional).
/// The resulting T-matrices (optional). /// The resulting T-matrices (optional).
qpms_tmatrix_t **tmatrices_array, 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. /// Tells whether qpms_load_scuff_tmatrix should crash if fopen() fails.
@ -161,7 +161,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
* This is desirable e.g. when used in Python (so that proper exception can * This is desirable e.g. when used in Python (so that proper exception can
* be thrown). * 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. /// Loads a scuff-tmatrix generated file.
/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a
@ -189,7 +189,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* is accessed via * is accessed via
* (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci]. * (*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. /// In-place application of point group elements on raw T-matrix data.
@ -200,7 +200,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* TODO more doc. * TODO more doc.
*/ */
qpms_errno_t qpms_symmetrise_tmdata_irot3arr( 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, const qpms_vswf_set_spec_t *bspec,
size_t n_symops, size_t n_symops,
const qpms_irot3_t *symops const qpms_irot3_t *symops
@ -213,7 +213,7 @@ qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
* TODO more doc. * TODO more doc.
*/ */
qpms_errno_t qpms_symmetrise_tmdata_finite_group( 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_vswf_set_spec_t *bspec,
const qpms_finite_group_t *pointgroup const qpms_finite_group_t *pointgroup
); );
@ -242,9 +242,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$. /// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$.
complex double *qpms_apply_tmatrix( _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. _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 _Complex double *a, ///< Incident field coefficient array of size T->spec->n.
const qpms_tmatrix_t *T ///< T-matrix \a T to apply. const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
); );
@ -257,7 +257,7 @@ complex double *qpms_apply_tmatrix(
*/ */
typedef struct qpms_tmatrix_generator_t { typedef struct qpms_tmatrix_generator_t {
qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill. 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 ///< Implementation dependent parameters.
); );
const void *params; ///< Parameter pointer passed to the function. const void *params; ///< Parameter pointer passed to the function.
@ -267,7 +267,7 @@ typedef struct qpms_tmatrix_generator_t {
qpms_tmatrix_t *qpms_tmatrix_init_from_generator( qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
const qpms_vswf_set_spec_t *bspec, const qpms_vswf_set_spec_t *bspec,
qpms_tmatrix_generator_t gen, qpms_tmatrix_generator_t gen,
complex double omega); _Complex double omega);
/// Implementation of qpms_matrix_generator_t that just copies a constant matrix. /// Implementation of qpms_matrix_generator_t that just copies a constant matrix.
@ -275,7 +275,7 @@ qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
* the same base spec. * the same base spec.
*/ */
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
complex double omega, _Complex double omega,
/// Source T-matrix, real type is (const qpms_tmatrix_t*). /// Source T-matrix, real type is (const qpms_tmatrix_t*).
const void *tmatrix_orig const void *tmatrix_orig
); );
@ -295,7 +295,7 @@ extern const gsl_interp_type * gsl_interp_steffen;
// struct gsl_interp_accel; // use if lookup proves to be too slow // struct gsl_interp_accel; // use if lookup proves to be too slow
typedef struct qpms_tmatrix_interpolator_t { typedef struct qpms_tmatrix_interpolator_t {
const qpms_vswf_set_spec_t *bspec; 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_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_spline **splines_imag; ///< There will be a spline object for each nonzero element
// gsl_interp_accel **accel_real; // gsl_interp_accel **accel_real;
@ -317,7 +317,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. 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 double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
const gsl_interp_type *iptype 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 +325,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! /** 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. 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 *. const void *interpolator ///< Parameter of type qpms_tmatrix_interpolator_t *.
); );
@ -342,14 +342,14 @@ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matri
* *
* TODO better doc. * TODO better doc.
*/ */
complex double *qpms_mie_coefficients_reflection( _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 *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. const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
double a, ///< Radius of the sphere. double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material 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 k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material. _Complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e, ///< Relative permeability of the surrounding medium. _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_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. qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS.
); );
@ -357,10 +357,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. /// 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. 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. double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material 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 k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material. _Complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium. _Complex double mu_e ///< Relative permeability of the surrounding medium.
); );
/// Parameter structure for qpms_tmatrix_generator_sphere(). /// Parameter structure for qpms_tmatrix_generator_sphere().
@ -372,7 +372,7 @@ typedef struct qpms_tmatrix_generator_sphere_param_t {
/// T-matrix generator for spherical particles using Lorentz-Mie solution. /// T-matrix generator for spherical particles using Lorentz-Mie solution.
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, 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. const void *params ///< Of type qpms_tmatrix_generator_sphere_param_t.
); );
@ -380,10 +380,10 @@ qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
static inline qpms_tmatrix_t *qpms_tmatrix_spherical( static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
const qpms_vswf_set_spec_t *bspec, const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere. double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material 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 k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material. _Complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium. _Complex double mu_e ///< Relative permeability of the surrounding medium.
) { ) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e); qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
@ -395,8 +395,8 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
double a, ///< Radius of the sphere. double a, ///< Radius of the sphere.
double omega, ///< Angular frequency. double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere. _Complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium. _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. /// Convenience function to calculate T-matrix of a non-magnetic spherical particle using the permittivity values.
@ -404,8 +404,8 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
const qpms_vswf_set_spec_t *bspec, const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere. double a, ///< Radius of the sphere.
double omega, ///< Angular frequency. double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere. _Complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium. _Complex double epsilon_bg ///< Relative permittivity of the background medium.
) { ) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg); qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
@ -460,7 +460,7 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta,
*/ */
qpms_errno_t qpms_tmatrix_axialsym_fill( qpms_errno_t qpms_tmatrix_axialsym_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. 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 outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material. qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation. qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -474,7 +474,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill(
/// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry. /// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry.
static inline qpms_tmatrix_t *qpms_tmatrix_axialsym( static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
const qpms_vswf_set_spec_t *bspec, 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 outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material. qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation. qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -501,22 +501,22 @@ typedef struct qpms_tmatrix_generator_axialsym_param_t {
/// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t /// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t
qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to fill. 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. 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). /// 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, qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(_Complex double *target,
complex double omega, _Complex double omega,
const qpms_tmatrix_generator_axialsym_param_t *param, const qpms_tmatrix_generator_axialsym_param_t *param,
qpms_normalisation_t norm, qpms_normalisation_t norm,
qpms_bessel_t J qpms_bessel_t J
); );
/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging). /// 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, qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(_Complex double *target,
complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside, _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_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$. 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 +537,7 @@ typedef struct qpms_tmatrix_function_t {
/// Convenience function to create a new T-matrix from qpms_tmatrix_function_t. /// Convenience function to create a new T-matrix from qpms_tmatrix_function_t.
// FIXME the name is not very intuitive. // 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); return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
} }
@ -558,9 +558,9 @@ typedef enum {
struct qpms_tmatrix_operation_lrmatrix { struct qpms_tmatrix_operation_lrmatrix {
/// Raw matrix data of \a M in row-major order. /// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */ /** The matrix must be taylored for the given bspec! */
complex double *m; _Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). 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; _Bool owns_m; ///< Whether \a m is owned by this;
}; };
struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations. struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
@ -591,9 +591,9 @@ struct qpms_tmatrix_operation_compose_chain {
struct qpms_tmatrix_operation_scmulz { struct qpms_tmatrix_operation_scmulz {
/// Raw matrix data of \a M in row-major order. /// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */ /** The matrix must be taylored for the given bspec! */
complex double *m; _Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). 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. _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. /// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t.
@ -602,7 +602,7 @@ struct qpms_tmatrix_operation_scmulz {
struct qpms_tmatrix_operation_irot3arr { struct qpms_tmatrix_operation_irot3arr {
size_t n; ///< Number of rotoreflections; size_t n; ///< Number of rotoreflections;
qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. 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. /// A generic T-matrix transformation operator.

1
qpms/translations.c
View File

@ -1,5 +1,6 @@
#include <math.h> #include <math.h>
#include "qpms_types.h" #include "qpms_types.h"
#include <complex.h>
#include "qpms_specfunc.h" #include "qpms_specfunc.h"
#include "gaunt.h" #include "gaunt.h"
#include "translations.h" #include "translations.h"

63
qpms/translations.h
View File

@ -28,7 +28,6 @@
#define QPMS_TRANSLATIONS_H #define QPMS_TRANSLATIONS_H
#include "vectors.h" #include "vectors.h"
#include "qpms_types.h" #include "qpms_types.h"
#include <complex.h>
#include <stdbool.h> #include <stdbool.h>
#include <stddef.h> #include <stddef.h>
@ -37,10 +36,10 @@
#endif #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); 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); bool r_ge_d, qpms_bessel_t J);
/// Structure holding the constant factors in normalisation operators. /// Structure holding the constant factors in normalisation operators.
@ -65,8 +64,8 @@ typedef struct qpms_trans_calculator {
qpms_normalisation_t normalisation; qpms_normalisation_t normalisation;
qpms_l_t lMax; qpms_l_t lMax;
qpms_y_t nelem; qpms_y_t nelem;
complex double **A_multipliers; _Complex double **A_multipliers;
complex double **B_multipliers; _Complex double **B_multipliers;
#if 0 #if 0
// Normalised values of the Legendre functions and derivatives // Normalised values of the Legendre functions and derivatives
// for θ == π/2, i.e. for the 2D case. // for θ == π/2, i.e. for the 2D case.
@ -90,44 +89,44 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation
/// Destructor for qpms_trans_calculator_t. /// Destructor for qpms_trans_calculator_t.
void qpms_trans_calculator_free(qpms_trans_calculator *); 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, 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); 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, 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); bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, 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, 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); bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, 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, size_t deststride, size_t srcstride,
csph_t kdlj, bool r_ge_d, qpms_bessel_t J); csph_t kdlj, bool r_ge_d, qpms_bessel_t J);
// TODO update the types later // TODO update the types later
complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, _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, int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J); 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, _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, int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J); 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, int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest, _Complex double *Adest, _Complex double *Bdest,
int m, int n, int mu, int nu, complex double kdlj_r, int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J); 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, 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, 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); int r_ge_d, int J);
// Convenience functions using VSWF base specs // Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c, 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. /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
@ -138,12 +137,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
/// and with automatic \a r_ge_d = `false`. /// and with automatic \a r_ge_d = `false`.
qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
const qpms_trans_calculator *c, const qpms_trans_calculator *c,
complex double *target, _Complex double *target,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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 qpms_bessel_t J
/// Workspace has to be at least 2 * c->neleme**2 long /// Workspace has to be at least 2 * c->neleme**2 long
); );
@ -156,10 +155,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr, _Complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr, _Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride, 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, cart2_t b1, cart2_t b2,
const cart2_t beta, const cart2_t beta,
const cart3_t particle_shift, const cart3_t particle_shift,
@ -167,10 +166,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, int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr, _Complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr, _Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride, 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, cart2_t b1, cart2_t b2,
const cart2_t beta, const cart2_t beta,
const cart3_t particle_shift, const cart3_t particle_shift,
@ -180,12 +179,12 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
// Convenience functions using VSWF base specs // Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c, 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. /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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, cart2_t b1, cart2_t b2,
const cart2_t beta, const cart2_t beta,
const cart3_t particle_shift, const cart3_t particle_shift,
@ -193,12 +192,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, 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. /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride, const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride, 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, cart2_t b1, cart2_t b2,
const cart2_t beta, const cart2_t beta,
const cart3_t particle_shift, const cart3_t particle_shift,
@ -212,8 +211,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; // e31z means that the particles are positioned along the z-axis;
// their positions and K-values are then denoted by a single z-coordinate // 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, int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr, _Complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr, _Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride, const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const double k, const double eta, const double k,
const double unitcell_area, // just lattice period const double unitcell_area, // just lattice period

8
qpms/translations_dbg.h
View File

@ -8,16 +8,16 @@
#include "translations.h" #include "translations.h"
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c, 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, 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 #ifdef LATTICESUMS32
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c, int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
complex double * const sigmas_total, double * serr_total, _Complex double * const sigmas_total, double * serr_total,
const double eta, const complex double k, const double eta, const _Complex double k,
const cart2_t b1, const cart2_t b2, const cart2_t b1, const cart2_t b2,
const cart2_t beta, const cart2_t beta,
const cart3_t particle_shift, 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. /// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix.
// TODO DOC // TODO DOC
static inline void qpms_trans_array_from_AB( static inline void qpms_trans_array_from_AB(
complex double *t, _Complex double *t,
const qpms_vswf_set_spec_t *const t_destspec, const qpms_vswf_set_spec_t *const t_destspec,
const size_t t_deststride, const size_t t_deststride,
const qpms_vswf_set_spec_t *const t_srcspec, const qpms_vswf_set_spec_t *const t_srcspec,
const size_t t_srcstride, 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. /// A and B matrices' lMax.
/** This also determines their size and stride: they are assumed to /** This also determines their size and stride: they are assumed to
* be square matrices of size `nelem * nelem` where * be square matrices of size `nelem * nelem` where

1
qpms/trivialgroup.c
View File

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

17
qpms/vectors.h
View File

@ -4,6 +4,7 @@
#ifndef VECTORS_H #ifndef VECTORS_H
#define VECTORS_H #define VECTORS_H
#include <math.h> #include <math.h>
#include <complex.h>
#ifndef M_PI_2 #ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487) #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
#endif #endif
@ -198,12 +199,12 @@ static inline double cart3_dist(const cart3_t a, const cart3_t b) {
return cart3norm(cart3_substract(a,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; return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
} }
/// Complex 3D vector scaling. /// 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}; ccart3_t res = {c * v.x, c * v.y, c * v.z};
return res; 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. /// 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; 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. /// 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}; csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
return res; return res;
} }
/// Complex 3D vector (geographic coordinates) "dot product" without conjugation. /// 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 //N.B. no complex conjugation done here
return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic; 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. /// "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}; csph_t res = {c * s.r, s.theta, s.phi};
return res; return res;
} }
@ -899,6 +900,6 @@ static inline cart2_t cart2_from_double_array(const double a[]) {
typedef double matrix3d[3][3]; typedef double matrix3d[3][3];
typedef double matrix2d[2][2]; typedef double matrix2d[2][2];
typedef complex double cmatrix3d[3][3]; typedef _Complex double cmatrix3d[3][3];
typedef complex double cmatrix2d[2][2]; typedef _Complex double cmatrix2d[2][2];
#endif //VECTORS_H #endif //VECTORS_H

5
qpms/vswf.c
View File

@ -9,6 +9,7 @@
#include <string.h> #include <string.h>
#include "qpms_error.h" #include "qpms_error.h"
#include "normalisation.h" #include "normalisation.h"
#include <stdbool.h>
qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() { 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; 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) { const qpms_vswf_set_spec_t *b) {
if (a == b) return true; if (a == b) return true;
if (a->norm != b->norm) return false; if (a->norm != b->norm) return false;
@ -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, 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; QPMS_UNTESTED;
const qpms_incfield_planewave_params_t *p = args; const qpms_incfield_planewave_params_t *p = args;

14
qpms/vswf.h
View File

@ -16,7 +16,7 @@
/// Calculates the (regular VSWF) expansion coefficients of an external incident field. /// Calculates the (regular VSWF) expansion coefficients of an external incident field.
typedef qpms_errno_t (*qpms_incfield_t)( typedef qpms_errno_t (*qpms_incfield_t)(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n. /// 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 qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made. const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
const void *args, ///< Pointer to additional function-specific arguments. const void *args, ///< Pointer to additional function-specific arguments.
@ -82,7 +82,7 @@ qpms_errno_t qpms_uvswf_fill(
/** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist. /** 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, 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. csph_t kr, ///< Evaluation point.
qpms_bessel_t btyp); qpms_bessel_t btyp);
@ -118,7 +118,7 @@ typedef struct qpms_incfield_planewave_params_t {
*/ */
qpms_errno_t qpms_incfield_planewave( qpms_errno_t qpms_incfield_planewave(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n. /// 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 qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made. 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 *)). const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)).
@ -227,19 +227,19 @@ 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_l_t lMax, sph_t dir, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir, ccart3_t amplitude, 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_l_t lMax, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude, 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); qpms_l_t lMax, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf(sph_t where, 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); qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf_csph(csph_t where, 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_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_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,

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, 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 #ifndef NDEBUG
qpms_irot3_checkdet(q); qpms_irot3_checkdet(q);
#endif #endif