diff --git a/CMakeLists.txt b/CMakeLists.txt index 7566bb8..6b6e18c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,6 +48,7 @@ add_subdirectory(faddeeva) add_subdirectory (qpms) +add_subdirectory (tests/catch EXCLUDE_FROM_ALL) enable_testing() add_subdirectory (tests EXCLUDE_FROM_ALL) diff --git a/qpms/beyn.h b/qpms/beyn.h index 13a43aa..9b05a3a 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -3,18 +3,20 @@ */ #ifndef BEYN_H #define BEYN_H +#ifdef __cplusplus +extern "C" { +#endif #include -#include /// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found. /** Pure C array version */ -typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params); +typedef int (*beyn_function_M_t)(_Complex double *target_M, size_t m, _Complex double z, void *params); /// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$. /** Pure C array version */ -typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l, - const complex double *Vhat, complex double z, void *params); +typedef int (*beyn_function_M_inv_Vhat_t)(_Complex double *target_M_inv_Vhat, size_t m, size_t l, + const _Complex double *Vhat, _Complex double z, void *params); /// Complex plane integration contour structure. typedef struct beyn_contour_t { @@ -26,15 +28,15 @@ typedef struct beyn_contour_t { * It does not have to be a centre in some strictly defined sense, * but it should be "somewhere around" where the contour is. */ - complex double centre; + _Complex double centre; /// Function testing that a point \a z lies inside the contour (optional). - _Bool (*inside_test)(struct beyn_contour_t *, complex double z); - complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. + _Bool (*inside_test)(struct beyn_contour_t *, _Complex double z); + _Complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. } beyn_contour_t; /// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes. /** Free using free(). */ -beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints); +beyn_contour_t *beyn_contour_ellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints); typedef enum { @@ -47,23 +49,23 @@ typedef enum { /// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes. /** Free using free(). */ -beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints, - beyn_contour_halfellipse_orientation or); +beyn_contour_t *beyn_contour_halfellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints, + beyn_contour_halfellipse_orientation o); /// Similar to halfellipse but with rounded corners. -beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im, +beyn_contour_t *beyn_contour_kidney(_Complex double centre, double halfax_re, double halfax_im, double rounding, ///< Must be in interval [0, 0.5) - size_t n, beyn_contour_halfellipse_orientation or); + size_t n, beyn_contour_halfellipse_orientation o); /// Beyn algorithm result structure (pure C array version). typedef struct beyn_result_t { size_t neig; ///< Number of eigenvalues found. size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec). - complex double *eigval; - complex double *eigval_err; + _Complex double *eigval; + _Complex double *eigval_err; double *residuals; - complex double *eigvec; // Rows are the eigenvectors + _Complex double *eigvec; // Rows are the eigenvectors double *ranktest_SV; } beyn_result_t; @@ -83,4 +85,7 @@ beyn_result_t *beyn_solve( double res_tol ///< (default: `0.0`) TODO DOC. ); +#ifdef __cplusplus +} +#endif #endif // BEYN_H diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index 91a22af..b447449 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -286,6 +286,10 @@ cdef class TMatrixGenerator: self.holder = what self.g.function = qpms_tmatrix_generator_interpolator self.g.params = (self.holder).rawpointer() + elif what == 0: + self.holder = 0 + self.g.function = qpms_tmatrix_generator_zero + self.g.params = 0 elif callable(what): warnings.warn("Custom python T-matrix generators are an experimental feature. Also expect it to be slow.") self.holder = what @@ -403,3 +407,15 @@ cdef class TMatrixGenerator: EpsMuGenerator(outside), EpsMuGenerator(inside), ArcFunction(__ArcCylinder(r, h)), *args, **kwargs)) + @staticmethod + def dummy(): + """Returns a dummy T-matrix generator (returns a zero T-matrix). + + Returns + ------- + tmgen_dummy : TMatrixGenerator + """ + return tmgen_dummy + +# pre-generate a dummy TMatrixGenerator (which returns zero T-matrix) +tmgen_dummy = TMatrixGenerator(0) diff --git a/qpms/ewald.c b/qpms/ewald.c index 8c727a1..cebf428 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -7,6 +7,7 @@ #include #include "tiny_inlines.h" #include "qpms_error.h" +#include "lattices.h" #include #include #include diff --git a/qpms/ewald.h b/qpms/ewald.h index 8881586..3d4ab7d 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -26,14 +26,18 @@ #ifndef EWALD_H #define EWALD_H +#ifdef __cplusplus +extern "C" { +#endif + #include #include #include #include #include // for inlined lilgamma -#include #include "qpms_types.h" -#include "lattices.h" +#include "lattices_types.h" +#include /// Use this handler to ignore underflows of incomplete gamma. @@ -51,7 +55,7 @@ typedef struct qpms_ewald3_constants_t { /// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`. qpms_l_t *s1_jMaxes; /// The constant factors for the long range part of a 2D Ewald sum. - complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] + _Complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version) * for m + n EVEN: * @@ -65,7 +69,7 @@ typedef struct qpms_ewald3_constants_t { * * s1_constfacs[y(m,n)][j] = 0 */ - complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors. + _Complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors. // similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0) /// The constant factors for the long range part of a 1D Ewald sum along the \a z axis. /** If the summation points lie along a different direction, use the formula for @@ -78,7 +82,7 @@ typedef struct qpms_ewald3_constants_t { // TODO indexing mechanisms /// The constant factors for the long range part of a 2D Ewald sum. - complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] + _Complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version) * for m + n EVEN: * @@ -92,7 +96,7 @@ typedef struct qpms_ewald3_constants_t { * * S1_constfacs[y(m,n)][j] = 0 */ - complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors. + _Complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors. /// The constant factors for the long range part of a 1D Ewald sum along the \a z axis. /** If the summation points lie along a different direction, use the formula for * 2D sum with additional factor of @@ -100,7 +104,7 @@ typedef struct qpms_ewald3_constants_t { */ - complex double **s1_constfacs_1Dz; + _Complex double **s1_constfacs_1Dz; /* These are the actual numbers now: * s1_constfacs_1Dz[n][j] = * @@ -108,7 +112,7 @@ typedef struct qpms_ewald3_constants_t { * -------------------------- * j! * 2**(2*j) * (n - 2*j)! */ - complex double *s1_constfacs_1Dz_base; ///= 1) return sqrt(t*t - 1); @@ -146,8 +150,8 @@ static inline complex double lilgamma(double t) { } // [1, (A.8)], complex version of lilgamma() -static inline complex double clilgamma(complex double z) { - complex double a1 = z - 1, a2 = z + 1; +static inline _Complex double clilgamma(_Complex double z) { + _Complex double a1 = z - 1, a2 = z + 1; // ensure -pi/2 < arg(z + 1) < 3*pi/2 if (creal(a2) < 0 && cimag(a2) <= 0) a2 = -csqrt(a2); @@ -172,7 +176,7 @@ static inline complex double clilgamma(complex double z) { * even if `z1 == z2`, because `-0 == 0` according to IEEE 754. * The side of the branch cut can be determined using `signbit(creal(z))`. */ -int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result); +int cx_gamma_inc_series_e(double a, _Complex double z, qpms_csf_result * result); /// Incomplete Gamma function as continued fractions. /** @@ -184,7 +188,7 @@ int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result); * even if `z1 == z2`, because `-0 == 0` according to IEEE 754. * The side of the branch cut can be determined using `signbit(creal(z))`. */ -int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result); +int cx_gamma_inc_CF_e(double a, _Complex double z, qpms_csf_result * result); /// Incomplete gamma for complex second argument. /** @@ -201,7 +205,7 @@ int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result); * Another than principal branch can be selected using non-zero \a m * argument. */ -int complex_gamma_inc_e(double a, complex double x, +int complex_gamma_inc_e(double a, _Complex double x, /// Branch index. /** If zero, the principal value is calculated. * Other branches might be chosen using non-zero \a m. @@ -218,7 +222,7 @@ int complex_gamma_inc_e(double a, complex double x, /// Exponential integral for complex second argument. /** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */ -int complex_expint_n_e(int n, complex double x, qpms_csf_result *result); +int complex_expint_n_e(int n, _Complex double x, qpms_csf_result *result); /// Hypergeometric 2F2, used to calculate some errors. @@ -237,15 +241,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result * unsuitable especially for big values of \a maxn. * */ -void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x, - int xbranch, complex double z); +void ewald3_2_sigma_long_Delta(_Complex double *target, double *target_err, int maxn, _Complex double x, + int xbranch, _Complex double z); /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /** This function always uses Kambe's (corrected) recurrent formula. * For production, use ewald3_2_sigma_long_Delta() instead. */ -void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x, - int xbranch, complex double z, _Bool bigimz); +void ewald3_2_sigma_long_Delta_recurrent(_Complex double *target, double *target_err, int maxn, _Complex double x, + int xbranch, _Complex double z, _Bool bigimz); /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /** This function always uses Taylor expansion in \a z. @@ -255,26 +259,26 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_ * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth * of the error. */ -void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x, - int xbranch, complex double z); +void ewald3_2_sigma_long_Delta_series(_Complex double *target, double *target_err, int maxn, _Complex double x, + int xbranch, _Complex double z); // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space /// The Ewald sum "self-interaction" term that appears in the lattice sums with zero (direct-space) Bravais lattice displacement. -int ewald3_sigma0(complex double *result, ///< Pointer to save the result (single complex double). +int ewald3_sigma0(_Complex double *result, ///< Pointer to save the result (single _Complex double). double *err, ///< Pointer to save the result error estimate (single double). const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init(). double eta, ///< Ewald parameter. - complex double wavenumber ///< Wavenumber of the background medium. + _Complex double wavenumber ///< Wavenumber of the background medium. ); /// Short-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{S}(\vect k,\vect s)\f$. int ewald3_sigma_short( - complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long. + _Complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long. double *target_sigmasr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`. const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init(). double eta, ///< Ewald parameter. - complex double wavenumber, ///< Wavenumber of the background medium. + _Complex double wavenumber, ///< Wavenumber of the background medium. /// Lattice dimensionality. /** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */ LatticeDimensionality latdim, @@ -285,7 +289,7 @@ int ewald3_sigma_short( * In such case, it is the responsibility of the caller to deallocate * the generator. */ - PGen *pgen_R, + struct PGen *pgen_R, /// Indicates whether pgen_R already generates shifted points. /** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(), * so the function assumes that the generated points correspond to the unshifted Bravais lattice, @@ -302,11 +306,11 @@ int ewald3_sigma_short( /// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$. int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim - complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long. + _Complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long. double *target_sigmalr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`. const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init(). double eta, ///< Ewald parameter. - complex double wavenumber, ///< Wavenumber of the background medium. + _Complex double wavenumber, ///< Wavenumber of the background medium. double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality). /// Lattice dimensionality. LatticeDimensionality latdim, @@ -317,7 +321,7 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep * In such case, it is the responsibility of the caller to deallocate * the generator. */ - PGen *pgen_K, + struct PGen *pgen_K, /// Indicates whether pgen_K already generates shifted points. /** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(), * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice, @@ -336,4 +340,8 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep extern int ewald_factor_ipow_l; // If nonzero, adds an additional factor \f$ i^{nm} \f$ to the Ewald sum result (for debubbing). extern int ewald_factor_ipow_m; + +#ifdef __cplusplus +} +#endif #endif //EWALD_H diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index c2ac92d..1688287 100644 --- a/qpms/ewaldsf.c +++ b/qpms/ewaldsf.c @@ -14,6 +14,7 @@ #include #include #include "tiny_inlines.h" +#include "qpms_error.h" // Some magic constants diff --git a/qpms/gaunt.h b/qpms/gaunt.h index 425c660..d6e2e81 100644 --- a/qpms/gaunt.h +++ b/qpms/gaunt.h @@ -3,6 +3,10 @@ */ #ifndef GAUNT_H #define GAUNT_H +#ifdef __cplusplus +extern "C" { +#endif + #include #define _GAUNT_H_MIN(x,y) (((x) > (y)) ? (y) : (x)) @@ -30,4 +34,7 @@ double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu); int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu); #endif //GAUNT_PRECOMPILED +#ifdef __cplusplus +} +#endif #endif //GAUNT_H diff --git a/qpms/groups.h b/qpms/groups.h index e78e462..36142d9 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -23,6 +23,9 @@ */ #ifndef QPMS_GROUPS_H #define QPMS_GROUPS_H +#ifdef __cplusplus +extern "C" { +#endif #include "qpms_types.h" #include @@ -35,7 +38,7 @@ struct qpms_finite_group_irrep_t { /** The r-th row, c-th column of the representation of the i'th element is retrieved as * m[i * dim * dim + r * dim + c] */ - complex double *m; + _Complex double *m; }; /// A point group with its irreducible representations and some metadata. @@ -93,4 +96,7 @@ qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *na extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL; extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL_G; +#ifdef __cplusplus +} +#endif #endif // QPMS_GROUPS_H diff --git a/qpms/indexing.h b/qpms/indexing.h index aab7478..e8ef012 100644 --- a/qpms/indexing.h +++ b/qpms/indexing.h @@ -26,6 +26,9 @@ */ #ifndef QPMS_INDEXING_H #define QPMS_INDEXING_H +#ifdef __cplusplus +extern "C" { +#endif #include "qpms_types.h" #include @@ -106,7 +109,7 @@ static const qpms_uvswfi_t QPMS_UI_L00 = 0; /** Returns a non-zero value if the u value is invalid. */ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) { - *t = u & 3; + *t = (qpms_vswf_type_t)(u & 3); qpms_y_sc_t y_sc = u / 4; qpms_y2mn_sc_p(y_sc, m, n); // Test validity @@ -119,7 +122,7 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, /** Does *not* allow for longitudinal waves. */ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u, qpms_vswf_type_t *t, qpms_y_t *y) { - *t = u & 3; + *t = (qpms_vswf_type_t)(u & 3); *y = u / 4 - 1; if (QPMS_UNLIKELY(*t == 0 || *t == 3)) return QPMS_ERROR; if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR; @@ -132,7 +135,7 @@ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u, */ static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u, qpms_vswf_type_t *t, qpms_y_t *y) { - *t = u & 3; + *t = (qpms_vswf_type_t)(u & 3); *y = u / 4 - 1; if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR; if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR; @@ -147,4 +150,7 @@ static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) { } +#ifdef __cplusplus +} +#endif #endif //QPMS_INDEXING_H diff --git a/qpms/kahansum.h b/qpms/kahansum.h index 80ce4a7..c2ef963 100644 --- a/qpms/kahansum.h +++ b/qpms/kahansum.h @@ -3,8 +3,9 @@ */ #ifndef KAHANSUM_H #define KAHANSUM_H - -#include +#ifdef __cplusplus +extern "C" { +#endif static inline void kahaninit(double * const sum, double * const compensation) { *sum = 0; @@ -19,16 +20,19 @@ static inline void kahanadd(double *sum, double *compensation, double input) { } -static inline void ckahaninit(complex double * const sum, complex double * const compensation) { +static inline void ckahaninit(_Complex double * const sum, _Complex double * const compensation) { *sum = 0; *compensation = 0; } -static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) { - complex double compensated_input = input - *compensation; - complex double nsum = *sum + compensated_input; +static inline void ckahanadd(_Complex double *sum, _Complex double *compensation, _Complex double input) { + _Complex double compensated_input = input - *compensation; + _Complex double nsum = *sum + compensated_input; *compensation = (nsum - *sum) - compensated_input; *sum = nsum; } +#ifdef __cplusplus +} +#endif #endif //KAHANSUM_H diff --git a/qpms/lattices.h b/qpms/lattices.h index d7daa3c..395a184 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -1,9 +1,15 @@ /*! \file lattices.h * \brief Lattice point generators and lattice vector analysis / transformation. * + * \bug Header file not C++ compatible. */ #ifndef LATTICES_H #define LATTICES_H +#ifdef __cplusplus // FIXME Not C++ compatible. Include "lattices_types.h" for minimal necessary enum decls. +extern "C" { +#endif + +#include "lattices_types.h" #include #include @@ -30,22 +36,6 @@ * */ -typedef enum LatticeDimensionality { - LAT1D = 1, - LAT2D = 2, - LAT3D = 4, - SPACE1D = 8, - SPACE2D = 16, - SPACE3D = 32, - LAT_1D_IN_3D = 33, - LAT_2D_IN_3D = 34, - LAT_3D_IN_3D = 40, - // special coordinate arrangements (indicating possible optimisations) - LAT_ZONLY = 64, - LAT_XYONLY = 128, - LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64 - LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128 -} LatticeDimensionality; inline static bool LatticeDimensionality_checkflags( LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) { @@ -945,4 +935,7 @@ int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxste int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r); void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g); +#ifdef __cplusplus +} +#endif #endif // LATTICES_H diff --git a/qpms/materials.c b/qpms/materials.c index 6b84949..5e5397e 100644 --- a/qpms/materials.c +++ b/qpms/materials.c @@ -5,6 +5,7 @@ #include #include "materials.h" #include "qpms_error.h" +#include #define SQ(x) ((x)*(x)) diff --git a/qpms/materials.h b/qpms/materials.h index 912e15d..c14e041 100644 --- a/qpms/materials.h +++ b/qpms/materials.h @@ -3,7 +3,12 @@ */ #ifndef QPMS_MATERIALS_H #define QPMS_MATERIALS_H +#ifdef __cplusplus +extern "C" { +#endif + #include "qpms_types.h" +#include #include #ifndef SPEED_OF_LIGHT @@ -20,36 +25,36 @@ typedef struct qpms_epsmu_generator_t { * qpms_permittivity_interpolator_epsmu_g(), * qpms_lorentzdrude_epsmu_g(). */ - qpms_epsmu_t (*function) (complex double omega, const void *params); + qpms_epsmu_t (*function) (_Complex double omega, const void *params); const void *params; } qpms_epsmu_generator_t; /// Convenience function for generating material properties at given frequency. static inline qpms_epsmu_t qpms_epsmu_generator_eval( - qpms_epsmu_generator_t gen, complex double omega) { + qpms_epsmu_generator_t gen, _Complex double omega) { return gen.function(omega, gen.params); } /// Constant optical property "generator" for qpms_epsmu_generator_t. -qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored. +qpms_epsmu_t qpms_epsmu_const_g(_Complex double omega, ///< Frequency ignored. const void *epsmu ///< Points to the qpms_epsmu_t to be returned. ); /// Gets refractive index of a material from its permeability and permittivity. /** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */ -static inline complex double qpms_refindex(qpms_epsmu_t em) { +static inline _Complex double qpms_refindex(qpms_epsmu_t em) { return csqrt(em.eps * em.mu); } /// Gets wave number \a k from angular frequency and material permeability and permittivity. /** \f[ k = \frac{n\omega}{c_0} = \frac{\omega\sqrt{\mu_r \varepsilon_r}}{c_0} \f] */ -static inline complex double qpms_wavenumber(complex double omega, qpms_epsmu_t em) { +static inline _Complex double qpms_wavenumber(_Complex double omega, qpms_epsmu_t em) { return qpms_refindex(em)*omega/SPEED_OF_LIGHT; } /// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity. /** \eta_r = \sqrt{\mu_r / \varepsilon_r} \f] */ -static inline complex double qpms_waveimpedance(qpms_epsmu_t em) { +static inline _Complex double qpms_waveimpedance(qpms_epsmu_t em) { return csqrt(em.mu / em.eps); } @@ -67,7 +72,7 @@ typedef struct qpms_ldparams_triple_t { * \f] */ typedef struct qpms_ldparams_t { - complex double eps_inf; ///< Permittivity at infinity. + _Complex double eps_inf; ///< Permittivity at infinity. double omega_p; ///< Plasma frequency. size_t n; ///< Number of "oscillators". qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters. @@ -86,14 +91,14 @@ extern const qpms_ldparams_t *const QPMS_LDPARAMS_PT; ///< Lorentz-Drude paramet extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten. /// Lorentz-Drude permittivity. -complex double qpms_lorentzdrude_eps(complex double omega, const qpms_ldparams_t *); +_Complex double qpms_lorentzdrude_eps(_Complex double omega, const qpms_ldparams_t *); /// Lorentz-Drude optical properties, with relative permeability set always to one. -qpms_epsmu_t qpms_lorentzdrude_epsmu(complex double omega, const qpms_ldparams_t *); +qpms_epsmu_t qpms_lorentzdrude_epsmu(_Complex double omega, const qpms_ldparams_t *); /// Lorentz-Drude optical properties, with relative permeability set always to one, compatible with qpms_epsmu_generator_t. qpms_epsmu_t qpms_lorentzdrude_epsmu_g( - complex double omega, + _Complex double omega, const void *ldparams ///< Lorentz-Drude parameters, in reality const qpms_ldparams_t *. ); @@ -125,14 +130,14 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml( ); /// Evaluates interpolated material permittivity at a given angular frequency. -complex double qpms_permittivity_interpolator_eps_at_omega( +_Complex double qpms_permittivity_interpolator_eps_at_omega( const qpms_permittivity_interpolator_t *interp, double omega_SI); /// Evaluates interpolated material permittivity at a given angular frequency, qpms_epsmu_generator_t compatible version. /** Permeability is always set to one. Imaginary part of omega is discarded. */ qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g( - complex double omega, ///< Angular frequency. The imaginary part is ignored! + _Complex double omega, ///< Angular frequency. The imaginary part is ignored! const void * interpolator ///< Interpolator of type qpms_permittivity_interpolator_t ); @@ -148,14 +153,17 @@ double qpms_permittivity_interpolator_omega_max( void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp); /// Relative permittivity from the Drude model. -static inline complex double qpms_drude_epsilon( - complex double eps_inf, ///< Relative permittivity "at infinity". - complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material. - complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material. - complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated. +static inline _Complex double qpms_drude_epsilon( + _Complex double eps_inf, ///< Relative permittivity "at infinity". + _Complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material. + _Complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material. + _Complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated. ) { return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p)); } +#ifdef __cplusplus +} +#endif #endif //QPMS_MATERIALS_H diff --git a/qpms/normalisation.h b/qpms/normalisation.h index f84f546..a282231 100644 --- a/qpms/normalisation.h +++ b/qpms/normalisation.h @@ -2,14 +2,18 @@ * \brief Convention-dependent coefficients for VSWFs. * * See also @ref qpms_normalisation_t and @ref vswf_conventions. + * + * \bug Header file not C++ compatible. */ #ifndef NORMALISATION_H #define NORMALISATION_H +#ifdef __cplusplus //FIXME not C++ compatible yet (enum bit operations) +extern "C" { +#endif #include "qpms_types.h" #include "qpms_error.h" #include -#include #include "indexing.h" #include "optim.h" @@ -36,8 +40,8 @@ static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qp * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.) */ -static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_normfactor(norm, l, m); +static inline _Complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_normfactor(norm, l, m); if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; @@ -51,8 +55,8 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * a `gsl_sf_legendre_*_e()` call. */ -static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m); +static inline _Complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m); return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; } @@ -62,8 +66,8 @@ static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t no * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.) */ -static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_normfactor(norm, l, m); +static inline _Complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_normfactor(norm, l, m); if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; @@ -77,14 +81,14 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * a `gsl_sf_legendre_*_e()` call. */ -static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m); +static inline _Complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m); return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; } /// Returns the factors of a electric basis VSWF divided by the factor of a magnetic VWFS of a given convention, compared to the reference one. -static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { +static inline _Complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { return qpms_normalisation_factor_N_noCS(norm, l, m) / qpms_normalisation_factor_M_noCS(norm, l, m); } @@ -95,8 +99,8 @@ static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * phase is already taken into account in a `gsl_sf_legendre_*_e()` call.) */ -static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_normfactor(norm, l, m); +static inline _Complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_normfactor(norm, l, m); if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I; if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac; @@ -109,13 +113,13 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * a `gsl_sf_legendre_*_e()` call. */ -static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { - complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m); +static inline _Complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { + _Complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m); return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; } /// Returns the factors of a basis VSWF of a given convention compared to the reference convention. -static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) { +static inline _Complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) { qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l; qpms_uvswfi2tmn(ui, &t, &m, &l); switch(t) { @@ -156,7 +160,7 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t * 0 \quad \mbox{if } m>0. \\ * \f] */ -static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) { +static inline _Complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) { switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { case 0: @@ -194,7 +198,7 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor * * */ -static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) { +static inline _Complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) { if(m==0) return 0; switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { @@ -283,4 +287,7 @@ static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t } #endif +#ifdef __cplusplus +} +#endif #endif //NORMALISATION_H diff --git a/qpms/parsing.h b/qpms/parsing.h index 3f60daa..aba1f30 100644 --- a/qpms/parsing.h +++ b/qpms/parsing.h @@ -3,6 +3,9 @@ */ #ifndef QPMS_PARSING_H #define QPMS_PARSING_H +#ifdef __cplusplus +extern "C" { +#endif #include @@ -57,4 +60,7 @@ size_t qpms_parse_doubles_fromfile( const char *filepath //< File to read from, or NULL, "", "-" to read from stdin. ); +#ifdef __cplusplus +} +#endif #endif // QPMS_PARSING_H diff --git a/qpms/pointgroups.c b/qpms/pointgroups.c index 0125655..a2b1cda 100644 --- a/qpms/pointgroups.c +++ b/qpms/pointgroups.c @@ -1,6 +1,7 @@ #include "pointgroups.h" #include #include +#include #include double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL; diff --git a/qpms/pointgroups.h b/qpms/pointgroups.h index 0fc92fe..2ddcac2 100644 --- a/qpms/pointgroups.h +++ b/qpms/pointgroups.h @@ -3,6 +3,9 @@ */ #ifndef POINTGROUPS_H #define POINTGROUPS_H +#ifdef __cplusplus +extern "C" { +#endif #include "qpms_error.h" #include "quaternions.h" @@ -18,9 +21,9 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) { case QPMS_PGS_DN: case QPMS_PGS_DND: case QPMS_PGS_DNH: - return true; + return 1; default: - return false; + return 0; } } @@ -81,4 +84,7 @@ _Bool qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b); +#ifdef __cplusplus +} +#endif #endif //POINTGROUPS_H diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 447041b..64121be 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -493,6 +493,7 @@ cdef extern from "tmatrices.h": qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, cdouble omega, const void *params) qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, cdouble omega, const void *params) qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, cdouble omega, const void *params) + qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, cdouble omega, const void *params) struct qpms_tmatrix_generator_sphere_param_t: qpms_epsmu_generator_t outside qpms_epsmu_generator_t inside diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index f554035..bb30c5c 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -5,6 +5,9 @@ */ #ifndef QPMS_ERROR_H #define QPMS_ERROR_H +#ifdef __cplusplus +extern "C" { +#endif #include "optim.h" @@ -266,4 +269,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); }\ } +#ifdef __cplusplus +} +#endif #endif diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 3c07a91..213a38d 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -3,6 +3,10 @@ */ #ifndef QPMS_SPECFUNC_H #define QPMS_SPECFUNC_H +#ifdef __cplusplus +extern "C" { +#endif + #include "qpms_types.h" #include @@ -11,25 +15,25 @@ ******************************************************************************/ // TODO unify types -qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array); +qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, _Complex double x, _Complex double *result_array); typedef struct { qpms_l_t lMax; double *akn; // coefficients as in DLMF 10.49.1 - //complex double *bkn; // coefficients of the derivatives + //_Complex double *bkn; // coefficients of the derivatives } qpms_sbessel_calculator_t; qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(void); void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c); qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax, - double x, complex double *result_array); + double x, _Complex double *result_array); -complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x); +_Complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, _Complex double x); qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax, - complex double x, complex double *result_array); + _Complex double x, _Complex double *result_array); /****************************************************************************** @@ -107,4 +111,7 @@ double qpms_legendre0(int m, int n); // Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2) double qpms_legendred0(int m, int n); +#ifdef __cplusplus +} +#endif #endif // QPMS_SPECFUNC_H diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 7b7bd3b..209d07e 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -3,10 +3,13 @@ */ #ifndef QPMS_TYPES_H #define QPMS_TYPES_H -#include -#include +#ifdef __cplusplus +extern "C" { +#endif + #include #include +#include #ifndef M_PI_2 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L) @@ -191,7 +194,7 @@ typedef struct cart3_t { /// 3D complex (actually 6D) coordinates. See also vectors.h. typedef struct ccart3_t { - complex double x, y, z; + _Complex double x, y, z; } ccart3_t; /// 3D complex vector pair (represents the E, H fields). @@ -212,13 +215,13 @@ typedef struct sph_t { /// Spherical coordinates with complex radial component. See also vectors.h. typedef struct csph_t { // Do I really need this??? - complex double r; + _Complex double r; double theta, phi; } csph_t; /// 3D complex vector components in local spherical basis. See also vectors.h. typedef struct csphvec_t { - complex double rc, thetac, phic; + _Complex double rc, thetac, phic; } csphvec_t; /// 2D polar coordinates. See also vectors.h. @@ -261,7 +264,7 @@ typedef enum { * See quaternions.h for "methods". */ typedef struct qpms_quat_t { - complex double a, b; + _Complex double a, b; } qpms_quat_t; /// Quaternion type as four doubles. @@ -346,7 +349,7 @@ typedef struct qpms_tmatrix_t { * there is also one with order \a −m. */ const qpms_vswf_set_spec_t *spec; - complex double *m; ///< Matrix elements in row-major order. + _Complex double *m; ///< Matrix elements in row-major order. bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free() } qpms_tmatrix_t; @@ -422,8 +425,8 @@ typedef struct qpms_abstract_tmatrix_t { /// A type holding electric permittivity and magnetic permeability of a material. typedef struct qpms_epsmu_t { - complex double eps; ///< Relative permittivity. - complex double mu; ///< Relative permeability. + _Complex double eps; ///< Relative permittivity. + _Complex double mu; ///< Relative permeability. } qpms_epsmu_t; struct qpms_tolerance_spec_t; // See tolerances.h @@ -436,4 +439,8 @@ typedef enum { } qpms_ewald_part; #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l)) + +#ifdef __cplusplus +} +#endif #endif // QPMS_TYPES diff --git a/qpms/qpmsblas.h b/qpms/qpmsblas.h index c32282b..3981b82 100644 --- a/qpms/qpmsblas.h +++ b/qpms/qpmsblas.h @@ -14,6 +14,10 @@ */ #ifndef QPMSBLAS_H #define QPMSBLAS_H +#ifdef __cplusplus +extern "C" { +#endif + #define QPMS_BLAS_INDEX_T long long int #ifndef CBLAS_H @@ -31,4 +35,7 @@ void qpms_zgemm(CBLAS_LAYOUT Order, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE Tran const _Complex double *B, const QPMS_BLAS_INDEX_T ldb, const _Complex double *beta, _Complex double *C, const QPMS_BLAS_INDEX_T ldc); +#ifdef __cplusplus +} +#endif #endif //QPMSBLAS_H diff --git a/qpms/quaternions.h b/qpms/quaternions.h index eed1ce8..474b06b 100644 --- a/qpms/quaternions.h +++ b/qpms/quaternions.h @@ -3,6 +3,9 @@ */ #ifndef QPMS_WIGNER_H #define QPMS_WIGNER_H +#ifdef __cplusplus +extern "C" { +#endif #include "qpms_types.h" #include "vectors.h" @@ -93,7 +96,7 @@ static inline double qpms_quat_norm(const qpms_quat_t q) { } /// Test approximate equality of quaternions. -static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) { +static inline _Bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) { return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol; } @@ -119,7 +122,7 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) { } /// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$. -static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) { +static inline _Bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) { return qpms_quat_norm(qpms_quat_sub( qpms_quat_standardise(p, atol), qpms_quat_standardise(q, atol))) <= atol; @@ -202,17 +205,17 @@ static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) { * The D matrix are calculated using formulae (3), (4), (6), (7) from * http://moble.github.io/spherical_functions/WignerDMatrices.html */ -complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l, +_Complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l, qpms_m_t mp, qpms_m_t m); /// A VSWF representation element of the O(3) group. /** * TODO more doc. */ -complex double qpms_vswf_irot_elem_from_irot3( +_Complex double qpms_vswf_irot_elem_from_irot3( const qpms_irot3_t q, ///< The O(3) element in the quaternion representation. qpms_l_t l, qpms_m_t mp, qpms_m_t m, - bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise. + _Bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise. ); @@ -251,7 +254,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) { } /// Test approximate equality of irot3. -static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) { +static inline _Bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) { return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det; } @@ -303,4 +306,7 @@ static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) { return qpms_irot3_zrot_angle(2 * M_PI * k / N); } +#ifdef __cplusplus +} +#endif #endif //QPMS_WIGNER_H diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 6cc0e88..c3affc2 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -24,6 +24,7 @@ #include "tolerances.h" #include "beyn.h" #include "tiny_inlines.h" +#include "lattices.h" #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #include "qpmsblas.h" diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 6f1ee65..eb3337d 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -10,6 +10,10 @@ */ #ifndef QPMS_SCATSYSTEM_H #define QPMS_SCATSYSTEM_H +#ifdef __cplusplus +extern "C" { +#endif + #include "qpms_types.h" #include "vswf.h" #include "tmatrices.h" @@ -107,7 +111,7 @@ typedef struct qpms_ss_orbit_type_t { * * TODO doc. */ - complex double *irbases; + _Complex double *irbases; /// TODO doc. size_t instance_count; /// Cumulative sum of the preceding ot->siza * ot->instance_count; @@ -243,9 +247,9 @@ typedef struct qpms_scatsys_at_omega_t { * i.e in the order corresponding to \a ss->tm. */ qpms_tmatrix_t **tm; - complex double omega; ///< Angular frequency + _Complex double omega; ///< Angular frequency qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency - complex double wavenumber; ///< Background medium wave number + _Complex double wavenumber; ///< Background medium wave number } qpms_scatsys_at_omega_t; @@ -291,7 +295,7 @@ typedef struct qpms_scatsys_at_omega_t { * */ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym, - complex double omega, const struct qpms_tolerance_spec_t *tol); + _Complex double omega, const struct qpms_tolerance_spec_t *tol); /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); @@ -303,50 +307,50 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw); /// Evaluates scattering system T-matrices at a given frequency. /** Free the result using qpms_scatsys_at_omega_free() when done. */ qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, - complex double omega); + _Complex double omega); /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. * * TODO doc about shape etc. */ -complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U, +_Complex double *qpms_scatsys_irrep_transform_matrix(_Complex double *target_U, const qpms_scatsys_t *ss, qpms_iri_t iri); /// Projects a "big" matrix onto an irrep (slow reference implementation). /** TODO doc */ -complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed, - const complex double *orig_full, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_pack_matrix_stupid(_Complex double *target_packed, + const _Complex double *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri); /// Transforms a big "packed" matrix into the full basis (slow reference implementation). /** TODO doc */ -complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full, - const complex double *orig_packed, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_unpack_matrix_stupid(_Complex double *target_full, + const _Complex double *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bool add); /// Projects a "big" matrix onto an irrep. /** TODO doc */ -complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, - const complex double *orig_full, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_pack_matrix(_Complex double *target_packed, + const _Complex double *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri); /// Transforms a big "packed" matrix into the full basis. /** TODO doc */ -complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full, - const complex double *orig_packed, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_unpack_matrix(_Complex double *target_full, + const _Complex double *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bool add); /// Projects a "big" vector onto an irrep. /** TODO doc */ -complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed, - const complex double *orig_full, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_pack_vector(_Complex double *target_packed, + const _Complex double *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri); /// Transforms a big "packed" vector into the full basis. /** TODO doc */ -complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, - const complex double *orig_packed, const qpms_scatsys_t *ss, +_Complex double *qpms_scatsys_irrep_unpack_vector(_Complex double *target_full, + const _Complex double *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bool add); /// Global translation matrix. @@ -354,31 +358,31 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, * The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves). * This may change in the future. */ -complex double *qpms_scatsys_build_translation_matrix_full( +_Complex double *qpms_scatsys_build_translation_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_t *ss, - complex double k ///< Wave number to use in the translation matrix. + _Complex double k ///< Wave number to use in the translation matrix. ); /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. /** * \returns \a target on success, NULL on error. */ -complex double *qpms_scatsysw_build_modeproblem_matrix_full( +_Complex double *qpms_scatsysw_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_t *ssw ); /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type. /** Might be useful for evaluation of cross sections and testing. */ -complex double *qpms_scatsys_build_translation_matrix_e_full( +_Complex double *qpms_scatsys_build_translation_matrix_e_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_t *ss, - complex double k, ///< Wave number to use in the translation matrix. + _Complex double k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ); @@ -387,12 +391,12 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( * The diagonal (particle self-) blocks are currently filled with zeros. * This may change in the future. */ -complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( +_Complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k, ///< Wave number to use in the translation matrix. + _Complex double k, ///< Wave number to use in the translation matrix. qpms_bessel_t J ); @@ -400,9 +404,9 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( /** * \returns \a target on success, NULL on error. */ -complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( +_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); @@ -410,9 +414,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( /** * \returns \a target on success, NULL on error. */ -complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( +_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); @@ -420,9 +424,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( /** * \returns \a target on success, NULL on error. */ -complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( +_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym ); @@ -435,7 +439,7 @@ typedef struct qpms_ss_LU { bool full; ///< true if full matrix; false if irrep-packed. qpms_iri_t iri; ///< Irrep index if `full == false`. /// LU decomposition array. - complex double *a; + _Complex double *a; /// Pivot index array, size at least max(1,min(m, n)). int *ipiv; } qpms_ss_LU; @@ -443,14 +447,14 @@ void qpms_ss_LU_free(qpms_ss_LU); /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only. qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( - complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). + _Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw ); /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( - complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). + _Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri @@ -458,7 +462,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( - complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. + _Complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems. const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL. @@ -466,16 +470,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise( - complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. + _Complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ); /// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$. -complex double *qpms_scatsys_scatter_solve( - complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated. - const complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`). +_Complex double *qpms_scatsys_scatter_solve( + _Complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated. + const _Complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`). qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data. ); @@ -495,33 +499,33 @@ typedef struct qpms_scatsys_at_omega_k_t { /** * \returns \a target on success, NULL on error. */ -complex double *qpms_scatsyswk_build_modeproblem_matrix_full( +_Complex double *qpms_scatsyswk_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_k_t *sswk ); /// Global translation matrix. -complex double *qpms_scatsys_periodic_build_translation_matrix_full( +_Complex double *qpms_scatsys_periodic_build_translation_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_t *ss, - complex double wavenumber, ///< Wave number to use in the translation matrix. + _Complex double wavenumber, ///< Wave number to use in the translation matrix. const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates. double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss. ); /// Global translation matrix. -complex double *qpms_scatsyswk_build_translation_matrix_full( +_Complex double *qpms_scatsyswk_build_translation_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, + _Complex double *target, const qpms_scatsys_at_omega_k_t *sswk ); /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only. qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU( - complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). + _Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_k_t *sswk ); @@ -539,7 +543,7 @@ struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes( const qpms_scatsys_t *ss, /// Wavevector in cartesian coordinates (must lie in the lattice plane). const double k[3], - complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. + _Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for. double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for. size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number) @@ -561,12 +565,12 @@ struct qpms_finite_group_t; /// Constructs a "full matrix action" of a point group element for an orbit type. /** TODO detailed doc */ -complex double *qpms_orbit_action_matrix( +_Complex double *qpms_orbit_action_matrix( /// Target array. If NULL, a new one is allocated. /** The size of the array is (orbit->size * bspec->n)**2 * (it makes sense to assume all the T-matrices share their spec). */ - complex double *target, + _Complex double *target, /// The orbit (type). const qpms_ss_orbit_type_t *orbit, /// Base spec of the t-matrices (we don't know it from orbit, as it has @@ -579,12 +583,12 @@ complex double *qpms_orbit_action_matrix( /// Constructs a dense matrix representation of a irrep projector for an orbit type. /** TODO detailed doc */ -complex double *qpms_orbit_irrep_projector_matrix( +_Complex double *qpms_orbit_irrep_projector_matrix( /// Target array. If NULL, a new one is allocated. /** The size of the array is (orbit->size * bspec->n)**2 * (it makes sense to assume all the T-matrices share their spec). */ - complex double *target, + _Complex double *target, /// The orbit (type). const qpms_ss_orbit_type_t *orbit, /// Base spec of the t-matrices (we don't know it from orbit, as it has @@ -596,14 +600,14 @@ complex double *qpms_orbit_irrep_projector_matrix( const qpms_iri_t iri); /// TODO DOC!!!!! -complex double *qpms_orbit_irrep_basis( +_Complex double *qpms_orbit_irrep_basis( /// Here theh size of theh basis shall be saved, size_t *basis_size, /// Target array. If NULL, a new one is allocated. /** The size of the array is basis_size * (orbit->size * bspec->n) * (it makes sense to assume all the T-matrices share their spec). */ - complex double *target, + _Complex double *target, /// The orbit (type). const qpms_ss_orbit_type_t *orbit, /// Base spec of the t-matrices (we don't know it from orbit, as it has @@ -618,10 +622,10 @@ complex double *qpms_orbit_irrep_basis( /// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points. /** TODO detailed doc! * \returns target_full if target_full was not NULL, otherwise the newly allocated array. */ -complex double *qpms_scatsys_incident_field_vector_full( +_Complex double *qpms_scatsys_incident_field_vector_full( /// Target array. If NULL, a new one is allocated. /** The length of the array is ss->fecv_size. */ - complex double *target_full, + _Complex double *target_full, const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, const void *args, ///< Pointer passed as the last argument to (*field_at_point)() @@ -629,9 +633,9 @@ complex double *qpms_scatsys_incident_field_vector_full( ); /// Applies T-matrices onto an incident field vector in the full basis. -complex double *qpms_scatsysw_apply_Tmatrices_full( - complex double *target_full, /// Target vector array. If NULL, a new one is allocated. - const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL. +_Complex double *qpms_scatsysw_apply_Tmatrices_full( + _Complex double *target_full, /// Target vector array. If NULL, a new one is allocated. + const _Complex double *inc_full, /// Incident field coefficient vector. Must not be NULL. const qpms_scatsys_at_omega_t *ssw ); @@ -650,7 +654,7 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes( const qpms_scatsys_t *ss, /// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system. qpms_iri_t iri, - complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. + _Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for. double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for. size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number) @@ -673,7 +677,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes( const qpms_scatsys_t *ss, double eta, ///< Ewald sum parameter const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system. - complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. + _Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for. double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for. double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for. size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number) @@ -687,10 +691,10 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes( #if 0 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points. /** TODO detailed doc! */ -complex double *qpms_scatsys_incident_field_vector_irrep_packed( +_Complex double *qpms_scatsys_incident_field_vector_irrep_packed( /// Target array. If NULL, a new one is allocated. /** The length of the array is ss->fecv_size. */ - complex double *target_full, + _Complex double *target_full, const qpms_scatsys_t *ss, const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym. qpms_incfield_t field_at_point, @@ -715,8 +719,8 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( ccart3_t qpms_scatsys_scattered_E( const qpms_scatsys_t *ss, qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - complex double wavenumber, ///< Wavenumber of the background medium. - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + _Complex double wavenumber, ///< Wavenumber of the background medium. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -734,7 +738,7 @@ ccart3_t qpms_scatsys_scattered_E( ccart3_t qpms_scatsysw_scattered_E( const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -753,7 +757,7 @@ qpms_errno_t qpms_scatsys_scattered_field_basis( ccart3_t *target, ///< Target array of length \a ss->fecv_size const qpms_scatsys_t *ss, qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - complex double wavenumber, ///< Wavenumber of the background medium. + _Complex double wavenumber, ///< Wavenumber of the background medium. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -773,7 +777,7 @@ qpms_errno_t qpms_scatsys_scattered_field_basis_pi( const qpms_scatsys_t *ss, qpms_ss_pi_t pi, ///< Particle index qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - complex double wavenumber, ///< Wavenumber of the background medium + _Complex double wavenumber, ///< Wavenumber of the background medium cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -832,8 +836,8 @@ qpms_errno_t qpms_scatsysw_scattered_field_basis_pi( ccart3_t qpms_scatsys_scattered_E__alt( const qpms_scatsys_t *ss, qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - complex double wavenumber, ///< Wavenumber of the background medium. - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + _Complex double wavenumber, ///< Wavenumber of the background medium. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -849,7 +853,7 @@ ccart3_t qpms_scatsys_scattered_E__alt( ccart3_t qpms_scatsysw_scattered_E__alt( const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); @@ -868,14 +872,14 @@ ccart3_t qpms_scatsysw_scattered_E__alt( ccart3_t qpms_scatsyswk_scattered_E( const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported). - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); ccart3_t qpms_scatsyswk_scattered_E_e( const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported). - const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. + const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated. qpms_ewald_part parts ); @@ -903,7 +907,7 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis_e( /// Adjusted Ewadl parameter to avoid high-frequency breakdown. // TODO DOC -double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]); +double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, _Complex double wavenumber, const double wavevector[3]); #if 0 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) @@ -913,10 +917,13 @@ double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, */ ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, ///< Irreducible representation - const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri. + const _Complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri. cart3_t where, ///< Evaluation point. - complex double k ///< Wave number. + _Complex double k ///< Wave number. ); #endif +#ifdef __cplusplus +} +#endif #endif //QPMS_SCATSYSTEM_H diff --git a/qpms/scatsystem_dbg.h b/qpms/scatsystem_dbg.h index a39707d..e7ad0dc 100644 --- a/qpms/scatsystem_dbg.h +++ b/qpms/scatsystem_dbg.h @@ -1,7 +1,7 @@ #include "scatsystem.h" qpms_errno_t qpms_scatsyswk_test_sswf_basis_e( - complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size + _Complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted). cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated. diff --git a/qpms/symmetries.h b/qpms/symmetries.h index 97e15cd..4ef1f4d 100644 --- a/qpms/symmetries.h +++ b/qpms/symmetries.h @@ -21,40 +21,44 @@ */ #ifndef SYMMETRIES_H #define SYMMETRIES_H +#ifdef __cplusplus +extern "C" { +#endif + #include "qpms_types.h" #include /// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring). -complex double *qpms_zflip_uvswi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_zflip_uvswi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec); /// Dense matrix representation of the y coordinate sign flip operation (xz-plane mirroring). -complex double *qpms_yflip_uvswi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_yflip_uvswi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec); /// Dense matrix representation of the x coordinate sign flip operation (yz-plane mirroring). -complex double *qpms_xflip_uvswi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_xflip_uvswi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec); /// Dense matrix representation of a rotation around the z-axis -complex double *qpms_zrot_uvswi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_zrot_uvswi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec, double phi ///< Rotation angle ); /// Dense matrix representation of a "rational" rotation around the z-axis /** Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$. */ -complex double *qpms_zrot_rational_uvswi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_zrot_rational_uvswi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec, int N, int w ); /// Dense matrix (uvswi-indexed) representation of any O(3) transformation. -complex double *qpms_irot3_uvswfi_dense( - complex double *target, ///< If NULL, a new array is allocated. +_Complex double *qpms_irot3_uvswfi_dense( + _Complex double *target, ///< If NULL, a new array is allocated. const qpms_vswf_set_spec_t *bspec, const qpms_irot3_t transf); @@ -74,5 +78,9 @@ size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol); * Works on real and imaginary parts separately. * TODO doc. */ -size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol); +size_t qpms_czero_roundoff_clean(_Complex double *arr, size_t nmemb, double atol); + +#ifdef __cplusplus +} +#endif #endif // SYMMETRIES_H diff --git a/qpms/tiny_inlines.h b/qpms/tiny_inlines.h index ebde555..71ab686 100644 --- a/qpms/tiny_inlines.h +++ b/qpms/tiny_inlines.h @@ -4,6 +4,7 @@ #ifndef TINY_INLINES_H #define TINY_INLINES_H #include +#include static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; } @@ -19,8 +20,8 @@ static inline int min1pow_m_neg(int m) { #if 0 #ifdef __GSL_SF_LEGENDRE_H__ -static inline complex double -spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, complex double exp_imf) { +static inline _Complex double +spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, _Complex double exp_imf) { return; } @@ -28,9 +29,9 @@ spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, #endif // this has shitty precision: -// static inline complex double ipow(int x) { return cpow(I, x); } +// static inline _Complex double ipow(int x) { return cpow(I, x); } -static inline complex double ipow(int x) { +static inline _Complex double ipow(int x) { x = ((x % 4) + 4) % 4; switch(x) { case 0: diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 5ae6beb..7ac501f 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -988,6 +988,12 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, return QPMS_SUCCESS; } +qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, + complex double omega_ignored, const void *params_ignored) { + memset(t->m, 0, SQ(t->spec->n)*sizeof(*(t->m))); + return QPMS_SUCCESS; +} + void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { switch(f->typ) { case QPMS_TMATRIX_OPERATION_NOOP: diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index f9e3fd8..ee95d05 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -3,6 +3,10 @@ */ #ifndef TMATRICES_H #define TMATRICES_H +#ifdef __cplusplus +extern "C" { +#endif + // #include "qpms_types.h" // included via materials.h // #include // included via materials.h #include "materials.h" @@ -14,7 +18,7 @@ struct qpms_finite_group_t; typedef struct qpms_finite_group_t qpms_finite_group_t; /// Returns a pointer to the beginning of the T-matrix row \a rowno. -static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){ +static inline _Complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){ return t->m + (t->spec->n * rowno); } @@ -42,7 +46,7 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t); * This function actually checks for identical vswf specs. * TODO define constants with "default" atol, rtol for this function. */ -bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, +_Bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, const double rtol, const double atol); /// Creates a T-matrix from another matrix and a symmetry operation. @@ -53,7 +57,7 @@ bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, */ qpms_tmatrix_t *qpms_tmatrix_apply_symop( const qpms_tmatrix_t *T, ///< the original T-matrix - const complex double *M ///< the symmetry op matrix in row-major format + const _Complex double *M ///< the symmetry op matrix in row-major format ); /// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data. @@ -64,7 +68,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop( */ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( qpms_tmatrix_t *T, ///< the original T-matrix - const complex double *M ///< the symmetry op matrix in row-major format + const _Complex double *M ///< the symmetry op matrix in row-major format ); /// Symmetrizes a T-matrix with an involution symmetry operation. @@ -75,7 +79,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( */ qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution( const qpms_tmatrix_t *T, ///< the original T-matrix - const complex double *M ///< the symmetry op matrix in row-major format + const _Complex double *M ///< the symmetry op matrix in row-major format ); /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix. @@ -105,7 +109,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N( */ qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace( qpms_tmatrix_t *T, ///< the original T-matrix - const complex double *M ///< the symmetry op matrix in row-major format + const _Complex double *M ///< the symmetry op matrix in row-major format ); /// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one. @@ -149,7 +153,7 @@ qpms_errno_t qpms_load_scuff_tmatrix( double **freqs_su, ///< Frequencies in SCUFF units (optional). /// The resulting T-matrices (optional). qpms_tmatrix_t **tmatrices_array, - complex double **tmdata ///< The T-matrices raw contents + _Complex double **tmdata ///< The T-matrices raw contents ); /// Tells whether qpms_load_scuff_tmatrix should crash if fopen() fails. @@ -161,7 +165,7 @@ qpms_errno_t qpms_load_scuff_tmatrix( * This is desirable e.g. when used in Python (so that proper exception can * be thrown). */ -extern bool qpms_load_scuff_tmatrix_crash_on_failure; +extern _Bool qpms_load_scuff_tmatrix_crash_on_failure; /// Loads a scuff-tmatrix generated file. /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a @@ -189,7 +193,7 @@ qpms_errno_t qpms_read_scuff_tmatrix( * is accessed via * (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci]. */ - complex double ** tmdata + _Complex double ** tmdata ); /// In-place application of point group elements on raw T-matrix data. @@ -200,7 +204,7 @@ qpms_errno_t qpms_read_scuff_tmatrix( * TODO more doc. */ qpms_errno_t qpms_symmetrise_tmdata_irot3arr( - complex double *tmdata, const size_t tmcount, + _Complex double *tmdata, const size_t tmcount, const qpms_vswf_set_spec_t *bspec, size_t n_symops, const qpms_irot3_t *symops @@ -213,7 +217,7 @@ qpms_errno_t qpms_symmetrise_tmdata_irot3arr( * TODO more doc. */ qpms_errno_t qpms_symmetrise_tmdata_finite_group( - complex double *tmdata, const size_t tmcount, + _Complex double *tmdata, const size_t tmcount, const qpms_vswf_set_spec_t *bspec, const qpms_finite_group_t *pointgroup ); @@ -242,9 +246,9 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( ); /// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$. -complex double *qpms_apply_tmatrix( - complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated. - const complex double *a, ///< Incident field coefficient array of size T->spec->n. +_Complex double *qpms_apply_tmatrix( + _Complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated. + const _Complex double *a, ///< Incident field coefficient array of size T->spec->n. const qpms_tmatrix_t *T ///< T-matrix \a T to apply. ); @@ -254,10 +258,11 @@ complex double *qpms_apply_tmatrix( * qpms_tmatrix_generator_interpolator() * qpms_tmatrix_generator_sphere() * qpms_tmatrix_generator_constant() + * qpms_tmatrix_generator_zero() */ typedef struct qpms_tmatrix_generator_t { qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill. - complex double omega, ///< Angular frequency. + _Complex double omega, ///< Angular frequency. const void *params ///< Implementation dependent parameters. ); const void *params; ///< Parameter pointer passed to the function. @@ -267,7 +272,7 @@ typedef struct qpms_tmatrix_generator_t { qpms_tmatrix_t *qpms_tmatrix_init_from_generator( const qpms_vswf_set_spec_t *bspec, qpms_tmatrix_generator_t gen, - complex double omega); + _Complex double omega); /// Implementation of qpms_matrix_generator_t that just copies a constant matrix. @@ -275,11 +280,20 @@ qpms_tmatrix_t *qpms_tmatrix_init_from_generator( * the same base spec. */ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, - complex double omega, + _Complex double omega_ignored, /// Source T-matrix, real type is (const qpms_tmatrix_t*). const void *tmatrix_orig ); +/// Dummy implementation of qpms_tmatrix_generator_t, yielding a zero matrix. +/** + * This effectively represents a part of the background medium without + * any scatterer present. + */ +qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, + _Complex double omega_ignored, + const void *params_ignored); + /* Fuck this, include the whole typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct typedef struct gsl_interp_type gsl_interp_type; @@ -295,7 +309,7 @@ extern const gsl_interp_type * gsl_interp_steffen; // struct gsl_interp_accel; // use if lookup proves to be too slow typedef struct qpms_tmatrix_interpolator_t { const qpms_vswf_set_spec_t *bspec; - //bool owns_bspec; + //_Bool owns_bspec; gsl_spline **splines_real; ///< There will be a spline object for each nonzero element gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element // gsl_interp_accel **accel_real; @@ -317,7 +331,7 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided. const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers! const gsl_interp_type *iptype - //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix. + //, _Bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix. /*, ...? */); @@ -325,7 +339,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num /** As in qpms_tmatrix_interpolator_eval(), the imaginary part of frequency is discarded! */ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matrix to fill. - complex double omega, ///< Angular frequency. + _Complex double omega, ///< Angular frequency. const void *interpolator ///< Parameter of type qpms_tmatrix_interpolator_t *. ); @@ -342,14 +356,14 @@ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matri * * TODO better doc. */ -complex double *qpms_mie_coefficients_reflection( - complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated. +_Complex double *qpms_mie_coefficients_reflection( + _Complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated. const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated. double a, ///< Radius of the sphere. - complex double k_i, ///< Wave number of the internal material of the sphere. - complex double k_e, ///< Wave number of the surrounding medium. - complex double mu_i, ///< Relative permeability of the sphere material. - complex double mu_e, ///< Relative permeability of the surrounding medium. + _Complex double k_i, ///< Wave number of the internal material of the sphere. + _Complex double k_e, ///< Wave number of the surrounding medium. + _Complex double mu_i, ///< Relative permeability of the sphere material. + _Complex double mu_e, ///< Relative permeability of the surrounding medium. qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR. qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS. ); @@ -357,10 +371,10 @@ complex double *qpms_mie_coefficients_reflection( /// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory. qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. double a, ///< Radius of the sphere. - complex double k_i, ///< Wave number of the internal material of the sphere. - complex double k_e, ///< Wave number of the surrounding medium. - complex double mu_i, ///< Relative permeability of the sphere material. - complex double mu_e ///< Relative permeability of the surrounding medium. + _Complex double k_i, ///< Wave number of the internal material of the sphere. + _Complex double k_e, ///< Wave number of the surrounding medium. + _Complex double mu_i, ///< Relative permeability of the sphere material. + _Complex double mu_e ///< Relative permeability of the surrounding medium. ); /// Parameter structure for qpms_tmatrix_generator_sphere(). @@ -372,7 +386,7 @@ typedef struct qpms_tmatrix_generator_sphere_param_t { /// T-matrix generator for spherical particles using Lorentz-Mie solution. qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, - complex double omega, + _Complex double omega, const void *params ///< Of type qpms_tmatrix_generator_sphere_param_t. ); @@ -380,10 +394,10 @@ qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, static inline qpms_tmatrix_t *qpms_tmatrix_spherical( const qpms_vswf_set_spec_t *bspec, double a, ///< Radius of the sphere. - complex double k_i, ///< Wave number of the internal material of the sphere. - complex double k_e, ///< Wave number of the surrounding medium. - complex double mu_i, ///< Relative permeability of the sphere material. - complex double mu_e ///< Relative permeability of the surrounding medium. + _Complex double k_i, ///< Wave number of the internal material of the sphere. + _Complex double k_e, ///< Wave number of the surrounding medium. + _Complex double mu_i, ///< Relative permeability of the sphere material. + _Complex double mu_e ///< Relative permeability of the surrounding medium. ) { qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e); @@ -395,8 +409,8 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill( qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. double a, ///< Radius of the sphere. double omega, ///< Angular frequency. - complex double epsilon_fg, ///< Relative permittivity of the sphere. - complex double epsilon_bg ///< Relative permittivity of the background medium. + _Complex double epsilon_fg, ///< Relative permittivity of the sphere. + _Complex double epsilon_bg ///< Relative permittivity of the background medium. ); /// Convenience function to calculate T-matrix of a non-magnetic spherical particle using the permittivity values. @@ -404,8 +418,8 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0( const qpms_vswf_set_spec_t *bspec, double a, ///< Radius of the sphere. double omega, ///< Angular frequency. - complex double epsilon_fg, ///< Relative permittivity of the sphere. - complex double epsilon_bg ///< Relative permittivity of the background medium. + _Complex double epsilon_fg, ///< Relative permittivity of the sphere. + _Complex double epsilon_bg ///< Relative permittivity of the background medium. ) { qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg); @@ -460,7 +474,7 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta, */ qpms_errno_t qpms_tmatrix_axialsym_fill( qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. - complex double omega, ///< Angular frequency. + _Complex double omega, ///< Angular frequency. qpms_epsmu_t outside, ///< Optical properties of the outside medium. qpms_epsmu_t inside, ///< Optical properties of the particle's material. qpms_arc_function_t shape, ///< Particle surface parametrisation. @@ -474,7 +488,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill( /// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry. static inline qpms_tmatrix_t *qpms_tmatrix_axialsym( const qpms_vswf_set_spec_t *bspec, - complex double omega, ///< Angular frequency. + _Complex double omega, ///< Angular frequency. qpms_epsmu_t outside, ///< Optical properties of the outside medium. qpms_epsmu_t inside, ///< Optical properties of the particle's material. qpms_arc_function_t shape, ///< Particle surface parametrisation. @@ -501,22 +515,22 @@ typedef struct qpms_tmatrix_generator_axialsym_param_t { /// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to fill. - complex double omega, ///< Angular frequency. + _Complex double omega, ///< Angular frequency. const void *params ///< Parameters of type qpms_tmatrix_generator_axialsym_param_t. ); /// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful for debugging). -qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target, - complex double omega, +qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(_Complex double *target, + _Complex double omega, const qpms_tmatrix_generator_axialsym_param_t *param, qpms_normalisation_t norm, qpms_bessel_t J ); /// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging). -qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target, - complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside, +qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(_Complex double *target, + _Complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside, qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm, qpms_bessel_t J ///< Use QPMS_BESSEL_REGULAR to calculate \f$ R^T\f$ or QPMS_HANKEL_PLUS to calculate \f$ Q^T\f$. ); @@ -537,7 +551,7 @@ typedef struct qpms_tmatrix_function_t { /// Convenience function to create a new T-matrix from qpms_tmatrix_function_t. // FIXME the name is not very intuitive. -static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) { +static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, _Complex double omega) { return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega); } @@ -558,9 +572,9 @@ typedef enum { struct qpms_tmatrix_operation_lrmatrix { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - complex double *m; - size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). - bool owns_m; ///< Whether \a m is owned by this; + _Complex double *m; + size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double). + _Bool owns_m; ///< Whether \a m is owned by this; }; struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations. @@ -591,9 +605,9 @@ struct qpms_tmatrix_operation_compose_chain { struct qpms_tmatrix_operation_scmulz { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - complex double *m; - size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). - bool owns_m; ///< Whether \a m is owned by this. + _Complex double *m; + size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double). + _Bool owns_m; ///< Whether \a m is owned by this. }; /// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t. @@ -602,7 +616,7 @@ struct qpms_tmatrix_operation_scmulz { struct qpms_tmatrix_operation_irot3arr { size_t n; ///< Number of rotoreflections; qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. - bool owns_ops; ///< Whether \a ops array is owned by this. + _Bool owns_ops; ///< Whether \a ops array is owned by this. }; /// A generic T-matrix transformation operator. @@ -689,4 +703,7 @@ typedef struct qpms_abstract_particle_t { typedef qpms_particle_tid_t qpms_abstract_particle_tid_t; #endif // 0 +#ifdef __cplusplus +} +#endif #endif //TMATRICES_H diff --git a/qpms/tolerances.h b/qpms/tolerances.h index 14939a6..eac65df 100644 --- a/qpms/tolerances.h +++ b/qpms/tolerances.h @@ -1,6 +1,9 @@ /*! \file tolerances.h */ #ifndef QPMS_TOLERANCES_H #define QPMS_TOLERANCES_H +#ifdef __cplusplus +extern "C" { +#endif // TODO DOC @@ -12,4 +15,7 @@ typedef struct qpms_tolerance_spec_t { /// A rather arbitrary default tolerance. static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-9, .rtol = 1e-8}; +#ifdef __cplusplus +} +#endif #endif // QPMS_TOLERANCES_H diff --git a/qpms/translations.c b/qpms/translations.c index 3ae667f..e502f77 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1,5 +1,6 @@ #include #include "qpms_types.h" +#include #include "qpms_specfunc.h" #include "gaunt.h" #include "translations.h" @@ -15,6 +16,10 @@ #include "normalisation.h" #include "translations_inlines.h" +#if defined LATTICESUMS31 || defined LATTICESUMS32 +#include "lattices.h" +#endif + /* * Define macros with additional factors that "should not be there" according * to the "original" formulae but are needed to work with my vswfs. diff --git a/qpms/translations.h b/qpms/translations.h index 008ff04..492a9c9 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -26,9 +26,12 @@ */ #ifndef QPMS_TRANSLATIONS_H #define QPMS_TRANSLATIONS_H +#ifdef __cplusplus +extern "C" { +#endif + #include "vectors.h" #include "qpms_types.h" -#include #include #include @@ -37,10 +40,10 @@ #endif -complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, +_Complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); -complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, +_Complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); /// Structure holding the constant factors in normalisation operators. @@ -65,8 +68,8 @@ typedef struct qpms_trans_calculator { qpms_normalisation_t normalisation; qpms_l_t lMax; qpms_y_t nelem; - complex double **A_multipliers; - complex double **B_multipliers; + _Complex double **A_multipliers; + _Complex double **B_multipliers; #if 0 // Normalised values of the Legendre functions and derivatives // for θ == π/2, i.e. for the 2D case. @@ -90,44 +93,44 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation /// Destructor for qpms_trans_calculator_t. void qpms_trans_calculator_free(qpms_trans_calculator *); -complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, +_Complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); -complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, +_Complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, + _Complex double *Adest, _Complex double *Bdest, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, + _Complex double *Adest, _Complex double *Bdest, size_t deststride, size_t srcstride, csph_t kdlj, bool r_ge_d, qpms_bessel_t J); // TODO update the types later -complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, complex double kdlj_r, +_Complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, _Complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); -complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, complex double kdlj_r, +_Complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, _Complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, complex double kdlj_r, + _Complex double *Adest, _Complex double *Bdest, + int m, int n, int mu, int nu, _Complex double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, + _Complex double *Adest, _Complex double *Bdest, size_t deststride, size_t srcstride, - complex double kdlj_r, double kdlj_theta, double kdlj_phi, + _Complex double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J); // Convenience functions using VSWF base specs qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c, - complex double *target, + _Complex double *target, /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. @@ -138,12 +141,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator * /// and with automatic \a r_ge_d = `false`. qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( const qpms_trans_calculator *c, - complex double *target, + _Complex double *target, /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - complex double k, cart3_t destpos, cart3_t srcpos, + _Complex double k, cart3_t destpos, cart3_t srcpos, qpms_bessel_t J /// Workspace has to be at least 2 * c->neleme**2 long ); @@ -156,10 +159,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, - complex double *Adest, double *Aerr, - complex double *Bdest, double *Berr, + _Complex double *Adest, double *Aerr, + _Complex double *Bdest, double *Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, - const double eta, const complex double wavenumber, + const double eta, const _Complex double wavenumber, cart2_t b1, cart2_t b2, const cart2_t beta, const cart3_t particle_shift, @@ -167,10 +170,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, ); int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c, - complex double *Adest, double *Aerr, - complex double *Bdest, double *Berr, + _Complex double *Adest, double *Aerr, + _Complex double *Bdest, double *Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, - const double eta, const complex double wavenumber, + const double eta, const _Complex double wavenumber, cart2_t b1, cart2_t b2, const cart2_t beta, const cart3_t particle_shift, @@ -180,12 +183,12 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c, // Convenience functions using VSWF base specs qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c, - complex double *target, double *err, + _Complex double *target, double *err, /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - const double eta, const complex double wavenumber, + const double eta, const _Complex double wavenumber, cart2_t b1, cart2_t b2, const cart2_t beta, const cart3_t particle_shift, @@ -193,12 +196,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat ); qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c, - complex double *target, double *err, + _Complex double *target, double *err, /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. const qpms_vswf_set_spec_t *destspec, size_t deststride, /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. const qpms_vswf_set_spec_t *srcspec, size_t srcstride, - const double eta, const complex double wavenumber, + const double eta, const _Complex double wavenumber, cart2_t b1, cart2_t b2, const cart2_t beta, const cart3_t particle_shift, @@ -212,8 +215,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul // e31z means that the particles are positioned along the z-axis; // their positions and K-values are then denoted by a single z-coordinate int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c, - complex double *Adest, double *Aerr, - complex double *Bdest, double *Berr, + _Complex double *Adest, double *Aerr, + _Complex double *Bdest, double *Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, const double eta, const double k, const double unitcell_area, // just lattice period @@ -244,4 +247,7 @@ int qpms_cython_trans_calculator_get_AB_arrays_loop( #endif //QPMS_COMPILE_PYTHON_EXTENSIONS +#ifdef __cplusplus +} +#endif #endif // QPMS_TRANSLATIONS_H diff --git a/qpms/translations_dbg.c b/qpms/translations_dbg.c index a2f998c..4334008 100644 --- a/qpms/translations_dbg.c +++ b/qpms/translations_dbg.c @@ -5,6 +5,7 @@ #include "tiny_inlines.h" #include "qpms_error.h" #include "normalisation.h" +#include "lattices.h" int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c, diff --git a/qpms/translations_dbg.h b/qpms/translations_dbg.h index 13aa0fc..262c07e 100644 --- a/qpms/translations_dbg.h +++ b/qpms/translations_dbg.h @@ -8,16 +8,16 @@ #include "translations.h" int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c, - complex double *dest, const csph_t kdlj, const qpms_bessel_t J); + _Complex double *dest, const csph_t kdlj, const qpms_bessel_t J); int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c, - complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J); + _Complex double *dest, const _Complex double k, const cart3_t r, const qpms_bessel_t J); #ifdef LATTICESUMS32 int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c, - complex double * const sigmas_total, double * serr_total, - const double eta, const complex double k, + _Complex double * const sigmas_total, double * serr_total, + const double eta, const _Complex double k, const cart2_t b1, const cart2_t b2, const cart2_t beta, const cart3_t particle_shift, diff --git a/qpms/translations_inlines.h b/qpms/translations_inlines.h index 854bdc3..1bfaa19 100644 --- a/qpms/translations_inlines.h +++ b/qpms/translations_inlines.h @@ -6,12 +6,12 @@ /// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix. // TODO DOC static inline void qpms_trans_array_from_AB( - complex double *t, + _Complex double *t, const qpms_vswf_set_spec_t *const t_destspec, const size_t t_deststride, const qpms_vswf_set_spec_t *const t_srcspec, const size_t t_srcstride, - const complex double *const A, const complex double *const B, + const _Complex double *const A, const _Complex double *const B, /// A and B matrices' lMax. /** This also determines their size and stride: they are assumed to * be square matrices of size `nelem * nelem` where diff --git a/qpms/trivialgroup.c b/qpms/trivialgroup.c index f5dd286..ca8dec2 100644 --- a/qpms/trivialgroup.c +++ b/qpms/trivialgroup.c @@ -1,4 +1,5 @@ #include "groups.h" +#include /// Trivial group, with one (reduntant) generator. /** diff --git a/qpms/vectors.h b/qpms/vectors.h index 7267873..30a84d7 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -4,6 +4,7 @@ #ifndef VECTORS_H #define VECTORS_H #include +#include #ifndef M_PI_2 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487) #endif @@ -198,12 +199,12 @@ static inline double cart3_dist(const cart3_t a, const cart3_t b) { return cart3norm(cart3_substract(a,b)); } -static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) { +static inline _Bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) { return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5; } /// Complex 3D vector scaling. -static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) { +static inline ccart3_t ccart3_scale(const _Complex double c, const ccart3_t v) { ccart3_t res = {c * v.x, c * v.y, c * v.z}; return res; } @@ -221,7 +222,7 @@ static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) { } /// Complex 3D cartesian vector "dot product" without conjugation. -static inline complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) { +static inline _Complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) { return a.x * b.x + a.y * b.y + a.z * b.z; } @@ -244,13 +245,13 @@ static inline csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) } /// Complex 3D vector (geographic coordinates) scaling. -static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) { +static inline csphvec_t csphvec_scale(_Complex double c, const csphvec_t v) { csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic}; return res; } /// Complex 3D vector (geographic coordinates) "dot product" without conjugation. -static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) { +static inline _Complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) { //N.B. no complex conjugation done here return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic; } @@ -262,7 +263,7 @@ static inline sph_t sph_scale(double c, const sph_t s) { } /// "Complex spherical" coordinate system scaling. -static inline csph_t sph_cscale(complex double c, const sph_t s) { +static inline csph_t sph_cscale(_Complex double c, const sph_t s) { csph_t res = {c * s.r, s.theta, s.phi}; return res; } @@ -358,8 +359,9 @@ static inline csph_t ccart2csph(const ccart3_t cart) { /// Real 3D cartesian to spherical (complex r) coordinates conversion. See @ref coord_conversions. static inline csph_t cart2csph(const cart3_t cart) { csph_t sph; - sph.r = cart3norm(cart); - sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2; + const double r = cart3norm(cart); + sph.r = r; + sph.theta = sph.r ? acos(cart.z / r) : M_PI_2; sph.phi = atan2(cart.y, cart.x); return sph; } @@ -578,7 +580,6 @@ static inline double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) { " to 1D not allowed."); } - /// Coordinate conversion of point arrays (something to something). /** The dest and src arrays must not overlap */ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest, @@ -589,31 +590,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td sph_t *d = (sph_t *) dest; switch (tsrc & QPMS_COORDS_BITRANGE) { case QPMS_COORDS_SPH: { - const sph_t *s = src; + const sph_t *s = (const sph_t*) src; for(size_t i = 0; i < nmemb; ++i) d[i] = s[i]; return; } break; case QPMS_COORDS_CART3: { - const cart3_t *s = src; + const cart3_t *s = (const cart3_t*) src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart2sph(s[i]); return; } break; case QPMS_COORDS_POL: { - const pol_t *s = src; + const pol_t *s = (const pol_t*) src; for(size_t i = 0; i < nmemb; ++i) d[i] = pol2sph_equator(s[i]); return; } break; case QPMS_COORDS_CART2: { - const cart2_t *s = src; + const cart2_t *s = (const cart2_t*) src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart22sph(s[i]); return; } break; case QPMS_COORDS_CART1: { - const double *s = src; + const double *s = (const double *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart12sph_zaxis(s[i]); return; @@ -627,31 +628,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td cart3_t *d = (cart3_t *) dest; switch (tsrc & QPMS_COORDS_BITRANGE) { case QPMS_COORDS_SPH: { - const sph_t *s = src; + const sph_t *s = (const sph_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = sph2cart(s[i]); return; } break; case QPMS_COORDS_CART3: { - const cart3_t *s = src; + const cart3_t *s = (const cart3_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = s[i]; return; } break; case QPMS_COORDS_POL: { - const pol_t *s = src; + const pol_t *s = (const pol_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = pol2cart3_equator(s[i]); return; } break; case QPMS_COORDS_CART2: { - const cart2_t *s = src; + const cart2_t *s = (const cart2_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart22cart3xy(s[i]); return; } break; case QPMS_COORDS_CART1: { - const double *s = src; + const double *s = (const double *)src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart12cart3z(s[i]); return; @@ -669,13 +670,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed"); break; case QPMS_COORDS_POL: { - const pol_t *s = src; + const pol_t *s = (const pol_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = s[i]; return; } break; case QPMS_COORDS_CART2: { - const cart2_t *s = src; + const cart2_t *s = (const cart2_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = cart2pol(s[i]); return; @@ -696,13 +697,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed"); break; case QPMS_COORDS_POL: { - const pol_t *s = src; + const pol_t *s = (const pol_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = pol2cart(s[i]); return; } break; case QPMS_COORDS_CART2: { - const cart2_t *s = src; + const cart2_t *s = (const cart2_t *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = s[i]; return; @@ -727,7 +728,7 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td QPMS_PR_ERROR("Implicit conversion from 3D to 1D coordinates not allowed"); break; case QPMS_COORDS_CART1: { - const double *s = src; + const double *s = (const double *) src; for(size_t i = 0; i < nmemb; ++i) d[i] = s[i]; return; @@ -899,6 +900,6 @@ static inline cart2_t cart2_from_double_array(const double a[]) { typedef double matrix3d[3][3]; typedef double matrix2d[2][2]; -typedef complex double cmatrix3d[3][3]; -typedef complex double cmatrix2d[2][2]; +typedef _Complex double cmatrix3d[3][3]; +typedef _Complex double cmatrix2d[2][2]; #endif //VECTORS_H diff --git a/qpms/vswf.c b/qpms/vswf.c index 37e728a..2bc6660 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -9,6 +9,7 @@ #include #include "qpms_error.h" #include "normalisation.h" +#include qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() { @@ -52,7 +53,7 @@ qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswf return QPMS_SUCCESS; } -bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a, +_Bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a, const qpms_vswf_set_spec_t *b) { if (a == b) return true; if (a->norm != b->norm) return false; @@ -307,7 +308,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t * qpms_y_t nelem = qpms_lMax2nelem(lMax); csphvec_t *a; QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t)) - csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem; + csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a1 + 2 * nelem; QPMS_ENSURE_SUCCESS(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm)); const csphvec_t *p1 = a1; const csphvec_t *p2 = a2; @@ -478,7 +479,7 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex } qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec, - const cart3_t evalpoint, const void *args, bool add) { + const cart3_t evalpoint, const void *args, _Bool add) { QPMS_UNTESTED; const qpms_incfield_planewave_params_t *p = args; diff --git a/qpms/vswf.h b/qpms/vswf.h index d2d01ac..86e68aa 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -7,6 +7,10 @@ */ #ifndef QPMS_VSWF_H #define QPMS_VSWF_H +#ifdef __cplusplus +extern "C" { +#endif + #include // ssize_t #include "qpms_types.h" #include @@ -16,7 +20,7 @@ /// Calculates the (regular VSWF) expansion coefficients of an external incident field. typedef qpms_errno_t (*qpms_incfield_t)( /// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n. - complex double *target, + _Complex double *target, const qpms_vswf_set_spec_t *bspec, const cart3_t evalpoint, ///< Point at which the VSWF expansion is made. const void *args, ///< Pointer to additional function-specific arguments. @@ -82,7 +86,7 @@ qpms_errno_t qpms_uvswf_fill( /** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist. */ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec, - const complex double *coeffs, ///< SVWF coefficient vector of size setspec->n. + const _Complex double *coeffs, ///< SVWF coefficient vector of size setspec->n. csph_t kr, ///< Evaluation point. qpms_bessel_t btyp); @@ -118,7 +122,7 @@ typedef struct qpms_incfield_planewave_params_t { */ qpms_errno_t qpms_incfield_planewave( /// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n. - complex double *target, + _Complex double *target, const qpms_vswf_set_spec_t *bspec, const cart3_t evalpoint, ///< Point at which the VSWF expansion is made. const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)). @@ -227,23 +231,26 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm); qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir, ccart3_t amplitude, - complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff, + _Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm); qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude, - complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff, + _Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm); csphvec_t qpms_eval_vswf(sph_t where, - complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs, + _Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs, qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm); csphvec_t qpms_eval_vswf_csph(csph_t where, - complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs, + _Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs, qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm); qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm);//NI void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI +#ifdef __cplusplus +} +#endif #endif // QPMS_VSWF_H diff --git a/qpms/wigner.c b/qpms/wigner.c index 07bd43c..3d6acf1 100644 --- a/qpms/wigner.c +++ b/qpms/wigner.c @@ -52,7 +52,7 @@ complex double qpms_wignerD_elem(const qpms_quat_t R, } complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q, - const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) { + const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, _Bool pseudo) { #ifndef NDEBUG qpms_irot3_checkdet(q); #endif diff --git a/tests/catch/CMakeLists.txt b/tests/catch/CMakeLists.txt new file mode 100644 index 0000000..6a82bca --- /dev/null +++ b/tests/catch/CMakeLists.txt @@ -0,0 +1,7 @@ +find_package(Catch2 REQUIRED) +add_executable(catchtest all_includes.C t_vswf.C catch_aux.C) +# TODO ensure that the tests are linked against the "working tree" version of qpms, +# not the installed one +target_link_libraries(catchtest PRIVATE Catch2::Catch2WithMain qpms lapacke) + + diff --git a/tests/catch/all_includes.C b/tests/catch/all_includes.C new file mode 100644 index 0000000..7147355 --- /dev/null +++ b/tests/catch/all_includes.C @@ -0,0 +1,30 @@ +// Just to test that all headers are C++ compatible +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// C++ type conversion issues: +#include +//#include +//#include +//#include diff --git a/tests/catch/t_vswf.C b/tests/catch/t_vswf.C new file mode 100644 index 0000000..118501e --- /dev/null +++ b/tests/catch/t_vswf.C @@ -0,0 +1,44 @@ +#include +#include +#include +#include "catch_aux.h" +using namespace qpmstest; + +constexpr qpms_l_t lMax = 6; + +static std::vector vswfvals_dvector(double r, double theta, double phi, + qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) +{ + std::vector vals(2 /*el and mg for now */ * (lMax * (lMax+2)) + * 3 /* 3D vector components */ * 2 /*for real and imag*/); + qpms_vswf_fill( + nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]), + lMax, sph_t{r, theta, phi}, + btyp, norm); + return vals; +} + +static std::vector vswfvals_dvector_alternative(double r, double theta, double phi, + qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) +{ + std::vector vals(2 /*el and mg for now */ * (lMax * (lMax+2)) + * 3 /* 3D vector components */ * 2 /*for real and imag*/); + qpms_vswf_fill_alternative( + nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]), + lMax, sph_t{r, theta, phi}, + btyp, norm); + return vals; +} + +TEST_CASE("VSWF alternative implementation") { + using Catch::Matchers::Approx; + auto norms = GENERATE(QPMS_NORMALISATION_DEFAULT, QPMS_NORMALISATION_NORM_POWER, + QPMS_NORMALISATION_NORM_SPHARM, QPMS_NORMALISATION_CONVENTION_SCUFF); + auto thetas = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666); + auto phis = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666); + auto rs_positive = GENERATE(0.0001, 0.001, 0.01, 0.1, 1., 10., 100.); + CHECK_THAT( + vswfvals_dvector(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms), + Approx(vswfvals_dvector_alternative(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms)) + ) ; +}