From 80c8faa26e64e6dda45b9fae1f30899329d71cf3 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 7 Mar 2024 23:19:24 +0100
Subject: [PATCH] Start transition to complex.h
---
src/include/Commons.h | 66 ++++++------
src/include/TransitionMatrix.h | 20 ++--
src/include/algebraic.h | 6 +-
src/include/lapack_calls.h | 4 +-
src/include/sph_subs.h | 78 +++++++-------
src/include/tfrfme.h | 16 +--
src/include/tra_subs.h | 84 +++++++--------
src/libnptm/Commons.cpp | 112 ++++++++++---------
src/libnptm/TransitionMatrix.cpp | 52 +++++----
src/libnptm/algebraic.cpp | 10 +-
src/libnptm/lapack_calls.cpp | 15 +--
src/libnptm/sph_subs.cpp | 178 ++++++++++++++++---------------
src/libnptm/tfrfme.cpp | 52 +++++----
src/libnptm/tra_subs.cpp | 131 +++++++++++------------
14 files changed, 423 insertions(+), 401 deletions(-)
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 6ad2a011..e9545444 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Commons.h
*
* \brief C++ porting of common data structures.
@@ -42,15 +44,15 @@ protected:
int nlmmt;
public:
//! \brief QUESTION: definition?
- std::complex<double> **rmi;
+ dcomplex **rmi;
//! \brief QUESTION: definition?
- std::complex<double> **rei;
+ dcomplex **rei;
//! \brief QUESTION: definition?
- std::complex<double> **w;
+ dcomplex **w;
//! \brief QUESTION: definition?
- std::complex<double> *fsas;
+ dcomplex *fsas;
//! \brief QUESTION: definition?
- std::complex<double> **vints;
+ dcomplex **vints;
//! \brief QUESTION: definition?
double *sscs;
//! \brief QUESTION: definition?
@@ -80,7 +82,7 @@ public:
//! \brief Vector of numbers of spherical layers.
int *nshl;
//! \brief QUESTION: definition?
- std::complex<double> ***sas;
+ dcomplex ***sas;
/*! \brief C1 instance constructor.
*
@@ -105,13 +107,13 @@ class C2 {
int nhspo;
public:
//! \brief QUESTION: definition?
- std::complex<double> *ris;
+ dcomplex *ris;
//! \brief QUESTION: definition?
- std::complex<double> *dlri;
+ dcomplex *dlri;
//! \brief QUESTION: definition?
- std::complex<double> *dc0;
+ dcomplex *dc0;
//! \brief QUESTION: definition?
- std::complex<double> *vkt;
+ dcomplex *vkt;
//! Vector of scaled sizes. QUESTION: correct?
double *vsz;
@@ -133,9 +135,9 @@ public:
class C3 {
public:
//! \brief QUESTION: definition?
- std::complex<double> tfsas;
+ dcomplex tfsas;
//! \brief QUESTION: definition?
- std::complex<double> **tsas;
+ dcomplex **tsas;
//! \brief QUESTION: definition?
double gcs;
//! \brief QUESTION: definition?
@@ -210,35 +212,35 @@ protected:
void allocate_vectors(C4 *c4);
public:
//! \brief QUESTION: definition?
- std::complex<double> *vh;
+ dcomplex *vh;
//! \brief QUESTION: definition?
- std::complex<double> *vj0;
+ dcomplex *vj0;
//! \brief QUESTION: definition?
- std::complex<double> *vj;
+ dcomplex *vj;
//! \brief QUESTION: definition?
- std::complex<double> *vyhj;
+ dcomplex *vyhj;
//! \brief QUESTION: definition?
- std::complex<double> *vyj0;
+ dcomplex *vyj0;
//! \brief QUESTION: definition?
- std::complex<double> **am0m;
+ dcomplex **am0m;
//! \brief QUESTION: definition?
- std::complex<double> *vint;
+ dcomplex *vint;
//! \brief QUESTION: definition?
- std::complex<double> *vintm;
+ dcomplex *vintm;
//! \brief QUESTION: definition?
- std::complex<double> **vints;
+ dcomplex **vints;
//! \brief QUESTION: definition?
- std::complex<double> *vintt;
+ dcomplex *vintt;
//! \brief QUESTION: definition?
- std::complex<double> **fsac;
+ dcomplex **fsac;
//! \brief QUESTION: definition?
- std::complex<double> **sac;
+ dcomplex **sac;
//! \brief QUESTION: definition?
- std::complex<double> **fsacm;
+ dcomplex **fsacm;
//! \brief QUESTION: definition?
double *scsc;
//! \brief QUESTION: definition?
- std::complex<double> *scscp;
+ dcomplex *scscp;
//! \brief QUESTION: definition?
double *ecsc;
//! \brief QUESTION: definition?
@@ -246,11 +248,11 @@ public:
//! \brief QUESTION: definition?
double *scscm;
//! \brief QUESTION: definition?
- std::complex<double> *ecscp;
+ dcomplex *ecscp;
//! \brief QUESTION: definition?
- std::complex<double> *scscpm;
+ dcomplex *scscpm;
//! \brief QUESTION: definition?
- std::complex<double> *ecscpm;
+ dcomplex *ecscpm;
//! \brief QUESTION: definition?
double *v3j0;
//! \brief QUESTION: definition?
@@ -296,11 +298,11 @@ protected:
int sam_size_0;
public:
//! \brief QUESTION: definition?
- std::complex<double> **gis;
+ dcomplex **gis;
//! \brief QUESTION: definition?
- std::complex<double> **gls;
+ dcomplex **gls;
//! \brief QUESTION: definition?
- std::complex<double> **sam;
+ dcomplex **sam;
/*! \brief C9 instance constructor.
*
diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
index 447457f4..3966abda 100644
--- a/src/include/TransitionMatrix.h
+++ b/src/include/TransitionMatrix.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file TransitionMatrix.h
*
* \brief Representation of the Transition Matrix.
@@ -19,7 +21,7 @@ class TransitionMatrix {
//! External medium refractive index.
double exri;
//! Vectorized matrix elements.
- std::complex<double> *elements;
+ dcomplex *elements;
//! Sphere radius.
double sphere_radius;
//! Matrix shape
@@ -66,11 +68,11 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _elements: `complex<double> *` Vectorized elements of the matrix.
+ * \param _elements: `complex double *` Vectorized elements of the matrix.
* \param _radius: `double` Radius for the single sphere case (defaults to 0.0).
*/
TransitionMatrix(
- int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ int _is, int _lm, double _vk, double _exri, dcomplex *_elements,
double _radius=0.0
);
@@ -82,13 +84,13 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _rmi: Matrix of complex.
- * \param _rei: Matrix of complex.
+ * \param _rmi: `complex double **`
+ * \param _rei: `complex double **`
* \param _sphere_radius: `double` Radius of the sphere.
*/
TransitionMatrix(
- int _lm, double _vk, double _exri, std::complex<double> **_rmi,
- std::complex<double> **_rei, double _sphere_radius
+ int _lm, double _vk, double _exri, dcomplex **_rmi,
+ dcomplex **_rei, double _sphere_radius
);
/*! \brief Transition Matrix instance constructor for a cluster of spheres.
@@ -100,9 +102,9 @@ class TransitionMatrix {
* \param _lm: `int` Maximum field expansion order.
* \param _vk: `double` Wave number in scale units.
* \param _exri: `double` External medium refractive index.
- * \param _am0m: Matrix of complex.
+ * \param _am0m: `complex double **`
*/
- TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, std::complex<double> **_am0m);
+ TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, dcomplex **_am0m);
/*! \brief Transition Matrix instance destroyer.
*/
diff --git a/src/include/algebraic.h b/src/include/algebraic.h
index a72d94da..ae256a7d 100644
--- a/src/include/algebraic.h
+++ b/src/include/algebraic.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file algebraic.h
*
* \brief Declaration of algebraic functions with different call-backs.
@@ -23,12 +25,12 @@
/*! \brief Perform in-place matrix inversion.
*
- * \param mat: Matrix of complex. The matrix to be inverted (must be a square matrix).
+ * \param mat: `complex double **` The matrix to be inverted (must be a square matrix).
* \param size: `lapack_int` The size of the matrix (i.e. the number of its rows or columns).
* \param ier: `int &` Reference to an integer variable for returning a result flag.
* \param max_size: `lapack_int` The maximum expected size (required by some call-backs,
* optional, defaults to 0).
*/
-void invert_matrix(std::complex<double> **mat, np_int size, int &ier, np_int max_size=0);
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size=0);
#endif
diff --git a/src/include/lapack_calls.h b/src/include/lapack_calls.h
index 430d9f57..6f2f365b 100644
--- a/src/include/lapack_calls.h
+++ b/src/include/lapack_calls.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file lapack_calss.h
*
* \brief C++ interface to LAPACK calls.
@@ -16,6 +18,6 @@
* \param n: `lapack_int` The number of rows and columns of the [n x n] matrix.
* \param jer: `int &` Reference to an integer return flag.
*/
-void zinvert(std::complex<double> **mat, lapack_int n, int &jer);
+void zinvert(dcomplex **mat, lapack_int n, int &jer);
#endif
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 2aafef22..427ae3a6 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sph_subs.h
*
* \brief C++ porting of SPH functions and subroutines.
@@ -35,11 +37,11 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
* backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
*
* \param n: `int` Order of the function.
- * \param z: `complex<double>` Argument of the function.
+ * \param z: `complex double` Argument of the function.
* \param nm: `int &` Highest computed order.
- * \param csj: Vector of complex. The desired function \f$j\f$.
+ * \param csj: `complex double *` The desired function \f$j\f$.
*/
-void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]);
+void cbf(int n, dcomplex z, int &nm, dcomplex *csj);
/*! \brief Clebsch-Gordan coefficients.
*
@@ -56,10 +58,10 @@ double cg1(int lmpml, int mu, int l, int m);
/*! \brief Conjugate of a double precision complex number
*
- * \param z: `complex<double>` The input complex number.
- * \return result: `complex<double>` The conjugate of the input number.
+ * \param z: `complex double` The input complex number.
+ * \return result: `complex double` The conjugate of the input number.
*/
-std::complex<double> dconjg(std::complex<double> z);
+dcomplex dconjg(dcomplex z);
/*! \brief Comute the continuous variation of the refractive index and of its radial derivative.
*
@@ -95,11 +97,11 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2);
* \param c2: `C2 *` Pointer to a `C2` data structure.
* \param jer: `int &` Reference to integer error code variable.
* \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
- * \param arg: `complex<double> &`
+ * \param arg: `complex double &`
*/
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
+ C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg
);
/*! \brief Bessel function calculation control parameters.
@@ -118,11 +120,11 @@ double envj(int n, double x);
* This function computes the Mueller Transformation Matrix, or Phase Matrix. See
* Sec. 2.8.1 of Borghese, Denti & Saija (2007).
*
- * \param vint: Vector of complex.
+ * \param vint: `complex double *`
* \param cmullr: `double **`
* \param cmul: `double **`
*/
-void mmulc(std::complex<double> *vint, double **cmullr, double **cmul);
+void mmulc(dcomplex *vint, double **cmullr, double **cmul);
/*! \brief Starting point for Bessel function magnitude.
*
@@ -168,14 +170,14 @@ void orunve(double *u1, double *u2, double *u3, int iorth, double torth);
*
* \param up: `double *`
* \param un: `double *`
- * \param ylm: Vector of complex. Field polar spherical harmonics.
+ * \param ylm: `complex double *` Field polar spherical harmonics.
* \param inpol: `int` Incident field polarization type (0 - linear; 1 - circular).
* \param lw: `int`
* \param isq: `int`
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
void pwma(
- double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+ double *up, double *un, dcomplex *ylm, int inpol, int lw,
int isq, C1 *c1
);
@@ -189,14 +191,14 @@ void pwma(
* \param li: `int` Maximum field expansion order.
* \param nsph: `int` Number of spheres.
* \param c1: `C1 *` Pointer to `C1` data structure.
- * \param tqse: Matrix of complex.
- * \param tqspe: Matrix of complex.
- * \param tqss: Matrix of complex.
- * \param tqsps: Matrix of complex.
+ * \param tqse: `double **`
+ * \param tqspe: `complex double **`
+ * \param tqss: `double **`
+ * \param tqsps: `complex double **`
*/
void rabas(
- int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
- double **tqss, std::complex<double> **tqsps
+ int inpol, int li, int nsph, C1 *c1, double **tqse, dcomplex **tqspe,
+ double **tqss, dcomplex **tqsps
);
/*! \brief Real Bessel Function.
@@ -219,18 +221,17 @@ void rbf(int n, double x, int &nm, double sj[]);
*
* \param npntmo: `int`
* \param step: `double`
- * \param dcc: `complex<double>`
+ * \param dcc: `complex double`
* \param x: `double &`
* \param lpo: `int`
- * \param y1: `complex<double> &`
- * \param y2: `complex<double> &`
- * \param dy1: `complex<double> &`
- * \param dy2: `complex<double> &`
+ * \param y1: `complex double &`
+ * \param y2: `complex double &`
+ * \param dy1: `complex double &`
+ * \param dy2: `complex double &`
*/
void rkc(
- int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
- std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
- std::complex<double> &dy2
+ int npntmo, double step, dcomplex dcc, double &x, int lpo,
+ dcomplex &y1, dcomplex &y2, dcomplex &dy1, dcomplex &dy2
);
/*! \brief Transition layer radial function and derivative.
@@ -242,16 +243,15 @@ void rkc(
* \param step: `double`
* \param x: `double &`
* \param lpo: `int`
- * \param y1: `complex<double> &`
- * \param y2: `complex<double> &`
- * \param dy1: `complex<double> &`
- * \param dy2: `complex<double> &`
+ * \param y1: `complex double &`
+ * \param y2: `complex double &`
+ * \param dy1: `complex double &`
+ * \param dy2: `complex double &`
* \param c2: `C2 *` Pointer to a `C2` data structure.
*/
void rkt(
- int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
- std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
- C2 *c2
+ int npntmo, double step, double &x, int lpo, dcomplex &y1,
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
);
/*! \brief Spherical Bessel functions.
@@ -262,9 +262,9 @@ void rkt(
* \param n: `int` Order of the function (from 0 up).
* \param x: `double` Argumento of the function (\f$x > 0\f$).
* \param nm: `int &` Highest computed order.
- * \param sy: `double[]` The desired function \f$y\f$.
+ * \param sy: `double *` The desired function \f$y\f$.
*/
-void rnf(int n, double x, int &nm, double sy[]);
+void rnf(int n, double x, int &nm, double *sy);
/*! \brief Spherical harmonics for given direction.
*
@@ -276,11 +276,11 @@ void rnf(int n, double x, int &nm, double sy[]);
* \param cosrph: `double` Cosine of direction's azimuth.
* \param sinrph: `double` Sine of direction's azimuth.
* \param ll: `int` L value expansion order.
- * \param ylm: Vector of complex. The requested spherical harmonics.
+ * \param ylm: `complex double *` The requested spherical harmonics.
*/
void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
+ int ll, dcomplex *ylm
);
/*! \brief Compute scattering, absorption and extinction cross-sections.
@@ -288,14 +288,14 @@ void sphar(
* This function computes the scattering, absorption and extinction cross-sections in terms
* of Forward Scattering Amplitudes. See Sec. 4.2.1 in Borghese, Denti & Saija (2007).
*
- * \param tfsas: `complex<double> &`
+ * \param tfsas: `complex double &`
* \param nsph: `int` Number of spheres.
* \param lm: `int` Maximum field expansion order.
* \param vk: `double` Wave number in scale units.
* \param exri: `double` External medium refractive index.
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
-void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1);
+void sscr0(dcomplex &tfsas, int nsph, int lm, double vk, double exri, C1 *c1);
/*! \brief C++ Compute the scattering amplitude and the scattered field intensity.
*
diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
index ba2954a7..1dcbb575 100644
--- a/src/include/tfrfme.h
+++ b/src/include/tfrfme.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tfrfme.h
*
* \brief Representation of the trapping calculation objects.
@@ -18,7 +20,7 @@ protected:
int nlmmt;
//! QUESTION: definition?
- std::complex<double> *wk;
+ dcomplex *wk;
/*! \brief Load a Swap1 instance from a HDF5 binary file.
*
@@ -60,9 +62,9 @@ public:
/*! \brief Append an element at the end of the vector.
*
- * \param value: `complex<double>` The value to be added to the vector.
+ * \param value: `complex double` The value to be added to the vector.
*/
- void append(std::complex<double> value) { wk[last_index++] = value; }
+ void append(dcomplex value) { wk[last_index++] = value; }
/*! \brief Load a Swap1 instance from binary file.
*
@@ -83,7 +85,7 @@ public:
/*! \brief Get the pointer to the WK vector.
*
- * \return value: `complex<double> *` Memory address of the WK vector.
+ * \return value: `complex double *` Memory address of the WK vector.
*/
std::complex<double> *get_vector() { return wk; }
@@ -311,7 +313,7 @@ protected:
//! Vector of computed z positions
double *zv;
//! QUESTION: definition?
- std::complex<double> **wsum;
+ dcomplex **wsum;
/*! \brief Load a configuration instance from a HDF5 binary file.
*
@@ -371,9 +373,9 @@ public:
/*! \brief Get the pointer to the WSUM matrix.
*
- * \return value: `complex<double> **` Pointer to the WSUM matrix.
+ * \return value: `complex double **` Pointer to the WSUM matrix.
*/
- std::complex<double> **get_matrix() { return wsum; }
+ dcomplex **get_matrix() { return wsum; }
/*! \brief Calculate the necessary amount of memory to create a new instance.
*
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index c32aa9ae..2260aeb3 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tra_subs.h
*
* \brief C++ porting of TRAPPING functions and subroutines.
@@ -48,17 +50,14 @@ struct CCR {
/*! C++ porting of CAMP
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param am0m: Matrix of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param am0m: `complex double **` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using AM0M and WS.
*/
-void camp(
- std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
- CIL *cil
-);
+void camp(dcomplex *ac, dcomplex **am0m, dcomplex *ws, CIL *cil);
/*! C++ porting of CZAMP
*
@@ -70,26 +69,33 @@ void camp(
*
* This function builds the AC vector using AMD, INDAM and WS.
*/
-void czamp(
- std::complex<double> *ac, std::complex<double> **amd, int **indam,
- std::complex<double> *ws, CIL *cil
-);
+void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex<double> *ws, CIL *cil);
/*! C++ porting of FFRF
*
* \param zpv: `double ****`. QUESTION: definition?
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param fffe: `double *`. QUESTION: definition?
* \param fffs: `double *`. QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
* \param ccr: `CCR *` Pointer to a CCR structure.
*/
void ffrf(
- double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
+ double ****zpv, dcomplex *ac, dcomplex *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
);
+/*! C++ porting of FFRT
+ *
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
+ * \param ffte: `double *`. QUESTION: definition?
+ * \param ffts: `double *`. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ */
+void ffrt(dcomplex *ac, dcomplex *ws, double *ffte, double *ffts, CIL *cil);
+
/*! C++ porting of FRFMER
*
* \param nkv: `int` QUESTION: definition?
@@ -105,68 +111,50 @@ void ffrf(
* \param pmf: `double` QUESTION: definition?
* \param tt1: `Swap1 *` Pointer to first swap object.
* \param tt2: `Swap2 *` Pointer to second swap object.
+ * \return wk: `complex double *` QUESTION: definition?
*/
-std::complex<double> *frfmer(
+dcomplex *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
);
-/*! C++ porting of FFRT
- *
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
- * \param ffte: `double *`. QUESTION: definition?
- * \param ffts: `double *`. QUESTION: definition?
- * \param cil: `CIL *` Pointer to a CIL structure.
- */
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-);
-
/*! C++ porting of PWMALP
*
- * \param w: Matrix of complex. QUESTION: definition?
+ * \param w: `complex double *` QUESTION: definition?
* \param up: `double *`
* \param un: `double *`
- * \param ylm: Vector of complex
+ * \param ylm: `complex double *` Field vector spherical harmonics.
* \param lw: `int`
*/
-void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw);
+void pwmalp(dcomplex **w, double *up, double *un, dcomplex *ylm, int lw);
/*! C++ porting of SAMP
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param tmsm: Vector of complex. QUESTION: definition?
- * \param tmse: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param tmsm: `complex double *` QUESTION: definition?
+ * \param tmse: `complex double *` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using TMSM, TMSE and WS.
*/
-void samp(
- std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
- std::complex<double> *ws, CIL *cil
-);
+void samp(dcomplex *ac, dcomplex *tmsm, dcomplex *tmse, dcomplex *ws, CIL *cil);
/*! C++ porting of SAMPOA
*
- * \param ac: Vector of complex. QUESTION: definition?
- * \param tms: Matrix of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
+ * \param ac: `complex double *` QUESTION: definition?
+ * \param tms: `complex double **` QUESTION: definition?
+ * \param ws: `complex double *` QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
*
* This function builds the AC vector using TMS and WS.
*/
-void sampoa(
- std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
- CIL *cil
-);
+void sampoa(dcomplex *ac, dcomplex **tms, dcomplex *ws, CIL *cil);
/*! C++ porting of WAMFF
*
- * \param wk: `complex<double> *` QUESTION: definition?
+ * \param wk: `complex double *` QUESTION: definition?
* \param x: `double`
* \param y: `double`
* \param z: `double`
@@ -180,6 +168,6 @@ void sampoa(
* \param pmf: `double` QUESTION: definition?
*/
void wamff(
- std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ dcomplex *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
);
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index d19e62b5..658c31c8 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Commons.cpp
*
* DEVELOPMENT NOTE:
@@ -10,38 +12,42 @@
* to the configuration objects. These, on their turn, need to
* expose methods to access the relevant data in read-only mode.
*/
-#include <complex>
+#include <complex.h>
+
+typedef __complex__ double dcomplex;
#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
#endif
-using namespace std;
+double real(dcomplex z) { return __real__ z; }
+
+double imag(dcomplex z) { return __imag__ z; }
C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
nlmmt = 2 * (l_max * (l_max + 2));
nsph = ns;
lm = l_max;
- rmi = new complex<double>*[lm];
- rei = new complex<double>*[lm];
+ rmi = new dcomplex*[lm];
+ rei = new dcomplex*[lm];
for (int ri = 0; ri < lm; ri++) {
- rmi[ri] = new complex<double>[nsph]();
- rei[ri] = new complex<double>[nsph]();
+ rmi[ri] = new dcomplex[nsph]();
+ rei[ri] = new dcomplex[nsph]();
}
- w = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
- vints = new complex<double>*[nsph];
+ vints = new dcomplex*[nsph];
rc = new double*[configurations];
nshl = new int[configurations]();
iog = new int[nsph]();
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
- vints[vi] = new complex<double>[16]();
+ vints[vi] = new dcomplex[16]();
iog[vi] = _iog[vi];
if (iog[vi] >= vi + 1) {
nshl[conf_index] = _nshl[conf_index];
@@ -49,7 +55,7 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
conf_index++;
}
}
- fsas = new complex<double>[nsph]();
+ fsas = new dcomplex[nsph]();
sscs = new double[nsph]();
sexs = new double[nsph]();
sabs = new double[nsph]();
@@ -62,11 +68,11 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
rzz = new double[nsph]();
ros = new double[nsph]();
- sas = new complex<double>**[nsph];
+ sas = new dcomplex**[nsph];
for (int si = 0; si < nsph; si++) {
- sas[si] = new complex<double>*[2];
- sas[si][0] = new complex<double>[2]();
- sas[si][1] = new complex<double>[2]();
+ sas[si] = new dcomplex*[2];
+ sas[si][0] = new dcomplex[2]();
+ sas[si][1] = new dcomplex[2]();
}
}
@@ -112,32 +118,32 @@ C1_AddOns::C1_AddOns(C4 *c4) {
nsph = c4->nsph;
lmpo = c4->lmpo;
nlemt = 2 * c4->nlem;
- vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
- vj0 = new complex<double>[nsph * c4->lmtpo]();
- vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
- vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
- vyj0 = new complex<double>[nsph * c4->lmtpos]();
- am0m = new complex<double>*[nlemt];
+ vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();
+ vj0 = new dcomplex[nsph * c4->lmtpo]();
+ vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
+ vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();
+ vyj0 = new dcomplex[nsph * c4->lmtpos]();
+ am0m = new dcomplex*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
- am0m[ai] = new complex<double>[nlemt]();
+ am0m[ai] = new dcomplex[nlemt]();
}
- vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
- vintm = new complex<double>[16]();
- vintt = new complex<double>[16]();
- vints = new complex<double>*[nsph];
- for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
- fsac = new complex<double>*[2];
- sac = new complex<double>*[2];
- fsacm = new complex<double>*[2];
+ vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
+ vintm = new dcomplex[16]();
+ vintt = new dcomplex[16]();
+ vints = new dcomplex*[nsph];
+ for (int vi = 0; vi < nsph; vi++) vints[vi] = new dcomplex[16]();
+ fsac = new dcomplex*[2];
+ sac = new dcomplex*[2];
+ fsacm = new dcomplex*[2];
for (int fi = 0; fi < 2; fi++) {
- fsac[fi] = new complex<double>[2]();
- sac[fi] = new complex<double>[2]();
- fsacm[fi] = new complex<double>[2]();
+ fsac[fi] = new dcomplex[2]();
+ sac[fi] = new dcomplex[2]();
+ fsacm[fi] = new dcomplex[2]();
}
- scscp = new complex<double>[2]();
- ecscp = new complex<double>[2]();
- scscpm = new complex<double>[2]();
- ecscpm = new complex<double>[2]();
+ scscp = new dcomplex[2]();
+ ecscp = new dcomplex[2]();
+ scscpm = new dcomplex[2]();
+ ecscpm = new dcomplex[2]();
allocate_vectors(c4);
sscs = new double[nsph]();
ecscm = new double[2]();
@@ -191,20 +197,20 @@ void C1_AddOns::allocate_vectors(C4 *c4) {
for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
// This calculates the size of vyhj
int ivy = (nsph * nsph - 1) * c4->litpos;
- vyhj = new complex<double>[ivy]();
+ vyhj = dcomplex[ivy]();
// This calculates the size of vyj0
ivy = c4->lmtpos * c4->nsph;
- vyj0 = new complex<double>[ivy]();
+ vyj0 = new dcomplex[ivy]();
}
C2::C2(int ns, int nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
- ris = new complex<double>[nhspo]();
- dlri = new complex<double>[nhspo]();
- vkt = new complex<double>[nsph]();
- dc0 = new complex<double>[nl]();
+ ris = new dcomplex[nhspo]();
+ dlri = new dcomplex[nhspo]();
+ vkt = new dcomplex[nsph]();
+ dc0 = new dcomplex[nl]();
vsz = new double[nsph]();
}
@@ -217,10 +223,10 @@ C2::~C2() {
}
C3::C3() {
- tsas = new complex<double>*[2];
- tsas[0] = new complex<double>[2];
- tsas[1] = new complex<double>[2];
- tfsas = complex<double>(0.0, 0.0);
+ tsas = new dcomplex*[2];
+ tsas[0] = new dcomplex[2];
+ tsas[1] = new dcomplex[2];
+ tfsas = 0.0 + 0.0 * I;
gcs = 0.0;
scs = 0.0;
ecs = 0.0;
@@ -242,14 +248,14 @@ C6::~C6() {
}
C9::C9(int ndi, int nlem, int ndit, int nlemt) {
- gis = new complex<double>*[ndi];
- gls = new complex<double>*[ndi];
+ gis = new dcomplex*[ndi];
+ gls = new dcomplex*[ndi];
for (int gi = 0; gi < ndi; gi++) {
- gis[gi] = new complex<double>[nlem]();
- gls[gi] = new complex<double>[nlem]();
+ gis[gi] = new dcomplex[nlem]();
+ gls[gi] = new dcomplex[nlem]();
}
- sam = new complex<double>*[ndit];
- for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
+ sam = new dcomplex*[ndit];
+ for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();
gis_size_0 = ndi;
sam_size_0 = ndit;
}
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index 30126ff1..2df1bc4f 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -1,13 +1,17 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file TransitionMatrix.cpp
*
* \brief Implementation of the Transition Matrix structure.
*/
-#include <complex>
+#include <complex.h>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
+typedef __complex__ double dcomplex;
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -26,13 +30,17 @@
using namespace std;
+double real(dcomplex z) { return __real__ z; }
+
+double imag(dcomplex z) { return __imag__ z; }
+
TransitionMatrix::~TransitionMatrix() {
if (elements != NULL) delete[] elements;
if (shape != NULL) delete[] shape;
}
TransitionMatrix::TransitionMatrix(
- int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ int _is, int _lm, double _vk, double _exri, dcomplex *_elements,
double _radius
) {
is = _is;
@@ -53,8 +61,8 @@ TransitionMatrix::TransitionMatrix(
}
TransitionMatrix::TransitionMatrix(
- int _lm, double _vk, double _exri, std::complex<double> **_rmi,
- std::complex<double> **_rei, double _sphere_radius
+ int _lm, double _vk, double _exri, dcomplex **_rmi,
+ dcomplex **_rei, double _sphere_radius
) {
is = 1111;
shape = new int[2];
@@ -64,7 +72,7 @@ TransitionMatrix::TransitionMatrix(
vk = _vk;
exri = _exri;
sphere_radius = _sphere_radius;
- elements = new complex<double>[2 * l_max]();
+ elements = new dcomplex[2 * l_max]();
for (int ei = 0; ei < l_max; ei++) {
elements[2 * ei] = -1.0 / _rmi[ei][0];
elements[2 * ei + 1] = -1.0 / _rei[ei][0];
@@ -73,7 +81,7 @@ TransitionMatrix::TransitionMatrix(
TransitionMatrix::TransitionMatrix(
int _nlemt, int _lm, double _vk, double _exri,
- std::complex<double> **_am0m
+ dcomplex **_am0m
) {
is = 1;
shape = new int[2];
@@ -83,7 +91,7 @@ TransitionMatrix::TransitionMatrix(
vk = _vk;
exri = _exri;
sphere_radius = 0.0;
- elements = new complex<double>[_nlemt * _nlemt]();
+ elements = new dcomplex[_nlemt * _nlemt]();
for (int ei = 0; ei < _nlemt; ei++) {
for (int ej = 0; ej < _nlemt; ej++) elements[_nlemt * ei + ej] = _am0m[ei][ej];
}
@@ -113,7 +121,7 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
double _vk;
double _exri;
// This vector will be passed to the new object. DO NOT DELETE HERE!
- complex<double> *_elements;
+ dcomplex *_elements;
double _radius = 0.0;
status = hdf_file->read("IS", "INT32", &_is);
status = hdf_file->read("L_MAX", "INT32", &_lm);
@@ -125,9 +133,9 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
hid_t file_id = hdf_file->get_file_id();
hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
- _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ _elements[ei] = file_vector[2 * ei] + file_vector[2 * ei + 1] * I;
}
status = H5Dclose(dset_id);
status = hdf_file->read("RADIUS", "FLOAT64", &_radius);
@@ -140,9 +148,9 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
hid_t file_id = hdf_file->get_file_id();
hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
- _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ _elements[ei] = file_vector[2 * ei] + file_vector[2 * ei + 1] * I;
}
status = H5Dclose(dset_id);
tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);
@@ -166,7 +174,7 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
double _vk;
double _exri;
// This vector will be passed to the new object. DO NOT DELETE HERE!
- complex<double> *_elements;
+ dcomplex *_elements;
double _radius = 0.0;
ttms.read(reinterpret_cast<char *>(&_is), sizeof(int));
ttms.read(reinterpret_cast<char *>(&_lm), sizeof(int));
@@ -174,12 +182,12 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
ttms.read(reinterpret_cast<char *>(&_exri), sizeof(double));
if (_is == 1111) {
num_elements = _lm * 2;
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- _elements[ei] = complex<double>(vreal, vimag);
+ _elements[ei] = vreal + vimag * I;
}
double _radius;
ttms.read(reinterpret_cast<char *>(&_radius), sizeof(double));
@@ -187,12 +195,12 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
} else if (_is == 1) {
int nlemt = 2 * _lm * (_lm + 2);
num_elements = nlemt * nlemt;
- _elements = new complex<double>[num_elements]();
+ _elements = new dcomplex[num_elements]();
for (int ei = 0; ei < num_elements; ei++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- _elements[ei] = complex<double>(vreal, vimag);
+ _elements[ei] = vreal + vimag * I;
}
tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements);
}
@@ -239,8 +247,8 @@ void TransitionMatrix::write_hdf5(string file_name) {
int num_elements = 2 * shape[0] * shape[1];
double *ptr_elements = new double[num_elements]();
for (int ei = 0; ei < num_elements / 2; ei++) {
- ptr_elements[2 * ei] = elements[ei].real();
- ptr_elements[2 * ei + 1] = elements[ei].imag();
+ ptr_elements[2 * ei] = real(elements[ei]);
+ ptr_elements[2 * ei + 1] = imag(elements[ei]);
}
rec_ptr_list.append(ptr_elements);
if (is == 1111) {
@@ -287,10 +295,10 @@ void TransitionMatrix::write_legacy(string file_name) {
if (ttms.is_open()) {
int num_elements = shape[0] * shape[1];
for (int ei = 0; ei < num_elements; ei++) {
- complex<double> element1 = elements[ei];
+ dcomplex element1 = elements[ei];
double vreal, vimag;
- vreal = element1.real();
- vimag = element1.imag();
+ vreal = real(element1);
+ vimag = imag(element1);
ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
diff --git a/src/libnptm/algebraic.cpp b/src/libnptm/algebraic.cpp
index 9f81015f..9a4b5f4b 100644
--- a/src/libnptm/algebraic.cpp
+++ b/src/libnptm/algebraic.cpp
@@ -1,9 +1,13 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file algebraic.cpp
*
* \brief Implementation of algebraic functions with different call-backs.
*/
-#include <complex>
+#include <complex.h>
+
+typedef __complex__ double dcomplex;
#ifndef INCLUDE_LAPACK_CALLS_H_
#include "lapacke.h"
@@ -15,12 +19,12 @@
#endif
// >>> FALL-BACK FUNCTIONS DECLARATION <<< //
-extern void lucin(std::complex<double> **mat, np_int max_size, np_int size, int &ier);
+extern void lucin(dcomplex **mat, np_int max_size, np_int size, int &ier);
// >>> END OF FALL-BACK FUNCTIONS <<< //
using namespace std;
-void invert_matrix(std::complex<double> **mat, np_int size, int &ier, np_int max_size) {
+void invert_matrix(dcomplex **mat, np_int size, int &ier, np_int max_size) {
ier = 0;
#ifdef USE_LAPACK
zinvert(mat, size, ier);
diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp
index a13ec001..e8c84de0 100644
--- a/src/libnptm/lapack_calls.cpp
+++ b/src/libnptm/lapack_calls.cpp
@@ -1,20 +1,23 @@
-#include <complex>
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
#include <complex.h>
+typedef __complex__ double dcomplex;
+
#ifndef INCLUDE_LAPACK_CALLS_H_
#include "lapacke.h"
#include "../include/lapack_calls.h"
#endif
#ifdef USE_LAPACK
-void zinvert(std::complex<double> **mat, lapack_int n, int &jer) {
+void zinvert(dcomplex **mat, lapack_int n, int &jer) {
jer = 0;
- __complex__ double *arr = new __complex__ double[n * n];
- const __complex__ double uim = 1.0*I;
+ __complex__ double *arr = new dcomplex[n * n];
+ const dcomplex uim = 0.0 + 1.0 * I;
for (lapack_int i = 0; i < n; i++) {
for (lapack_int j = 0; j < n; j++) {
lapack_int idx = i + n * j;
- arr[idx] = mat[j][i].real() + uim * mat[j][i].imag();
+ arr[idx] = mat[j][i];
}
}
@@ -25,7 +28,7 @@ void zinvert(std::complex<double> **mat, lapack_int n, int &jer) {
for (lapack_int i = 0; i < n; i++) {
for (lapack_int j = 0; j < n; j++) {
lapack_int idx = i + n * j;
- mat[j][i] = std::complex<double>(creal(arr[idx]), cimag(arr[idx]));
+ mat[j][i] = arr[idx];
}
}
delete[] IPIV;
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 30e31bf9..38d3ab77 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -1,8 +1,12 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sph_subs.cpp
*
* \brief C++ implementation of SPHERE subroutines.
*/
-#include <complex>
+#include <complex.h>
+
+typedef __complex__ double dcomplex;
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -12,11 +16,13 @@
#include "../include/sph_subs.h"
#endif
-using namespace std;
+double real(dcomplex z) { return __real__ z; }
+
+double imag(dcomplex z) { return __imag__ z; }
void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
- complex<double> cc0 = complex<double>(0.0, 0.0);
- complex<double> summ, sume, suem, suee, sum;
+ dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex summ, sume, suem, suee, sum;
double half_pi = acos(0.0);
double cofs = half_pi * 2.0 / sqk;
for (int i40 = 0; i40 < nsph; i40++) {
@@ -64,7 +70,7 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
}
}
-void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
+void cbf(int n, dcomplex z, int &nm, dcomplex *csj) {
/*
* FROM CSPHJY OF LIBRARY specfun
*
@@ -86,25 +92,25 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
for (int k = 2; k <= n + 1; k++) {
csj[k - 1] = 0.0;
}
- csj[0] = complex<double>(1.0, 0.0);
+ csj[0] = 1.0 + 0.0 * I;
return;
}
- csj[0] = std::sin(z) / z;
+ csj[0] = csin(z) / z;
if (n == 0) {
return;
}
- csj[1] = (csj[0] -std::cos(z)) / z;
+ csj[1] = (csj[0] -ccos(z)) / z;
if (n == 1) {
return;
}
- complex<double> csa = csj[0];
- complex<double> csb = csj[1];
+ dcomplex csa = csj[0];
+ dcomplex csb = csj[1];
int m = msta1(a0, 200);
if (m < n) nm = m;
else m = msta2(a0, n, 15);
- complex<double> cf0 = 0.0;
- complex<double> cf1 = 1.0e-100;
- complex<double> cf, cs;
+ dcomplex cf0 = 0.0;
+ dcomplex cf1 = 1.0e-100;
+ dcomplex cf, cs;
for (int k = m; k >= 0; k--) {
cf = (2.0 * k + 3.0) * cf1 / z - cf0;
if (k <= nm) csj[k] = cf;
@@ -177,17 +183,17 @@ double cg1(int lmpml, int mu, int l, int m) {
return result;
}
-complex<double> dconjg(std::complex<double> z) {
- double zreal = z.real();
- double zimag = z.imag();
- return complex<double>(zreal, -zimag);
+dcomplex dconjg(dcomplex z) {
+ double zreal = real(z);
+ double zimag = imag(z);
+ return (zreal - zimag * I);
}
void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
const double half_step = 0.5 * dif / npntmo;
double rr = c1->rc[i - 1][ns - 1];
- const complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ const dcomplex delta = c2->dc0[ic] - c2->dc0[ic - 1];
const int kpnt = npntmo + npntmo;
c2->ris[kpnt] = c2->dc0[ic];
c2->dlri[kpnt] = complex<double>(0.0, 0.0);
@@ -204,25 +210,25 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
+ C1 *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg
) {
const int lipo = li + 1;
const int lipt = li + 2;
double *rfj = new double[lipt];
double *rfn = new double[lipt];
- complex<double> cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
- complex<double> rmf[li], drmf[li], ref[li], dref[li];
- complex<double> dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
- complex<double> y1, dy1, y2, dy2, arin, cri, uim;
+ dcomplex cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
+ dcomplex rmf[li], drmf[li], ref[li], dref[li];
+ dcomplex dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
+ dcomplex y1, dy1, y2, dy2, arin, cri;
+ const dcomplex uim = 1.0 * I;
jer = 0;
- uim = complex<double>(0.0, 1.0);
int nstp = npnt - 1;
int nstpts = npntts - 1;
double sz = vk * c1->ros[i - 1];
c2->vsz[i - 1] = sz;
double vkr1 = vk * c1->rc[i - 1][0];
int nsh = c1->nshl[i - 1];
- c2->vkt[i - 1] = std::sqrt(c2->dc0[0]);
+ c2->vkt[i - 1] = csqrt(c2->dc0[0]);
arg = vkr1 * c2->vkt[i - 1];
arin = arg;
bool goto32 = false;
@@ -354,22 +360,22 @@ double envj(int n, double x) {
}
void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
- double sm2 = vint[0].real();
- double s24 = vint[1].real();
- double d24 = vint[1].imag();
- double sm4 = vint[5].real();
- double s23 = vint[8].real();
- double d32 = vint[8].imag();
- double s34 = vint[9].real();
- double d34 = vint[9].imag();
- double sm3 = vint[10].real();
- double s31 = vint[11].real();
- double d31 = vint[11].imag();
- double s21 = vint[12].real();
- double d12 = vint[12].imag();
- double s41 = vint[13].real();
- double d14 = vint[13].imag();
- double sm1 = vint[15].real();
+ double sm2 = real(vint[0]);
+ double s24 = real(vint[1]);
+ double d24 = imag(vint[1]);
+ double sm4 = real(vint[5]);
+ double s23 = real(vint[8]);
+ double d32 = imag(vint[8]);
+ double s34 = real(vint[9]);
+ double d34 = imag(vint[9]);
+ double sm3 = real(vint[10]);
+ double s31 = real(vint[11]);
+ double d31 = imag(vint[11]);
+ double s21 = real(vint[12]);
+ double d12 = imag(vint[12]);
+ double s41 = real(vint[13]);
+ double d14 = imag(vint[13]);
+ double sm1 = real(vint[15]);
cmullr[0][0] = sm2;
cmullr[0][1] = sm3;
cmullr[0][2] = -s23;
@@ -462,7 +468,7 @@ int msta2(double x, int n, int mp) {
return result;
}
-void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
+void orunve(double *u1, double *u2, double *u3, int iorth, double torth) {
if (iorth <= 0) {
double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
double abs_cp = cp;
@@ -481,7 +487,7 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
}
void pwma(
- double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+ double *up, double *un, dcomplex *ylm, int inpol, int lw,
int isq, C1 *c1
) {
const double four_pi = 8.0 * acos(0.0);
@@ -492,19 +498,19 @@ void pwma(
int nlwm = lw * (lw + 2);
int nlwmt = nlwm + nlwm;
const double sqrtwi = 1.0 / sqrt(2.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> cm1 = 0.5 * complex<double>(up[0], up[1]);
- complex<double> cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ const dcomplex uim = 1.0 * I;
+ dcomplex cm1 = 0.5 * (up[0] + up[1] * I);
+ dcomplex cp1 = 0.5 * (up[0] - up[1] * I);
double cz1 = up[2];
- complex<double> cm2 = 0.5 * complex<double>(un[0], un[1]);
- complex<double> cp2 = 0.5 * complex<double>(un[0], -un[1]);
+ dcomplex cm2 = 0.5 * (un[0] + un[1] * I);
+ dcomplex cp2 = 0.5 * (un[0] - un[1] * I);
double cz2 = un[2];
for (int l20 = 0; l20 < lw; l20++) {
int l = l20 + 1;
int lf = l + 1;
int lftl = lf * l;
double x = 1.0 * lftl;
- complex<double> cl = complex<double>(four_pi / sqrt(x), 0.0);
+ dcomplex cl = four_pi / sqrt(x) + 0.0 * I;
for (int powi = 1; powi <= l; powi++) cl *= uim;
int mv = l + lf;
int m = -lf;
@@ -536,8 +542,8 @@ void pwma(
if (inpol != 0) {
for (int k40 = 0; k40 < nlwm; k40++) {
int i = k40 + nlwm;
- complex<double> cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
- complex<double> cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
+ dcomplex cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
+ dcomplex cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
c1->w[k40][ispo - 1] = cc2;
c1->w[i][ispo - 1] = -cc2;
c1->w[k40][ispt - 1] = cc1;
@@ -560,11 +566,11 @@ void pwma(
}
void rabas(
- int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
- double **tqss, std::complex<double> **tqsps
+ int inpol, int li, int nsph, C1 *c1, double **tqse, dcomplex **tqspe,
+ double **tqss, dcomplex **tqsps
) {
- complex<double> cc0 = complex<double>(0.0, 0.0);
- complex<double> uim = complex<double>(0.0, 1.0);
+ dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex uim = 0.0 + 1.0 * I;
double two_pi = 4.0 * acos(0.0);
for (int i80 = 0; i80 < nsph; i80++) {
int i = i80 + 1;
@@ -580,13 +586,13 @@ void rabas(
for (int l70 = 0; l70 < li; l70++) {
int l = l70 + 1;
double fl = 1.0 * (l + l + 1);
- complex<double> rm = 1.0 / c1->rmi[l70][i80];
+ dcomplex rm = 1.0 / c1->rmi[l70][i80];
double rmm = (rm * dconjg(rm)).real();
- complex<double> re = 1.0 / c1->rei[l70][i80];
+ dcomplex re = 1.0 / c1->rei[l70][i80];
double rem = (re * dconjg(re)).real();
if (inpol == 0) {
- complex<double> pce = ((rm + re) * uim) * fl;
- complex<double> pcs = ((rmm + rem) * fl) * uim;
+ dcomplex pce = ((rm + re) * uim) * fl;
+ dcomplex pcs = ((rmm + rem) * fl) * uim;
tqspe[0][i80] -= pce;
tqspe[1][i80] += pce;
tqsps[0][i80] -= pcs;
@@ -672,12 +678,11 @@ void rbf(int n, double x, int &nm, double sj[]) {
}
void rkc(
- int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
- std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
- std::complex<double> &dy2
+ int npntmo, double step, dcomplex dcc, double &x, int lpo,
+ dcomplex &y1, dcomplex &y2, dcomplex &dy1, dcomplex &dy2
) {
- complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
- complex<double> cy4, yy, c14, c21, c22, c23, c24;
+ dcomplex cy1, cdy1, c11, cy23, yc2, c12, c13;
+ dcomplex cy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
@@ -710,12 +715,11 @@ void rkc(
}
void rkt(
- int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
- std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
- C2 *c2
+ int npntmo, double step, double &x, int lpo, dcomplex &y1,
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
) {
- complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
- complex<double> cy4, cdy4, yy, c14, c21, c22, c23, c24;
+ dcomplex cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
+ dcomplex cy4, cdy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
@@ -758,7 +762,7 @@ void rkt(
}
}
-void rnf(int n, double x, int &nm, double sy[]) {
+void rnf(int n, double x, int &nm, double *sy) {
/*
* FROM SPHJY OF LIBRARY specfun
*
@@ -804,7 +808,7 @@ void rnf(int n, double x, int &nm, double sy[]) {
void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
+ int ll, dcomplex *ylm
) {
const int rmp_size = ll;
const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
@@ -867,10 +871,10 @@ void sphar(
lmp = l0p + m;
save = pi4irs * plegn[lmp - 1];
lmy = l0y + m;
- ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
+ ylm[lmy - 1] = save * (cosrmp[m - 1] + sinrmp[m - 1] * I);
if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
lmy = l0y - m;
- ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
+ ylm[lmy - 1] = save * (cosrmp[m - 1] - sinrmp[m - 1 * I]);
label45:
if (m >= l) goto label47;
m += 1;
@@ -881,27 +885,27 @@ void sphar(
goto label40;
}
-void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
- complex<double> sum21, rm, re, csam;
- const complex<double> cc0 = complex<double>(0.0, 0.0);
+void sscr0(dcomplex &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
+ dcomplex sum21, rm, re, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
const double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * (0.0 + 0.5 * I);
tfsas = cc0;
for (int i12 = 0; i12 < nsph; i12++) {
int i = i12 + 1;
int iogi = c1->iog[i12];
if (iogi >= i) {
double sums = 0.0;
- complex<double> sum21 = cc0;
+ dcomplex sum21 = cc0;
for (int l10 = 0; l10 < lm; l10++) {
int l = l10 + 1;
double fl = 1.0 + l + l;
rm = 1.0 / c1->rmi[l10][i12];
re = 1.0 / c1->rei[l10][i12];
- complex<double> rm_cnjg = dconjg(rm);
- complex<double> re_cnjg = dconjg(re);
+ dcomplex rm_cnjg = dconjg(rm);
+ dcomplex re_cnjg = dconjg(re);
sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
sum21 += (rm + re) * fl;
}
@@ -923,10 +927,10 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
}
void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
- std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
- const complex<double> cc0(0.0, 0.0);
+ dcomplex s11, s21, s12, s22, rm, re, csam, cc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double ccs = 1.0 / (vk * vk);
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * (0.0 + 0.5 * I);
const double pigfsq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pigfsq * ccs * ccs);
int nlmm = lm * (lm + 2);
@@ -1111,9 +1115,9 @@ void wmamp(
double *un, C1 *c1
) {
const int ylm_size = (lm + 1) * (lm + 1) + 1;
- complex<double> *ylm = new complex<double>[ylm_size];
+ dcomplex *ylm = new dcomplex[ylm_size];
const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ ylm[nlmp - 1] = 0.0 + 0.0 * I;
if (idot != 0) {
if (idot != 1) {
for (int n40 = 0; n40 < nsph; n40++) {
@@ -1140,9 +1144,9 @@ void wmasp(
int nsph, double *argi, double *args, C1 *c1
) {
const int ylm_size = (lm + 1) * (lm + 1) + 1;
- complex<double> *ylm = new complex<double>[ylm_size];
+ dcomplex *ylm = new dcomplex[ylm_size];
const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ ylm[nlmp - 1] = 0.0 + 0.0 * I;
if (idot != 0) {
if (idot != 1) {
for (int n40 = 0; n40 < nsph; n40++) {
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
index b3478d3c..811264d4 100644
--- a/src/libnptm/tfrfme.cpp
+++ b/src/libnptm/tfrfme.cpp
@@ -1,14 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tfrfme.cpp
*
* \brief Implementation of the trapping calculation objects.
*/
-#include <complex>
+#include <complex.h>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
+typedef __complex__ double dcomplex;
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -27,6 +31,10 @@
using namespace std;
+double real(dcomplex z) { return __real__ z; }
+
+double imag(dcomplex z) { return __imag__ z; }
+
// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
Swap1::Swap1(int lm, int _nkv) {
nkv = _nkv;
@@ -57,8 +65,8 @@ Swap1* Swap1::from_hdf5(string file_name) {
double *elements;
string str_type;
int _nlmmt, _nkv, lm, num_elements, index;
- complex<double> value;
- complex<double> *_wk = NULL;
+ dcomplex value;
+ dcomplex *_wk = NULL;
if (status == 0) {
status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
status = hdf_file->read("NKV", "INT32", &_nkv);
@@ -71,7 +79,7 @@ Swap1* Swap1::from_hdf5(string file_name) {
status = hdf_file->read("WK", str_type, elements);
for (int wi = 0; wi < num_elements / 2; wi++) {
index = 2 * wi;
- value = complex<double>(elements[index], elements[index + 1]);
+ value = elements[index] + elements[index + 1] * I;
_wk[wi] = value;
} // wi loop
delete[] elements;
@@ -84,7 +92,7 @@ Swap1* Swap1::from_hdf5(string file_name) {
Swap1* Swap1::from_legacy(string file_name) {
fstream input;
Swap1 *instance = NULL;
- complex<double> *_wk = NULL;
+ dcomplex *_wk = NULL;
int _nlmmt, _nkv, lm;
double rval, ival;
input.open(file_name.c_str(), ios::in | ios::binary);
@@ -98,7 +106,7 @@ Swap1* Swap1::from_legacy(string file_name) {
for (int j = 0; j < num_elements; j++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- _wk[j] = complex<double>(rval, ival);
+ _wk[j] = rval + ival * I;
}
input.close();
} else {
@@ -109,7 +117,7 @@ Swap1* Swap1::from_legacy(string file_name) {
long Swap1::get_memory_requirement(int lm, int _nkv) {
long size = (long)(3 * sizeof(int));
- size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2) * _nkv * _nkv);
+ size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * _nkv * _nkv);
return size;
}
@@ -142,8 +150,8 @@ void Swap1::write_hdf5(string file_name) {
rec_type_list.append(str_type);
double *ptr_elements = new double[num_elements]();
for (int wi = 0; wi < num_elements / 2; wi++) {
- ptr_elements[2 * wi] = wk[wi].real();
- ptr_elements[2 * wi + 1] = wk[wi].imag();
+ ptr_elements[2 * wi] = real(wk[wi]);
+ ptr_elements[2 * wi + 1] = imag(wk[wi]);
}
rec_ptr_list.append(ptr_elements);
@@ -173,8 +181,8 @@ void Swap1::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
for (int j = 0; j < num_elements; j++) {
- rval = wk[j].real();
- ival = wk[j].imag();
+ rval = real(wk[j]);
+ ival = imag(wk[j]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
}
@@ -591,8 +599,8 @@ TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) {
zv = new double[nzv]();
nlmmt = lm * (lm + 2) * 2;
nrvc = nxv * nyv * nzv;
- wsum = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new complex<double>[nrvc]();
+ wsum = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new dcomplex[nrvc]();
}
TFRFME::~TFRFME() {
@@ -624,8 +632,8 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
double *elements;
string str_type;
int _nlmmt, _nrvc, num_elements;
- complex<double> value;
- complex<double> **_wsum = NULL;
+ dcomplex value;
+ dcomplex **_wsum = NULL;
if (status == 0) {
int lmode, lm, nkv, nxv, nyv, nzv;
double vk, exri, an, ff, tra, spd, frsh, exril;
@@ -668,7 +676,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
int index = 0;
for (int wj = 0; wj < _nrvc; wj++) {
for (int wi = 0; wi < _nlmmt; wi++) {
- value = complex<double>(elements[index], elements[index + 1]);
+ value = elements[index] + elements[index + 1] * I;
_wsum[wi][wj] = value;
index += 2;
} // wi loop
@@ -735,7 +743,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
for (int wi = 0; wi < _nlmmt; wi++) {
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- complex<double> value(rval, ival);
+ dcomplex value = rval + ival * I;
_wsum[wi][wj] = value;
} // wi loop
} // wj loop
@@ -754,7 +762,7 @@ long TFRFME::get_memory_requirement(
int _nrvc = _nxv * _nyv * _nzv;
long size = (long)(8 * sizeof(int) + 8 * sizeof(double));
size += (long)((_nxv + _nyv + _nzv) * sizeof(double));
- size += (long)(_nlmmt * _nrvc * sizeof(complex<double>));
+ size += (long)(_nlmmt * _nrvc * sizeof(dcomplex));
return size;
}
@@ -877,8 +885,8 @@ void TFRFME::write_hdf5(string file_name) {
int index = 0;
for (int wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
- ptr_elements[index++] = wsum[wi][wj].real();
- ptr_elements[index++] = wsum[wi][wj].imag();
+ ptr_elements[index++] = real(wsum[wi][wj]);
+ ptr_elements[index++] = imag(wsum[wi][wj]);
} // wi loop
} // wj loop
rec_ptr_list.append(ptr_elements);
@@ -926,8 +934,8 @@ void TFRFME::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
for (int wj = 0; wj < nrvc; wj++) {
for (int wi = 0; wi < nlmmt; wi++) {
- double rval = wsum[wi][wj].real();
- double ival = wsum[wi][wj].imag();
+ double rval = real(wsum[wi][wj]);
+ double ival = imag(wsum[wi][wj]);
output.write(reinterpret_cast<char *>(&rval), sizeof(double));
output.write(reinterpret_cast<char *>(&ival), sizeof(double));
} // wi loop
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index 47c2b28e..1740872f 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -1,11 +1,15 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file tra_subs.cpp
*
* \brief C++ implementation of TRAPPING subroutines.
*/
#include <cmath>
-#include <complex>
+#include <complex.h>
#include <fstream>
+typedef __complex__ double dcomplex;
+
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif
@@ -22,12 +26,11 @@
#include "../include/tra_subs.h"
#endif
-using namespace std;
+double real(dcomplex z) { return __real__ z; }
-void camp(
- std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
- CIL *cil
-) {
+double imag(dcomplex z) { return __imag__ z; }
+
+void camp(dcomplex *ac, dcomplex **am0m, dcomplex *ws, CIL *cil) {
for (int j = 0; j < cil->nlemt; j++) {
for (int i = 0; i < cil->nlemt; i++) {
ac[j] += (am0m[j][i] * ws[i]);
@@ -35,13 +38,10 @@ void camp(
} // j loop
}
-void czamp(
- std::complex<double> *ac, std::complex<double> **amd, int **indam,
- std::complex<double> *ws, CIL *cil
-) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> summ, sume;
+void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex *ws, CIL *cil) {
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex summ, sume;
for (int im20 = 1; im20 <= cil->mxim; im20++) {
int m = im20 - cil->mxmpo;
int abs_m = (m < 0) ? -m : m;
@@ -65,13 +65,13 @@ void czamp(
}
void ffrf(
- double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
+ double ****zpv, dcomplex *ac, dcomplex *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
complex<double> uimmp, summ, sume, suem, suee;
- complex<double> *gap = new complex<double>[3]();
+ dcomplex *gap = new dcomplex[3]();
for (int imu50 = 1; imu50 <= 3; imu50++) {
int mu = imu50 - 2;
@@ -123,9 +123,9 @@ void ffrf(
sume = -gap[0] * ccr->cimu;
suee = -gap[1] * ccr->cof;
suem = -gap[2] * ccr->cimu;
- fffe[0] = (sume - suem).real();
- fffe[1] = ((sume + suem) * uim).real();
- fffe[2] = suee.real();
+ fffe[0] = real(sume - suem);
+ fffe[1] = real((sume + suem) * uim);
+ fffe[2] = real(suee);
for (int imu90 = 1; imu90 <= 3; imu90++) {
int mu = imu90 - 2;
@@ -177,25 +177,22 @@ void ffrf(
sume = gap[0] * ccr->cimu;
suee = gap[1] * ccr->cof;
suem = gap[2] * ccr->cimu;
- fffs[0] = (sume - suem).real();
- fffs[1] = ((sume + suem) * uim).real();
- fffs[2] = suee.real();
+ fffs[0] = real(sume - suem);
+ fffs[1] = real((sume + suem) * uim);
+ fffs[2] = real(suee);
delete[] gap;
}
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+void ffrt(dcomplex *ac, dcomplex *ws, double *ffte, double *ffts, CIL *cil) {
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
- const complex<double> sq2iti = uim * sq2i;
- complex<double> aca, acw;
- complex<double> *ctqce, *ctqcs;
+ const dcomplex sq2iti = uim * sq2i;
+ dcomplex aca, acw;
+ dcomplex *ctqce, *ctqcs;
- ctqce = new complex<double>[3]();
- ctqcs = new complex<double>[3]();
+ ctqce = new dcomplex[3]();
+ ctqcs = new dcomplex[3]();
for (int l60 = 1; l60 <= cil->le; l60++) {
int lpo = l60 + 1;
int il = l60 * lpo;
@@ -236,27 +233,27 @@ void ffrt(
} // im60 loop
} // l60 loop
ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
- ffte[1] = (sq2iti * (ctqce[0] + ctqce[2])).real();
- ffte[2] = ctqce[1].real();
- ffts[0] = -sq2i * (ctqcs[0] - ctqcs[2]).real();
- ffts[1] = -1.0 * (sq2iti * (ctqcs[0] + ctqcs[2])).real();
- ffts[2] = -1.0 * ctqcs[1].real();
+ ffte[1] = real(sq2iti * (ctqce[0] + ctqce[2]));
+ ffte[2] = real(ctqce[1]);
+ ffts[0] = -sq2i * real(ctqcs[0] - ctqcs[2]);
+ ffts[1] = -1.0 * real(sq2iti * (ctqcs[0] + ctqcs[2]));
+ ffts[2] = -1.0 * real(ctqcs[1]);
delete[] ctqce;
delete[] ctqcs;
}
-complex<double> *frfmer(
+dcomplex *frfmer(
int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
Swap1 *tt1, Swap2 *tt2
) {
const int nlemt = le * (le + 2) * 2;
- const complex<double> cc0(0.0, 0.0);
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double sq = vkm * vkm;
double *_vkv = tt2->get_vector();
double **_vkzm = tt2->get_matrix();
- complex<double> *wk = new complex<double>[nlemt]();
+ dcomplex *wk = new dcomplex[nlemt]();
for (int jy90 = 0; jy90 < nkv; jy90++) {
double vky = _vkv[jy90];
double sqy = vky * vky;
@@ -282,22 +279,22 @@ complex<double> *frfmer(
return wk;
}
-void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
- complex<double> cp1, cm1, cp2, cm2, cl;
- const complex<double> uim(0.0, 1.0);
+void pwmalp(dcomplex **w, double *up, double *un, dcomplex *ylm, int lw) {
+ dcomplex cp1, cm1, cp2, cm2, cl;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double four_pi = acos(0.0) * 8.0;
const int nlwm = lw * (lw + 2);
- cm1 = 0.5 * complex<double>(up[0], up[1]);
- cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ cm1 = 0.5 * (up[0] + up[1] * I);
+ cp1 = 0.5 * (up[0] - up[1] * I);
double cz1 = up[2];
- cm2 = 0.5 * complex<double>(un[0], un[1]);
- cp2 = 0.5 * complex<double>(un[0], -un[1]);
+ cm2 = 0.5 * (un[0] + un[1] * I);
+ cp2 = 0.5 * (un[0] - un[1] * I);
double cz2 =un[2];
for (int l20 = 1; l20 <= lw; l20++) {
int lf = l20 + 1;
int lftl = lf * l20;
double x = 1.0 * lftl;
- complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, 1.0 * l20);
+ dcomplex cl = (four_pi / sqrt(x)) * cpow(uim, 1.0 * l20);
int mv = l20 + lf;
int m = -lf;
for (int mf20 = 1; mf20 <= mv; mf20++) {
@@ -319,10 +316,7 @@ void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<doubl
} // k30 loop
}
-void samp(
- std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
- std::complex<double> *ws, CIL *cil
-) {
+void samp(dcomplex *ac, dcomplex *tmsm, dcomplex *tmse, dcomplex *ws, CIL *cil) {
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
int l = l20 + 1;
@@ -336,13 +330,10 @@ void samp(
} // l20 loop
}
-void sampoa(
- std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
- CIL *cil
-) {
- complex<double> **tm = new complex<double>*[2];
- tm[0] = new complex<double>[2]();
- tm[1] = new complex<double>[2]();
+void sampoa(dcomplex *ac, dcomplex **tms, dcomplex *ws, CIL *cil) {
+ dcomplex **tm = new dcomplex*[2];
+ tm[0] = new dcomplex[2]();
+ tm[1] = new dcomplex[2]();
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
tm[0][0] = tms[l20][0];
@@ -364,23 +355,23 @@ void sampoa(
}
void wamff(
- std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ dcomplex *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);
const int nlmmt = 2 * nlmm;
const int nlmp = nlmm + 2;
- complex<double> **w, *ylm;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ dcomplex **w, *ylm;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
complex<double> cfam, cf1, cf2;
double rho, cph, sph, cth, sth, r;
double ftc1, ftc2;
double *up = new double[3];
double *un = new double[3];
- w = new complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[2]();
- ylm = new complex<double>[nlmp]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[2]();
+ ylm = new dcomplex[nlmp]();
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
@@ -446,7 +437,7 @@ void wamff(
double sthl = sth * rir;
double cthl = sqrt(1.0 - sthl * sthl);
double arg = spd * (z - (r / rir) * cthl);
- cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
+ cfam = (tra * cexp(-apfa * apfa) / sqrt(cthl)) * cexp(uim * arg);
if (lmode == 0) {
if (sth == 0.0) { // label 45
ftc1 = ftcn;
@@ -478,7 +469,7 @@ void wamff(
skip62 = false;
}
} else { // label 57
- double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
+ double fam = tra * cexp(-apfa * apfa) / sqrt(cth);
if (lmode == 0) {
double f1 = cph * fam;
double f2 = -sph * fam;
--
GitLab