diff --git a/qpms/beyn.h b/qpms/beyn.h index 13a43aa..9395265 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -5,16 +5,15 @@ #define BEYN_H #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 +25,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,11 +46,11 @@ 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_t *beyn_contour_halfellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints, beyn_contour_halfellipse_orientation or); /// 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); @@ -60,10 +59,10 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, dou 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; diff --git a/qpms/ewald.h b/qpms/ewald.h index 8881586..cb0acbe 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -31,7 +31,6 @@ #include #include #include // for inlined lilgamma -#include #include "qpms_types.h" #include "lattices.h" @@ -51,7 +50,7 @@ typedef struct qpms_ewald3_constants_t { /// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`. 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 +64,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 +77,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 +91,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 +99,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 +107,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 +145,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 +171,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 +183,7 @@ int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result); * even if `z1 == z2`, because `-0 == 0` according to IEEE 754. * 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 +200,7 @@ int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result); * Another than principal branch can be selected using non-zero \a m * 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 +217,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 +236,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 +254,26 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_ * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth * 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, @@ -302,11 +301,11 @@ int ewald3_sigma_short( /// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$. 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, diff --git a/qpms/groups.h b/qpms/groups.h index e78e462..771a7e9 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -35,7 +35,7 @@ struct qpms_finite_group_irrep_t { /** The r-th row, c-th column of the representation of the i'th element is retrieved as * m[i * dim * dim + r * dim + c] */ - complex double *m; + _Complex double *m; }; /// A point group with its irreducible representations and some metadata. diff --git a/qpms/kahansum.h b/qpms/kahansum.h index 80ce4a7..350bfa9 100644 --- a/qpms/kahansum.h +++ b/qpms/kahansum.h @@ -4,8 +4,6 @@ #ifndef KAHANSUM_H #define KAHANSUM_H -#include - static inline void kahaninit(double * const sum, double * const compensation) { *sum = 0; *compensation = 0; @@ -19,14 +17,14 @@ static inline void kahanadd(double *sum, double *compensation, double input) { } -static inline void ckahaninit(complex double * const sum, complex double * const compensation) { +static inline void ckahaninit(_Complex double * const sum, _Complex double * const compensation) { *sum = 0; *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; } 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..d01885f 100644 --- a/qpms/materials.h +++ b/qpms/materials.h @@ -4,6 +4,7 @@ #ifndef QPMS_MATERIALS_H #define QPMS_MATERIALS_H #include "qpms_types.h" +#include #include #ifndef SPEED_OF_LIGHT @@ -20,36 +21,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 +68,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 +87,14 @@ extern const qpms_ldparams_t *const QPMS_LDPARAMS_PT; ///< Lorentz-Drude paramet extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten. /// 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 +126,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,11 +149,11 @@ 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)); } diff --git a/qpms/normalisation.h b/qpms/normalisation.h index f84f546..6afa7ac 100644 --- a/qpms/normalisation.h +++ b/qpms/normalisation.h @@ -9,7 +9,6 @@ #include "qpms_types.h" #include "qpms_error.h" #include -#include #include "indexing.h" #include "optim.h" @@ -36,8 +35,8 @@ static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qp * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * 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 +50,8 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * 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 +61,8 @@ static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t no * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * 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 +76,14 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * 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 +94,8 @@ static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * 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 +108,13 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation * Do not use if the C.-S. has already been taken into account e.g. in * 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 +155,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 +193,7 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor * * */ -static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) { +static inline _Complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) { if(m==0) return 0; switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT) & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { 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..2823d47 100644 --- a/qpms/pointgroups.h +++ b/qpms/pointgroups.h @@ -18,9 +18,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; } } diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 3c07a91..6a9b32f 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -11,25 +11,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); /****************************************************************************** diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 7b7bd3b..75aeaca 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -3,10 +3,9 @@ */ #ifndef QPMS_TYPES_H #define QPMS_TYPES_H -#include -#include #include #include +#include #ifndef M_PI_2 #define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L) @@ -191,7 +190,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 +211,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 +260,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 +345,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 +421,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 diff --git a/qpms/quaternions.h b/qpms/quaternions.h index eed1ce8..520d821 100644 --- a/qpms/quaternions.h +++ b/qpms/quaternions.h @@ -93,7 +93,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 +119,7 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) { } /// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$. -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 +202,17 @@ static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) { * The D matrix are calculated using formulae (3), (4), (6), (7) from * 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 +251,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; } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 6f1ee65..7449ede 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -107,7 +107,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 +243,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 +291,7 @@ typedef struct qpms_scatsys_at_omega_t { * */ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym, - 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 +303,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 +354,31 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, * The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves). * 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 +387,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 +400,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 +410,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 +420,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 +435,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 +443,14 @@ void qpms_ss_LU_free(qpms_ss_LU); /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only. 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 +458,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. 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 +466,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. 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 +495,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 +539,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 +561,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 +579,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 +596,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 +618,10 @@ complex double *qpms_orbit_irrep_basis( /// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points. /** 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 +629,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 +650,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 +673,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 +687,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 +715,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 +734,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 +753,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 +773,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 +832,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 +849,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 +868,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 +903,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,9 +913,9 @@ double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, */ ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, 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 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..3e590a3 100644 --- a/qpms/symmetries.h +++ b/qpms/symmetries.h @@ -25,36 +25,36 @@ #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 +74,5 @@ 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); #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.h b/qpms/tmatrices.h index f9e3fd8..b46571e 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -14,7 +14,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 +42,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 +53,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 +64,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 +75,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 +105,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 +149,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 +161,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 +189,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 +200,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 +213,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 +242,9 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( ); /// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$. -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. ); @@ -257,7 +257,7 @@ complex double *qpms_apply_tmatrix( */ 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 +267,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,7 +275,7 @@ 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, /// Source T-matrix, real type is (const qpms_tmatrix_t*). const void *tmatrix_orig ); @@ -295,7 +295,7 @@ extern const gsl_interp_type * gsl_interp_steffen; // struct gsl_interp_accel; // use if lookup proves to be too slow 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 +317,7 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided. 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 +325,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num /** As in qpms_tmatrix_interpolator_eval(), the imaginary part of frequency is discarded! */ 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 +342,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 +357,10 @@ complex double *qpms_mie_coefficients_reflection( /// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory. 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 +372,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 +380,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 +395,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 +404,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 +460,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 +474,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 +501,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 +537,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 +558,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 +591,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 +602,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. diff --git a/qpms/translations.c b/qpms/translations.c index 3ae667f..dc15c97 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" diff --git a/qpms/translations.h b/qpms/translations.h index 008ff04..c5a3155 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -28,7 +28,6 @@ #define QPMS_TRANSLATIONS_H #include "vectors.h" #include "qpms_types.h" -#include #include #include @@ -37,10 +36,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 +64,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 +89,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 +137,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 +155,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, - 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 +166,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c, ); int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c, - 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 +179,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 +192,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat ); qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c, - 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 +211,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul // e31z means that the particles are positioned along the z-axis; // 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 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..b637727 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; } @@ -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..e526ce8 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; @@ -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..736e803 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -16,7 +16,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 +82,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 +118,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,19 +227,19 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm); qpms_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, 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