diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 06e27ab69fb08a81506685406c779f61ff6e9fd7..bcdb133fe2c3fda5c51ca20e71eff4828a7f401c 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -1,13 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file cluster.cpp
*
* \brief Implementation of the calculation for a cluster of spheres.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -47,10 +52,6 @@
#include "../include/algebraic.h"
#endif
-//#ifdef LAPACK_ILP64
-//#define USE_LAPACK
-//#endif
-
using namespace std;
/*! \brief C++ implementation of CLU
@@ -136,17 +137,16 @@ void cluster(string config_file, string data_file, string output_path) {
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0, lcalc = 0;
- complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);
+ dcomplex arg = 0.0 + 0.0 * I;
+ dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (sconf->iog_vec[ci -1] >= ci) configurations++;
}
C2 *c2 = new C2(nsph, configurations, npnt, npntts);
- complex<double> **am = new complex<double>*[mxndm];
- //complex<double> **tam = new complex<double>*[mxndm];
+ dcomplex **am = new dcomplex*[mxndm];
for (int ai = 0; ai < mxndm; ai++) {
- am[ai] = new complex<double>[mxndm]();
- //tam[ai] = new complex<double>[mxndm]();
+ am[ai] = new dcomplex[mxndm]();
}
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
@@ -154,35 +154,35 @@ void cluster(string config_file, string data_file, string output_path) {
double *tqev = new double[3]();
double *tqsv = new double[3]();
double **tqse, **tqss, **tqce, **tqcs;
- complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
+ dcomplex **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2];
- tqspe = new complex<double>*[2];
+ tqspe = new dcomplex*[2];
tqss = new double*[2];
- tqsps = new complex<double>*[2];
+ tqsps = new dcomplex*[2];
tqce = new double*[2];
- tqcpe = new complex<double>*[2];
+ tqcpe = new dcomplex*[2];
tqcs = new double*[2];
- tqcps = new complex<double>*[2];
+ tqcps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph]();
- tqspe[ti] = new complex<double>[nsph]();
+ tqspe[ti] = new dcomplex[nsph]();
tqss[ti] = new double[nsph]();
- tqsps[ti] = new complex<double>[nsph]();
+ tqsps[ti] = new dcomplex[nsph]();
tqce[ti] = new double[3]();
- tqcpe[ti] = new complex<double>[3]();
+ tqcpe[ti] = new dcomplex[3]();
tqcs[ti] = new double[3]();
- tqcps[ti] = new complex<double>[3]();
+ tqcps[ti] = new dcomplex[3]();
}
double *gapv = new double[3]();
- complex<double> **gapp, **gappm;
+ dcomplex **gapp, **gappm;
double **gap, **gapm;
- gapp = new complex<double>*[3];
- gappm = new complex<double>*[3];
+ gapp = new dcomplex*[3];
+ gappm = new dcomplex*[3];
gap = new double*[3];
gapm = new double*[3];
for (int gi = 0; gi < 3; gi++) {
- gapp[gi] = new complex<double>[2]();
- gappm[gi] = new complex<double>[2]();
+ gapp[gi] = new dcomplex[2]();
+ gappm[gi] = new dcomplex[2]();
gap[gi] = new double[2]();
gapm[gi] = new double[2]();
}
@@ -215,7 +215,7 @@ void cluster(string config_file, string data_file, string output_path) {
double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
- fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
+ fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
@@ -361,7 +361,6 @@ void cluster(string config_file, string data_file, string output_path) {
cms(am, c1, c1ao, c4, c6);
lapack_int ndit = 2 * nsph * c4->nlim;
invert_matrix(am, ndit, jer, mxndm);
- //zinvert(am, ndit, jer);
if (jer != 0) break; // jxi488 loop: goes to memory clean
ztm(am, c1, c1ao, c4, c6, c9);
if (sconf->idfc >= 0) {
@@ -383,9 +382,8 @@ void cluster(string config_file, string data_file, string output_path) {
// label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0);
double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
- std::complex<double> s0(0.0, 0.0);
+ dcomplex s0 = 0.0 + 0.0 * I;
scr0(vk, exri, c1, c1ao, c3, c4);
- //printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, c4->li, nsph, c1, sqk, gaps);
rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
@@ -401,19 +399,19 @@ void cluster(string config_file, string data_file, string output_path) {
if (c1->nshl[i] != 1) {
fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
} else { // label 162
- fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
+ fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], real(c2->vkt[i]), imag(c2->vkt[i]));
}
// label 164
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i]), imag(c1->fsas[i]));
csch = 2.0 * vk * sqsfi / c1->gcsv[i];
s0 = c1->fsas[i] * exri;
- qschu = s0.imag() * csch;
- pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
+ qschu = imag(s0) * csch;
+ pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
double rapr = c1->sexs[i] - gaps[i];
double cosav = gaps[i] / c1->sscs[i];
@@ -422,12 +420,12 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
}
} // i170 loop
- fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
+ fprintf(output, " FSAT=%15.7lE%15.7lE\n", real(c3->tfsas), imag(c3->tfsas));
csch = 2.0 * vk * sqsfi / c3->gcs;
s0 = c3->tfsas * exri;
- qschu = s0.imag() * csch;
- pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
+ qschu = imag(s0) * csch;
+ pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
@@ -537,24 +535,24 @@ void cluster(string config_file, string data_file, string output_path) {
for (int i = 0; i < 2; i++) {
double value = c1ao->scscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscpm[i].real();
+ value = real(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscpm[i].imag();
+ value = imag(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscpm[i].real();
+ value = real(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscpm[i].imag();
+ value = imag(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gapm[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gappm[i][j].real();
+ value = real(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gappm[i][j].imag();
+ value = imag(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -575,12 +573,12 @@ void cluster(string config_file, string data_file, string output_path) {
double absrm = abssm / c3->acs;
double acsecs = c3->acs / c3->ecs;
if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
- complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
- double qschum = s0m.imag() * csch;
- double pschum = s0m.real() * csch;
- double s0magm = abs(s0m) * cs0;
- double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real();
- double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag();
+ dcomplex s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
+ double qschum = imag(s0m) * csch;
+ double pschum = real(s0m) * csch;
+ double s0magm = cabs(s0m) * cs0;
+ double rfinrm = real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(c3->tfsas);
+ double extcrm = imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 206
@@ -595,9 +593,9 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr210, ilr210, c1ao->fsacm[ilr210 - 1][ilr210 - 1].real(),
- c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,
- c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()
+ ilr210, ilr210, real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
+ imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
+ real(c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(c1ao->fsacm[jlr - 1][ilr210 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
@@ -610,8 +608,8 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " Fk=%15.7lE\n", fz);
} // ilr210 loop
- double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();
- double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();
+ double rmbrif = (real(c1ao->fsacm[0][0]) - real(c1ao->fsacm[1][1])) / real(c1ao->fsacm[0][0]);
+ double rmdchr = (imag(c1ao->fsacm[0][0]) - imag(c1ao->fsacm[1][1])) / imag(c1ao->fsacm[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
}
@@ -641,13 +639,13 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " SPHERE %2d\n", i226);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
- c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
+ real(c1->sas[i226 - 1][0][0]), imag(c1->sas[i226 - 1][0][0]),
+ real(c1->sas[i226 - 1][1][0]), imag(c1->sas[i226 - 1][1][0])
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[i226 - 1][0][1].real(), c1->sas[i226 - 1][0][1].imag(),
- c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()
+ real(c1->sas[i226 - 1][0][1]), imag(c1->sas[i226 - 1][0][1]),
+ real(c1->sas[i226 - 1][1][1]), imag(c1->sas[i226 - 1][1][1])
);
for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
@@ -671,13 +669,13 @@ void cluster(string config_file, string data_file, string output_path) {
} // i226 loop
fprintf(
output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
- c3->tsas[0][0].real(), c3->tsas[0][0].imag(),
- c3->tsas[1][0].real(), c3->tsas[1][0].imag()
+ real(c3->tsas[0][0]), imag(c3->tsas[0][0]),
+ real(c3->tsas[1][0]), imag(c3->tsas[1][0])
);
fprintf(
output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
- c3->tsas[0][1].real(), c3->tsas[0][1].imag(),
- c3->tsas[1][1].real(), c3->tsas[1][1].imag()
+ real(c3->tsas[0][1]), imag(c3->tsas[0][1]),
+ real(c3->tsas[1][1]), imag(c3->tsas[1][1])
);
fprintf(output, " CLUSTER\n");
pcros(vk, exri, c1, c1ao, c4);
@@ -695,24 +693,24 @@ void cluster(string config_file, string data_file, string output_path) {
for (int i = 0; i < 2; i++) {
double value = c1ao->scsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscp[i].real();
+ value = real(c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->scscp[i].imag();
+ value = imag(c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscp[i].real();
+ value = real(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscp[i].imag();
+ value = imag(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gap[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gapp[i][j].real();
+ value = real(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = gapp[i][j].imag();
+ value = imag(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -720,9 +718,9 @@ void cluster(string config_file, string data_file, string output_path) {
for (int j = 0; j < 3; j++) {
double value = tqce[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcpe[i][j].real();
+ value = real(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcpe[i][j].imag();
+ value = imag(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -730,9 +728,9 @@ void cluster(string config_file, string data_file, string output_path) {
for (int j = 0; j < 3; j++) {
double value = tqcs[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcps[i][j].real();
+ value = real(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = tqcps[i][j].imag();
+ value = imag(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -747,9 +745,9 @@ void cluster(string config_file, string data_file, string output_path) {
}
// label 254
for (int i = 0; i < 16; i++) {
- double value = c1ao->vint[i].real();
+ double value = real(c1ao->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->vint[i].imag();
+ value = imag(c1ao->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
@@ -775,11 +773,11 @@ void cluster(string config_file, string data_file, string output_path) {
double ratio = c3->acs / c3->ecs;
if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
- double qschu = s0.imag() * csch;
- double pschu = s0.real() * csch;
- s0mag = abs(s0) * cs0;
- double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
- double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
+ double qschu = imag(s0) * csch;
+ double pschu = real(s0) * csch;
+ s0mag = cabs(s0) * cs0;
+ double refinr = real(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(c3->tfsas);
+ double extcor = imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 273
@@ -803,13 +801,13 @@ void cluster(string config_file, string data_file, string output_path) {
);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, c1ao->fsac[ilr290 - 1][ilr290 - 1].real(), c1ao->fsac[ilr290 - 1][ilr290 - 1].imag(),
- jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()
+ ilr290, ilr290, real(c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(c1ao->fsac[jlr - 1][ilr290 - 1]), imag(c1ao->fsac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, c1ao->sac[ilr290 - 1][ilr290 - 1].real(), c1ao->sac[ilr290 - 1][ilr290 - 1].imag(),
- jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()
+ ilr290, ilr290, real(c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(c1ao->sac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(c1ao->sac[jlr - 1][ilr290 - 1]), imag(c1ao->sac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
@@ -851,8 +849,8 @@ void cluster(string config_file, string data_file, string output_path) {
);
}
} //ilr290 loop
- double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();
- double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();
+ double rbirif = (real(c1ao->fsac[0][0]) - real(c1ao->fsac[1][1])) / real(c1ao->fsac[0][0]);
+ double rdichr = (imag(c1ao->fsac[0][0]) - imag(c1ao->fsac[1][1])) / imag(c1ao->fsac[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
fprintf(output, " MULC\n");
@@ -874,9 +872,9 @@ void cluster(string config_file, string data_file, string output_path) {
// Some implicit loops writing to binary.
for (int i = 0; i < 16; i++) {
double value;
- value = c1ao->vintm[i].real();
+ value = real(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->vintm[i].imag();
+ value = imag(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
diff --git a/src/cluster/np_cluster.cpp b/src/cluster/np_cluster.cpp
index 73caa518f9d330b8ad4114a30f123ba320035756..66fb783bc2e9d51e32b5dbc6fb0c033793348936 100644
--- a/src/cluster/np_cluster.cpp
+++ b/src/cluster/np_cluster.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file np_cluster.cpp
*
* \brief Cluster of spheres scattering problem handler.
@@ -17,9 +19,12 @@
*/
#include <cstdio>
-#include <complex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index f3c4df91bcc64e6673865b739b5f7dd9a06af5a3..e8e3b8951b4b01538cbf4f541b13b100d84e3e94 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Configuration.h
*
* \brief Configuration data structures.
@@ -168,7 +170,7 @@ class ScattererConfiguration {
friend void sphere(std::string, std::string, std::string);
protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
- std::complex<double> ***dc0_matrix;
+ dcomplex ***dc0_matrix;
//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
double *radii_of_spheres;
//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
@@ -258,7 +260,7 @@ public:
* \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,
* \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor
* for dimensions).
- * \param dc_matrix: Matrix of reference dielectric constants.
+ * \param dc_matrix: `complex double ***` Matrix of reference dielectric constants.
* \param has_external: `bool` Flag to set whether to add an external spherical layer.
* \param exdc: `double` External medium dielectric constant.
* \param wp: `double` wp
@@ -274,7 +276,7 @@ public:
int *nshl_vector,
double **rcf_vector,
int dielectric_func_type,
- std::complex<double> ***dc_matrix,
+ dcomplex ***dc_matrix,
bool has_external,
double exdc,
double wp,
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 6e0587d982b999b1fc907cb7434c21ac5d14718f..6852a46906998e878b0021dc5d606ac8500e8322 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clu_subs.h
*
* \brief C++ porting of CLU functions and subroutines.
@@ -31,15 +33,15 @@
*
* \param zpv: `double ****`
* \param le: `int`
- * \param am0m: Matrix of complex.
- * \param w: Matrix of complex.
+ * \param am0m: `complex double **`
+ * \param w: `complex double **`
* \param sqk: `double`
* \param gapr: `double **`
- * \param gapp: Matrix of complex.
+ * \param gapp: `complex double **`
*/
void apc(
- double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
- double sqk, double **gapr, std::complex<double> **gapp
+ double ****zpv, int le, dcomplex **am0m, dcomplex **w,
+ double sqk, double **gapr, dcomplex **gapp
);
/*! \brief Compute the asymmetry-corrected scattering cross-section under random average
@@ -50,32 +52,29 @@ void apc(
*
* \param zpv: `double ****`
* \param le: `int`
- * \param am0m: Matrix of complex.
+ * \param am0m: `complex double **`
* \param inpol: `int` Polarization type.
* \param sqk: `double`
* \param gaprm: `double **`
- * \param gappm: Matrix of complex.
+ * \param gappm: `complex double **`
*/
void apcra(
- double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
- double **gaprm, std::complex<double> **gappm
+ double ****zpv, const int le, dcomplex **am0m, int inpol, double sqk,
+ double **gaprm, dcomplex **gappm
);
/*! \brief Complex inner product.
*
* This function performs the complex inner product. It is used by `lucin()`.
*
- * \param z: `complex<double>`
- * \param am: Matrix of complex.
+ * \param z: `complex double`
+ * \param am: `complex double **`
* \param i: `int`
* \param jf: `int`
* \param k: `int`
* \param nj: `int`
*/
-std::complex<double> cdtp(
- std::complex<double> z, std::complex<double> **am, int i, int jf,
- int k, int nj
-);
+dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj);
/*! \brief C++ porting of CGEV. QUESTION: description?
*
@@ -92,13 +91,13 @@ double cgev(int ipamo, int mu, int l, int m);
* This function constructs the multi-centered M-matrix of the cluster, according
* to Eq. (5.28) of Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
+ * \param am: `complex double **`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
-void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
+void cms(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
/*! \brief Compute orientation-averaged scattered field intensity.
*
@@ -132,9 +131,9 @@ void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
-std::complex<double> ghit(
- int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
- C1_AddOns *c1ao, C4 *c4, C6 *c6
+dcomplex ghit(
+ int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
+ C1_AddOns *c1ao, C4 *c4, C6 *c6
);
/*! \brief Compute Hankel funtion and Bessel functions.
@@ -152,7 +151,7 @@ std::complex<double> ghit(
* \param c4: `C4 *`
*/
void hjv(
- double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
+ double exri, double vk, int &jer, int &lcalc, dcomplex &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
);
@@ -161,12 +160,12 @@ void hjv(
* This function performs the inversion of the multi-centered M-matrix through
* LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
+ * \param am: `complex double **`
* \param nddmst: `const int64_t`
* \param n: `int64_t`
* \param ier: `int &`
*/
-void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier);
+void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier);
/*! \brief Compute the average extinction cross-section.
*
@@ -180,7 +179,7 @@ void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier);
* \param cextlr: `double **`
* \param cext: `double **`
*/
-void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext);
+void mextc(double vk, double exri, dcomplex **fsac, double **cextlr, double **cext);
/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.
*
@@ -274,16 +273,16 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
* in Borghese, Denti & Saija (2007).
*
* \param le: `int`
- * \param am0m: Matrix of complex.
- * \param w: Matrix of complex.
+ * \param am0m: `complex double **`
+ * \param w: `complex double **`
* \param tqce: `double **`
- * \param tqcpe: Matrix of complex.
+ * \param tqcpe: `complex double **`
* \param tqcs: `double **`
- * \param tqcps: Matrix of complex.
+ * \param tqcps: `complex double **`
*/
void raba(
- int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
- std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
+ int le, dcomplex **am0m, dcomplex **w, double **tqce,
+ dcomplex **tqcpe, double **tqcs, dcomplex **tqcps
);
/*! \brief Compute the radiation force Cartesian components.
@@ -391,13 +390,13 @@ void tqr(
* This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
* of Borghese, Denti & Saija (2007).
*
- * \param am: Matrix of complex.
+ * \param am: `complex double **`
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
* \param c4: `C4 *` Pointer to a C4 structure.
* \param c6: `C6 *` Pointer to a C6 instance.
* \param c9: `C9 *` Pointer to a C9 instance.
*/
-void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9);
+void ztm(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9);
#endif
diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
index 1dcbb5759d0f99b9142dccdd9ea984f57c52085b..2e55c4822b1d186bb7a8094e6ba89aafa46f176c 100644
--- a/src/include/tfrfme.h
+++ b/src/include/tfrfme.h
@@ -87,7 +87,7 @@ public:
*
* \return value: `complex double *` Memory address of the WK vector.
*/
- std::complex<double> *get_vector() { return wk; }
+ dcomplex *get_vector() { return wk; }
/*! \brief Bring the pointer to the next element at the start of vector.
*/
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 2260aeb351dea239ac9570f25ccd11a1950bb294..37e728507fe71f9444d7f1fd405044b2aad6e3e4 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -69,7 +69,7 @@ void camp(dcomplex *ac, dcomplex **am0m, dcomplex *ws, CIL *cil);
*
* This function builds the AC vector using AMD, INDAM and WS.
*/
-void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex<double> *ws, CIL *cil);
+void czamp(dcomplex *ac, dcomplex **amd, int **indam, dcomplex *ws, CIL *cil);
/*! C++ porting of FFRF
*
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 658c31c80b2003a7d53fc62d471ed02e24dde7ef..fc9498a839ede48183ce665b3fd15c01d8e404cb 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -12,18 +12,14 @@
* to the configuration objects. These, on their turn, need to
* expose methods to access the relevant data in read-only mode.
*/
-#include <complex.h>
-
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
#endif
-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;
@@ -197,7 +193,7 @@ 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 = dcomplex[ivy]();
+ vyhj = new dcomplex[ivy]();
// This calculates the size of vyj0
ivy = c4->lmtpos * c4->nsph;
vyj0 = new dcomplex[ivy]();
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 04cdb5c03a0930affa172afb46d7db960166914b..1bf78b7db3e571858bee6ad893f0c575dbb33644 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -1,10 +1,11 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file Configuration.cpp
*
* \brief Implementation of the configuration classes.
*/
#include <cmath>
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
@@ -12,6 +13,10 @@
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -202,7 +207,7 @@ ScattererConfiguration::ScattererConfiguration(
int *nshl_vector,
double **rcf_vector,
int dielectric_func_type,
- complex<double> ***dc_matrix,
+ dcomplex ***dc_matrix,
bool is_external,
double ex,
double w,
@@ -477,12 +482,12 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
} // ns112 loop
} // i113 loop
- complex<double> ***dc0m = new complex<double>**[configurations];
- complex<double> *dc0 = new complex<double>[configurations]();
+ dcomplex ***dc0m = new dcomplex**[configurations];
+ dcomplex *dc0 = new dcomplex[configurations]();
int dim3 = (_idfc == 0) ? nxi : 1;
for (int di = 0; di < configurations; di++) {
- dc0m[di] = new complex<double>*[nsph];
- for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[dim3]();
+ dc0m[di] = new dcomplex*[nsph];
+ for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new dcomplex[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop
@@ -504,7 +509,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
double ival = stod(str_number);
- dc0[ic157] = complex<double>(rval, ival);
+ dc0[ic157] = rval + ival * I;
dc0m[ic157][i162 - 1][jxi468 - 1] = dc0[ic157];
} // ic157 loop
}
@@ -547,7 +552,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
double *xi_vec;
double *ros_vector;
double **rcf_vector;
- complex<double> ***dc0m;
+ dcomplex ***dc0m;
status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
status = hdf_file->read("IES", "INT32_(1)", &ies);
status = hdf_file->read("EXDC", "FLOAT64_(1)", &_exdc);
@@ -591,17 +596,17 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
str_name = "DC0M";
str_type = "FLOAT64_(" + to_string(element_size) + ")";
status = hdf_file->read(str_name, str_type, elements);
- dc0m = new complex<double>**[configuration_count];
+ dc0m = new dcomplex**[configuration_count];
int dc_index = 0;
for (int di = 0; di < configuration_count; di++) {
- dc0m[di] = new complex<double>*[nsph];
+ dc0m[di] = new dcomplex*[nsph];
for (int dj = 0; dj < nsph; dj++) {
- dc0m[di][dj] = new complex<double>[dim3]();
+ dc0m[di][dj] = new dcomplex[dim3]();
for (int dk = 0; dk < dim3; dk++) {
double rval = elements[2 * dc_index];
double ival = elements[2 * dc_index + 1];
dc_index++;
- dc0m[di][dj][dk] = complex<double>(rval, ival);
+ dc0m[di][dj][dk] = rval + ival * I;
}
} // dj loop
} // di loop
@@ -638,7 +643,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
double *_xi_vec;
double *_ros_vec;
double **_rcf_vec;
- complex<double> ***_dc0m;
+ dcomplex ***_dc0m;
fstream input;
input.open(file_name.c_str(), ios::in | ios::binary);
@@ -672,11 +677,11 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
for (int nsi = 0; nsi < nsh; nsi++)
input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));
}
- _dc0m = new complex<double>**[configurations];
+ _dc0m = new dcomplex**[configurations];
int dim3 = (_idfc == 0) ? _nxi : 1;
for (int di = 0; di < configurations; di++) {
- _dc0m[di] = new complex<double>*[_nsph];
- for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[dim3]();
+ _dc0m[di] = new dcomplex*[_nsph];
+ for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3]();
} // di loop
for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
if (_idfc != 0 && jxi468 > 1) continue;
@@ -689,7 +694,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
double rval, ival;
input.read(reinterpret_cast<char *>(&rval), sizeof(double));
input.read(reinterpret_cast<char *>(&ival), sizeof(double));
- _dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(rval, ival);
+ _dc0m[i157][i162 - 1][jxi468 - 1] = rval + ival * I;
}
}
}
@@ -771,7 +776,7 @@ void ScattererConfiguration::print() {
printf(" [");
for (int k = 0; k < number_of_scales; k++) {
if (idfc != 0 and k > 0) continue;
- printf("\t%lg + i(%lg)", dc0_matrix[i][j][k].real(), dc0_matrix[i][j][k].imag());
+ printf("\t%lg + i(%lg)", real(dc0_matrix[i][j][k]), imag(dc0_matrix[i][j][k]));
}
printf("\t]\n");
}
@@ -862,8 +867,8 @@ void ScattererConfiguration::write_hdf5(string file_name) {
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
double dc0_real, dc0_imag;
- dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
- dc0_imag = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ dc0_real = real(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
+ dc0_imag = imag(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
dc0m[2 * dc0_index] = dc0_real;
dc0m[2 * dc0_index + 1] = dc0_imag;
dc0_index++;
@@ -926,8 +931,8 @@ void ScattererConfiguration::write_legacy(string file_name) {
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
double dc0_real, dc0_img;
- dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
- dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+ dc0_real = real(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
+ dc0_img = imag(dc0_matrix[i157][i162 - 1][jxi468 - 1]);
// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
// real part and 8 bytes to represent the imaginary one.
@@ -1096,7 +1101,8 @@ void ScattererConfiguration::write_formatted(string file_name) {
if (i473 == 1) ici += ies;
fprintf(output, " SPHERE N.%4d\n", i473);
for (int ic472 = 0; ic472 < ici; ic472++) {
- double dc0_real = dc0_matrix[ic472][i473 - 1][0].real(), dc0_img = dc0_matrix[ic472][i473 - 1][0].imag();
+ double dc0_real = real(dc0_matrix[ic472][i473 - 1][0]);
+ double dc0_img = imag(dc0_matrix[ic472][i473 - 1][0]);
fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
}
}
@@ -1110,8 +1116,8 @@ void ScattererConfiguration::write_formatted(string file_name) {
for (int ic477 = 1; ic477 <= ici; ic477++) {
fprintf(output, " NONTRANSITION LAYER N.%2d, SCALE = %3s\n", ic477, reference_variable_name.c_str());
for (int jxi476 = 0; jxi476 < number_of_scales; jxi476++) {
- double dc0_real = dc0_matrix[ic477 - 1][i478 - 1][jxi476].real();
- double dc0_img = dc0_matrix[ic477 - 1][i478 - 1][jxi476].imag();
+ double dc0_real = real(dc0_matrix[ic477 - 1][i478 - 1][jxi476]);
+ double dc0_img = imag(dc0_matrix[ic477 - 1][i478 - 1][jxi476]);
fprintf(output, "%5d %12.4lE%12.4lE\n", (jxi476 + 1), dc0_real, dc0_img);
}
}
diff --git a/src/libnptm/Makefile b/src/libnptm/Makefile
index 61e447c42c97dbce95de53d651657a641e96d30b..24f96311d82573ecf01703372fd4fdeab9dce682 100644
--- a/src/libnptm/Makefile
+++ b/src/libnptm/Makefile
@@ -19,11 +19,11 @@ endif
include ../make.inc
-CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o $(OBJDIR)/lapack_calls.o $(OBJDIR)/algebraic.o
+CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o $(OBJDIR)/lapack_calls.o $(OBJDIR)/algebraic.o $(OBJDIR)/types.o
-CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o $(DYNOBJDIR)/lapack_calls.o $(DYNOBJDIR)/algebraic.o
+CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o $(DYNOBJDIR)/lapack_calls.o $(DYNOBJDIR)/algebraic.o $(DYNOBJDIR)/types.o
-CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g* $(OBJDIR)/lapack_calls.g* $(OBJDIR)/algebraic.g*
+CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g* $(OBJDIR)/lapack_calls.g* $(OBJDIR)/algebraic.g* $(OBJDIR)/types.g*
all: $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index 2df1bc4f6a737ffd330632e17aa14dd6f97c742d..fc6c92af6c77dca0aabce667de72f1e51ef15fec 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -4,13 +4,14 @@
*
* \brief Implementation of the Transition Matrix structure.
*/
-#include <complex.h>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
@@ -30,10 +31,6 @@ typedef __complex__ double dcomplex;
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;
diff --git a/src/libnptm/algebraic.cpp b/src/libnptm/algebraic.cpp
index 9456eaee80e119b1cb6f5f50d40396afa7461a15..48335792275ba7097b29f5b639fe828bbd1aa8d2 100644
--- a/src/libnptm/algebraic.cpp
+++ b/src/libnptm/algebraic.cpp
@@ -4,10 +4,9 @@
*
* \brief Implementation of algebraic functions with different call-backs.
*/
-
-#include <complex.h>
-
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifdef USE_LAPACK
#ifdef USE_MKL
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 377677d897495bf402b2002b7bb585cd3aa38c5c..72cac9914a8d1d61a6f67e93783871476f03ea21 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1,8 +1,13 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clu_subs.cpp
*
* \brief C++ implementation of CLUSTER subroutines.
*/
-#include <complex>
+
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -19,22 +24,22 @@
using namespace std;
void apc(
- double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
- double sqk, double **gapr, std::complex<double> **gapp
+ double ****zpv, int le, dcomplex **am0m, dcomplex **w,
+ double sqk, double **gapr, dcomplex **gapp
) {
- complex<double> **ac, **gap;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
- complex<double> suemp, sueep;
+ dcomplex **ac, **gap;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex uimmp, summ, sume, suem, suee, summp, sumep;
+ dcomplex suemp, sueep;
double cof = 1.0 / sqk;
double cimu = cof / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
- ac = new complex<double>*[nlemt];
- gap = new complex<double>*[3];
- for (int ai = 0; ai < nlemt; ai++) ac[ai] = new complex<double>[2]();
- for (int gi = 0; gi < 3; gi++) gap[gi] = new complex<double>[2]();
+ ac = new dcomplex*[nlemt];
+ gap = new dcomplex*[3];
+ for (int ai = 0; ai < nlemt; ai++) ac[ai] = new dcomplex[2]();
+ for (int gi = 0; gi < 3; gi++) gap[gi] = new dcomplex[2]();
for (int j45 = 1; j45 <= nlemt; j45++) {
int j = j45 - 1;
ac[j][0] = cc0;
@@ -119,9 +124,9 @@ void apc(
sume = gap[0][ipo95 - 1] * cimu;
suee = gap[1][ipo95 - 1] * cof;
suem = gap[2][ipo95 - 1] * cimu;
- gapr[0][ipo95 - 1] = (sume - suem).real();
- gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
- gapr[2][ipo95 - 1] = suee.real();
+ gapr[0][ipo95 - 1] = real(sume - suem);
+ gapr[1][ipo95 - 1] = real((sume + suem) * uim);
+ gapr[2][ipo95 - 1] = real(suee);
sumep = gapp[0][ipo95 - 1] * cimu;
sueep = gapp[1][ipo95 - 1] * cof;
suemp = gapp[2][ipo95 - 1] * cimu;
@@ -137,24 +142,24 @@ void apc(
}
void apcra(
- double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
- double **gaprm, std::complex<double> **gappm
+ double ****zpv, const int le, dcomplex **am0m, int inpol, double sqk,
+ double **gaprm, dcomplex **gappm
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
- complex<double> a11, a12, a21, a22, sum1, sum2, fc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex uimtl, uimtls, ca11, ca12, ca21, ca22;
+ dcomplex a11, a12, a21, a22, sum1, sum2, fc;
double ****svw = new double***[le];
- complex<double> ****svs = new complex<double>***[le];
+ dcomplex ****svs = new dcomplex***[le];
for (int i = 0; i < le; i++) {
svw[i] = new double**[3];
- svs[i] = new complex<double>**[3];
+ svs[i] = new dcomplex**[3];
for (int j = 0; j < 3; j++) {
svw[i][j] = new double*[2];
- svs[i][j] = new complex<double>*[2];
+ svs[i][j] = new dcomplex*[2];
for (int k = 0; k < 2; k++) {
svw[i][j][k] = new double[2]();
- svs[i][j][k] = new complex<double>[2]();
+ svs[i][j][k] = new dcomplex[2]();
}
}
}
@@ -206,7 +211,7 @@ void apcra(
int lmp = l58 + lmpml;
int impmm = lmp * (lmp + 1);
uimtl = uim * (1.0 * lmpml);
- if (lmpml == 0) uimtl = complex<double>(1.0, 0.0);
+ if (lmpml == 0) uimtl = 1.0 + 0.0 * I;
for (int im54 = 1; im54 <= ltpo; im54++) {
int m = im54 - lpo;
int i = imm + m;
@@ -229,7 +234,7 @@ void apcra(
int lsmp = ls + lsmpml;
int ismpmm = lsmp * (lsmp + 1);
uimtls = -uim * (1.0 * lsmpml);
- if (lsmpml == 0) uimtls = complex<double>(1.0, 0.0);
+ if (lsmpml == 0) uimtls = 1.0 + 0.0 * I;
for (int ims = 1; ims <= lstpo; ims++) {
int ms = ims - lspo;
int msmp = ms - mu;
@@ -328,14 +333,14 @@ void apcra(
if (inpol == 0) {
sum1 *= cofs;
sum2 *= cofs;
- gaprm[2][0] = sum1.real();
- gaprm[2][1] = sum1.real();
+ gaprm[2][0] = real(sum1);
+ gaprm[2][1] = real(sum1);
gappm[2][0] = sum2 * uim;
gappm[2][1] = -gappm[2][0];
} else { // label 72
cofs *= 2.0;
- gaprm[2][0] = sum1.real() * cofs;
- gaprm[2][1] = sum2.real() * cofs;
+ gaprm[2][0] = real(sum1) * cofs;
+ gaprm[2][1] = real(sum2) * cofs;
gappm[2][0] = cc0;
gappm[2][1] = cc0;
}
@@ -357,15 +362,12 @@ void apcra(
delete[] svs;
}
-complex<double> cdtp(
- std::complex<double> z, std::complex<double> **am, int i, int jf,
- int k, int nj
-) {
+dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj) {
/* NOTE: the original FORTRAN code treats the AM matrix as a
* vector. This is not directly allowed in C++ and it requires
* accounting for the different dimensions.
*/
- complex<double> result = z;
+ dcomplex result = z;
if (nj > 0) {
int jl = jf + nj - 1;
for (int j = jf; j <= jl; j++) {
@@ -409,9 +411,9 @@ double cgev(int ipamo, int mu, int l, int m) {
return result;
}
-void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
- complex<double> dm, de, cgh, cgk;
- const complex<double> cc0(0.0, 0.0);
+void cms(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
+ dcomplex dm, de, cgh, cgk;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
int ndi = c4->nsph * c4->nlim;
int nbl = 0;
int nsphmo = c4->nsph - 1;
@@ -492,19 +494,19 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
}
void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
- complex<double> ***svf, ***svw, **svs;
- const complex<double> cc0(0.0, 0.0);
- complex<double> cam(0.0, 0.0);
+ dcomplex ***svf, ***svw, **svs;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex cam = cc0;
const int le4po = 4 * c4->le + 1;
- svf = new complex<double>**[le4po];
- svw = new complex<double>**[le4po];
- svs = new complex<double>*[le4po];
+ svf = new dcomplex**[le4po];
+ svw = new dcomplex**[le4po];
+ svs = new dcomplex*[le4po];
for (int si = 0; si < le4po; si++) {
- svf[si] = new complex<double>*[le4po];
- svw[si] = new complex<double>*[4];
- svs[si] = new complex<double>[4]();
- for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new complex<double>[4]();
- for (int sj = 0; sj < 4; sj++) svw[si][sj] = new complex<double>[4]();
+ svf[si] = new dcomplex*[le4po];
+ svw[si] = new dcomplex*[4];
+ svs[si] = new dcomplex[4]();
+ for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new dcomplex[4]();
+ for (int sj = 0; sj < 4; sj++) svw[si][sj] = new dcomplex[4]();
}
double exdc = exri * exri;
double ccs = 1.0 / (vk * vk);
@@ -611,22 +613,22 @@ void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
delete[] svs;
}
-complex<double> ghit(
- int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
- C1_AddOns *c1ao, C4 *c4, C6 *c6
+dcomplex ghit(
+ int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
+ C1_AddOns *c1ao, C4 *c4, C6 *c6
) {
/* NBL identifies transfer vector going from N2 to N1;
* IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
* depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> csum(0.0, 0.0), cfun(0.0, 0.0);
- complex<double> result = cc0;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex csum = cc0, cfun = cc0;
+ dcomplex result = cc0;
if (ihi == 2) {
if (c1->rxx[nbl - 1] == 0.0 && c1->ryy[nbl - 1] == 0.0 && c1->rzz[nbl - 1] == 0.0) {
if (ipamo == 0) {
- if (l1 == l2 && m1 == m2) result = complex<double>(1.0, 0.0);
+ if (l1 == l2 && m1 == m2) result = 1.0 + 0.0 * I;
}
return result;
}
@@ -823,7 +825,7 @@ complex<double> ghit(
}
void hjv(
- double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
+ double exri, double vk, int &jer, int &lcalc, dcomplex &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
) {
int nsphmo = c4->nsph - 1;
@@ -844,7 +846,7 @@ void hjv(
double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
double rr = sqrt(rx * rx + ry * ry + rz * rz);
double rarg = rr * vk * exri;
- arg = complex<double>(rarg, 0.0);
+ arg = rarg + 0.0 * I;
rbf(lit, rarg, lcalc, rfj);
if (lcalc < lit) {
jer = 1;
@@ -862,7 +864,7 @@ void hjv(
for (int lpo38 = 1; lpo38 <= c4->litpo; lpo38++) {
double rpart = rfj[lpo38 - 1];
double ipart = rfn[lpo38 - 1];
- c1ao->vh[lpo38 + ivhb - 1] = complex<double>(rpart, ipart);
+ c1ao->vh[lpo38 + ivhb - 1] = rpart + ipart * I;
}
ivhb += c4->litpo;
} // ns40 loop
@@ -892,23 +894,23 @@ void hjv(
delete[] rfn;
}
-void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier) {
+void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) {
/* NDDMST FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
* STATEMENT.
* N NUMBER OF ROWS IN AM.
* IER IS REPLACED BY 1 FOR SINGULARITY.
*/
double *v = new double[nddmst];
- complex<double> ctemp, cfun;
- complex<double> cc0 = complex<double>(0.0, 0.0);
+ dcomplex ctemp, cfun;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
ier = 0;
int nminus = n - 1;
for (int64_t i = 1; i <= n; i++) {
double sum = 0.0;
for (int64_t j = 1; j <= n; j++) {
sum += (
- am[i - 1][j - 1].real() * am[i - 1][j - 1].real()
- + am[i - 1][j - 1].imag() * am[i - 1][j - 1].imag()
+ real(am[i - 1][j - 1]) * real(am[i - 1][j - 1])
+ + imag(am[i - 1][j - 1]) * imag(am[i - 1][j - 1])
);
} // j1319 loop
v[i - 1] = 1.0 / sum;
@@ -917,7 +919,6 @@ void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier) {
// REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.
// (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS
// ARE PLACED IN V.)
- /* >>> THERE APPEARS TO BE A BUG IN THE FOLLOWING LOOP <<< */
for (int64_t k = 1; k <= n; k++) {
int64_t kplus = k + 1;
int64_t kminus = k - 1;
@@ -927,7 +928,7 @@ void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier) {
cfun = cdtp(-am[i - 1][k - 1], am, i, 1, k, kminus);
ctemp = -cfun;
am[i - 1][k - 1] = ctemp;
- double psq = v[i - 1] * (ctemp.real() * ctemp.real() + ctemp.imag() * ctemp.imag());
+ double psq = v[i - 1] * (real(ctemp) * real(ctemp) + imag(ctemp) * imag(ctemp));
if (psq > psqmax) {
psqmax = psq;
l = i;
@@ -997,15 +998,15 @@ void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier) {
delete[] v;
}
-void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext) {
- double fa11r = fsac[0][0].real();
- double fa11i = fsac[0][0].imag();
- double fa21r = fsac[1][0].real();
- double fa21i = fsac[1][0].imag();
- double fa12r = fsac[0][1].real();
- double fa12i = fsac[0][1].imag();
- double fa22r = fsac[1][1].real();
- double fa22i = fsac[1][1].imag();
+void mextc(double vk, double exri, dcomplex **fsac, double **cextlr, double **cext) {
+ double fa11r = real(fsac[0][0]);
+ double fa11i = imag(fsac[0][0]);
+ double fa21r = real(fsac[1][0]);
+ double fa21i = imag(fsac[1][0]);
+ double fa12r = real(fsac[0][1]);
+ double fa12i = imag(fsac[0][1]);
+ double fa22r = real(fsac[1][1]);
+ double fa22i = imag(fsac[1][1]);
cextlr[0][0] = fa11i * 2.0;
cextlr[0][1] = 0.0;
cextlr[0][2] = -fa12i;
@@ -1048,12 +1049,12 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
}
void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
- const complex<double> cc0(0.0, 0.0);
- complex<double> sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ dcomplex sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
const double exdc = exri * exri;
double ccs = 1.0 / (vk * vk);
double cccs = ccs / exdc;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pi4sq *ccs * ccs);
const int nlemt = c4->nlem + c4->nlem;
@@ -1077,7 +1078,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
am += (c1ao->am0m[i][j] * c1->w[j][ipo18 - 1]);
amp += (c1ao->am0m[i][j] * c1->w[j][jpo - 1]);
} // j10 loop
- sum += (dconjg(am) * am).real();
+ sum += real(dconjg(am) * am);
sump += (dconjg(amp) * am);
sum1 += (dconjg(c1->w[i][ipo18 - 1]) * am);
sum2 += (dconjg(c1->w[i][jpo - 1]) * am);
@@ -1086,7 +1087,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
} // i12 loop
c1ao->scsc[ipo18 - 1] = cccs * sum;
c1ao->scscp[ipo18 - 1] = cccs * sump;
- c1ao->ecsc[ipo18 - 1] = -cccs * sum1.real();
+ c1ao->ecsc[ipo18 - 1] = -cccs * real(sum1);
c1ao->ecscp[ipo18 - 1] = -cccs * sum2;
c1ao->fsac[ipo18 - 1][ipo18 - 1] = csam * sum1;
c1ao->fsac[jpo - 1][ipo18 - 1] = csam * sum2;
@@ -1107,14 +1108,14 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
}
void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> sum1, sum2, sum3, sum4, sumpd;
- complex<double> sums1, sums2, sums3, sums4, csam;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex sum1, sum2, sum3, sum4, sumpd;
+ dcomplex sums1, sums2, sums3, sums4, csam;
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.5 * I;
sum2 = cc0;
sum3 = cc0;
for (int i14 = 1; i14 <= c4->nlem; i14++) { // GPU portable?
@@ -1128,10 +1129,10 @@ void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4)
for (int i16 = 1; i16 <= nlemt; i16++) {
for (int j16 = 1; j16 <= c4->nlem; j16++) {
int je = j16 + c4->nlem;
- double rvalue = (
- dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
- + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
- ).real();
+ double rvalue = real(
+ dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
+ + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
+ );
sumpi += rvalue;
sumpd += (
dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][je - 1]
@@ -1158,16 +1159,16 @@ void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4)
sums4 = cc0;
}
// label 20
- c1ao->ecscm[0] = -cccs * sum2.real();
- c1ao->ecscm[1] = -cccs * sum1.real();
+ c1ao->ecscm[0] = -cccs * real(sum2);
+ c1ao->ecscm[1] = -cccs * real(sum1);
c1ao->ecscpm[0] = -cccs * sum4;
c1ao->ecscpm[1] = -cccs * sum3;
c1ao->fsacm[0][0] = csam * sum2;
c1ao->fsacm[1][0] = csam * sum4;
c1ao->fsacm[1][1] = csam * sum1;
c1ao->fsacm[0][1] = csam * sum3;
- c1ao->scscm[0] = cccs * sums1.real();
- c1ao->scscm[1] = cccs * sums2.real();
+ c1ao->scscm[0] = cccs * real(sums1);
+ c1ao->scscm[1] = cccs * real(sums2);
c1ao->scscpm[0] = cccs * sums3;
c1ao->scscpm[1] = cccs * sums4;
}
@@ -1538,23 +1539,23 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
}
void raba(
- int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
- std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
+ int le, dcomplex **am0m, dcomplex **w, double **tqce,
+ dcomplex **tqcpe, double **tqcs, dcomplex **tqcps
) {
- complex<double> **a, **ctqce, **ctqcs;
- complex<double> acw, acwp, aca, acap, c1, c2, c3;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ dcomplex **a, **ctqce, **ctqcs;
+ dcomplex acw, acwp, aca, acap, c1, c2, c3;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
- a = new complex<double>*[nlemt];
- ctqce = new complex<double>*[2];
- ctqcs = new complex<double>*[2];
- for (int ai = 0; ai < nlemt; ai++) a[ai] = new complex<double>[2]();
+ a = new dcomplex*[nlemt];
+ ctqce = new dcomplex*[2];
+ ctqcs = new dcomplex*[2];
+ for (int ai = 0; ai < nlemt; ai++) a[ai] = new dcomplex[2]();
for (int ci = 0; ci < 2; ci++) {
- ctqce[ci] = new complex<double>[3]();
- ctqcs[ci] = new complex<double>[3]();
+ ctqce[ci] = new dcomplex[3]();
+ ctqcs[ci] = new dcomplex[3]();
}
for (int i20 = 1; i20 <= nlemt; i20++) {
int i = i20 - 1;
@@ -1638,18 +1639,18 @@ void raba(
} // ipo70 loop
for (int ipo78 = 1; ipo78 <= 2; ipo78++) {
int ipo = ipo78 - 1;
- tqce[ipo][0] = (ctqce[ipo][0] - ctqce[ipo][2]).real() * sq2i;
- tqce[ipo][1] = ((ctqce[ipo][0] + ctqce[ipo][2]) * uim).real() * sq2i;
- tqce[ipo][2] = ctqce[ipo][1].real();
+ tqce[ipo][0] = real(ctqce[ipo][0] - ctqce[ipo][2]) * sq2i;
+ tqce[ipo][1] = real((ctqce[ipo][0] + ctqce[ipo][2]) * uim) * sq2i;
+ tqce[ipo][2] = real(ctqce[ipo][1]);
c1 = tqcpe[ipo][0];
c2 = tqcpe[ipo][1];
c3 = tqcpe[ipo][2];
tqcpe[ipo][0] = (c1 - c3) * sq2i;
tqcpe[ipo][1] = (c1 + c3) * (uim * sq2i);
tqcpe[ipo][2] = c2;
- tqcs[ipo][0] = -sq2i * (ctqcs[ipo][0] - ctqcs[ipo][2]).real();
- tqcs[ipo][1] = -sq2i * ((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim).real();
- tqcs[ipo][2] = -1.0 * ctqcs[ipo][1].real();
+ tqcs[ipo][0] = -sq2i * real(ctqcs[ipo][0] - ctqcs[ipo][2]);
+ tqcs[ipo][1] = -sq2i * real((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim);
+ tqcs[ipo][2] = -1.0 * real(ctqcs[ipo][1]);
c1 = tqcps[ipo][0];
c2 = tqcps[ipo][1];
c3 = tqcps[ipo][2];
@@ -1685,12 +1686,12 @@ void rftr(
}
void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
- const complex<double> cc0(0.0, 0.0);
+ const dcomplex cc0 = 0.0 + 0.0 * I;
double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
- complex<double> sum21, rm, re, csam;
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ dcomplex sum21, rm, re, csam;
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
//double scs = 0.0, ecs = 0.0, acs = 0.0;
c3->scs = 0.0;
c3->ecs = 0.0;
@@ -1705,13 +1706,13 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
double fl = 1.0 * (l10 + l10 + 1);
rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
re = 1.0 / c1->rei[l10 - 1][i14 - 1];
- double rvalue = (dconjg(rm) * rm + dconjg(re) * re).real() * fl;
+ double rvalue = real(dconjg(rm) * rm + dconjg(re) * re) * fl;
sums += rvalue;
sum21 += ((rm + re) * fl);
} // l10 loop
sum21 *= -1.0;
double scasec = cccs * sums;
- double extsec = -cccs * sum21.real();
+ double extsec = -cccs * real(sum21);
double abssec = extsec - scasec;
c1->sscs[i14 - 1] = scasec;
c1->sexs[i14 - 1] = extsec;
@@ -1734,11 +1735,11 @@ void scr2(
double vk, double vkarg, double exri, double *duk, C1 *c1, C1_AddOns *c1ao,
C3 *c3, C4 *c4
) {
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
- complex<double> s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
+ dcomplex s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
double ccs = 1.0 / (vk * vk);
- csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ csam = -(ccs / (exri * vk)) * 0.5 * I;
const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pi4sq * ccs * ccs);
cph = uim * exri * vkarg;
@@ -1776,7 +1777,7 @@ void scr2(
c1->sas[i][1][1] = s22 * csam;
}
// label 12
- phas = exp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
+ phas = cexp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
c3->tsas[0][0] += (c1->sas[iogi - 1][0][0] * phas);
c3->tsas[1][0] += (c1->sas[iogi - 1][1][0] * phas);
c3->tsas[0][1] += (c1->sas[iogi - 1][0][1] * phas);
@@ -1814,7 +1815,7 @@ void scr2(
}
void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
- complex<double> *ylm;
+ dcomplex *ylm;
const double pi = acos(-1.0);
c3->gcs = 0.0;
double gcss = 0.0;
@@ -1831,7 +1832,7 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
c3->gcs += gcss;
} // i18 loop
int ylm_size = (c4->litpos > c4->lmtpos) ? c4->litpos : c4->lmtpos;
- ylm = new complex<double>[ylm_size]();
+ ylm = new dcomplex[ylm_size]();
int i = 0;
for (int l1po28 = 1; l1po28 <= c4->lmpo; l1po28++) {
int l1 = l1po28 - 1;
@@ -1901,9 +1902,9 @@ void tqr(
tsk = u[0] * tqsv[0] + u[1] * tqsv[1] + u[2] * tqsv[2];
}
-void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
- complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
- const complex<double> cc0(0.0, 0.0);
+void ztm(dcomplex **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
+ dcomplex gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
int ndi = c4->nsph * c4->nlim;
int i2 = 0;
for (int n2 = 1; n2 <= c4->nsph; n2++) { // GPU portable?
diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp
index 333871e0ae08702f4fcddfef68500ed454a77947..ddf8ddd42f832e06be6a55303d6dc00ea29798d9 100644
--- a/src/libnptm/lapack_calls.cpp
+++ b/src/libnptm/lapack_calls.cpp
@@ -1,8 +1,12 @@
/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
-#include <complex.h>
-
-typedef __complex__ double dcomplex;
+/*! \file lapack_calls.cpp
+ *
+ * \brief Implementation of the interface with LAPACK libraries.
+ */
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifdef USE_LAPACK
#ifdef USE_MKL
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 38d3ab778352797d85366fd1180fd8419f6b4aef..72f8cb1024912d78193dedbb3d71467e8e35de53 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -4,9 +4,9 @@
*
* \brief C++ implementation of SPHERE subroutines.
*/
-#include <complex.h>
-
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -16,10 +16,6 @@ typedef __complex__ double dcomplex;
#include "../include/sph_subs.h"
#endif
-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) {
dcomplex cc0 = 0.0 + 0.0 * I;
dcomplex summ, sume, suem, suee, sum;
@@ -66,7 +62,7 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
}
}
}
- gaps[i40] = sum.real() * cofs;
+ gaps[i40] = real(sum) * cofs;
}
}
@@ -85,7 +81,7 @@ void cbf(int n, dcomplex z, int &nm, dcomplex *csj) {
* point for backward recurrence
* ==========================================================
*/
- double zz = z.real() * z.real() + z.imag() * z.imag();
+ double zz = real(z) * real(z) + imag(z) * imag(z);
double a0 = sqrt(zz);
nm = n;
if (a0 < 1.0e-60) {
@@ -117,8 +113,8 @@ void cbf(int n, dcomplex z, int &nm, dcomplex *csj) {
cf0 = cf1;
cf1 = cf;
}
- double abs_csa = abs(csa);
- double abs_csb = abs(csb);
+ double abs_csa = cabs(csa);
+ double abs_csb = cabs(csb);
if (abs_csa > abs_csb) cs = csa / cf;
else cs = csb / cf0;
for (int k = 0; k <= nm; k++) {
@@ -196,7 +192,7 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
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);
+ c2->dlri[kpnt] = 0.0 + 0.0 * I;
const int i90 = i - 1;
const int ns90 = ns - 1;
const int ic90 = ic - 1;
@@ -232,7 +228,7 @@ void dme(
arg = vkr1 * c2->vkt[i - 1];
arin = arg;
bool goto32 = false;
- if (arg.imag() != 0.0) {
+ if (imag(arg) != 0.0) {
cbf(lipo, arg, lcalc, cfj);
if (lcalc < lipo) {
jer = 5;
@@ -244,7 +240,7 @@ void dme(
goto32 = true;
}
if (!goto32) {
- rbf(lipo, arg.real(), lcalc, rfj);
+ rbf(lipo, real(arg), lcalc, rfj);
if (lcalc < lipo) {
jer = 5;
delete[] rfj;
@@ -359,7 +355,7 @@ double envj(int n, double x) {
return result;
}
-void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
+void mmulc(dcomplex *vint, double **cmullr, double **cmul) {
double sm2 = real(vint[0]);
double s24 = real(vint[1]);
double d24 = imag(vint[1]);
@@ -587,9 +583,9 @@ void rabas(
int l = l70 + 1;
double fl = 1.0 * (l + l + 1);
dcomplex rm = 1.0 / c1->rmi[l70][i80];
- double rmm = (rm * dconjg(rm)).real();
+ double rmm = real(rm * dconjg(rm));
dcomplex re = 1.0 / c1->rei[l70][i80];
- double rem = (re * dconjg(re)).real();
+ double rem = real(re * dconjg(re));
if (inpol == 0) {
dcomplex pce = ((rm + re) * uim) * fl;
dcomplex pcs = ((rmm + rem) * fl) * uim;
@@ -598,7 +594,7 @@ void rabas(
tqsps[0][i80] -= pcs;
tqsps[1][i80] += pcs;
} else {
- double ce = (rm + re).real() * fl;
+ double ce = real(rm + re) * fl;
double cs = (rmm + rem) * fl;
tqse[0][i80] -= ce;
tqse[1][i80] += ce;
@@ -874,7 +870,7 @@ void sphar(
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 * (cosrmp[m - 1] - sinrmp[m - 1 * I]);
+ ylm[lmy - 1] = save * (cosrmp[m - 1] - sinrmp[m - 1] * I);
label45:
if (m >= l) goto label47;
m += 1;
@@ -906,12 +902,12 @@ void sscr0(dcomplex &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
re = 1.0 / c1->rei[l10][i12];
dcomplex rm_cnjg = dconjg(rm);
dcomplex re_cnjg = dconjg(re);
- sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
+ sums += real(rm_cnjg * rm + re_cnjg * re) * fl;
sum21 += (rm + re) * fl;
}
sum21 *= -1.0;
double scasec = cccs * sums;
- double extsec = -cccs * sum21.real();
+ double extsec = -cccs * real(sum21);
double abssec = extsec - scasec;
c1->sscs[i12] = scasec;
c1->sexs[i12] = extsec;
@@ -970,7 +966,7 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
int j = 0;
for (int ipo1 = 0; ipo1 < 2; ipo1++) {
for (int jpo1 = 0; jpo1 < 2; jpo1++) {
- complex<double> cc = dconjg(c1->sas[i24][jpo1][ipo1]);
+ dcomplex cc = dconjg(c1->sas[i24][jpo1][ipo1]);
for (int ipo2 = 0; ipo2 < 2; ipo2++) {
for (int jpo2 = 0; jpo2 < 2; jpo2++) {
c1->vints[i24][j++] = c1->sas[i24][jpo2][ipo2] * cc * cfsq;
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
index 811264d43c8f2b8370e60ee1ad854f93c7c69f45..a1a8e1aecf19841471964057b2d3442711638e82 100644
--- a/src/libnptm/tfrfme.cpp
+++ b/src/libnptm/tfrfme.cpp
@@ -5,13 +5,14 @@
* \brief Implementation of the trapping calculation objects.
*/
-#include <complex.h>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
@@ -31,16 +32,12 @@ typedef __complex__ double dcomplex;
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;
nlmmt = 2 * lm * (lm + 2);
const int size = nkv * nkv * nlmmt;
- wk = new complex<double>[size]();
+ wk = new dcomplex[size]();
last_index = 0;
}
@@ -691,7 +688,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
TFRFME* TFRFME::from_legacy(string file_name) {
fstream input;
TFRFME *instance = NULL;
- complex<double> **_wsum = NULL;
+ dcomplex **_wsum = NULL;
double *coord_vec = NULL;
input.open(file_name.c_str(), ios::in | ios::binary);
if (input.is_open()) {
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index 1740872f1d70f5b300df3c23eb7436f6885dc8f5..2591b8c335a245a64cd41ae26ad50e6720810ace 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -5,10 +5,11 @@
* \brief C++ implementation of TRAPPING subroutines.
*/
#include <cmath>
-#include <complex.h>
#include <fstream>
-typedef __complex__ double dcomplex;
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
@@ -26,10 +27,6 @@ typedef __complex__ double dcomplex;
#include "../include/tra_subs.h"
#endif
-double real(dcomplex z) { return __real__ z; }
-
-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++) {
@@ -70,7 +67,7 @@ void ffrf(
) {
const dcomplex cc0 = 0.0 + 0.0 * I;
const dcomplex uim = 0.0 + 1.0 * I;
- complex<double> uimmp, summ, sume, suem, suee;
+ dcomplex uimmp, summ, sume, suem, suee;
dcomplex *gap = new dcomplex[3]();
for (int imu50 = 1; imu50 <= 3; imu50++) {
@@ -232,7 +229,7 @@ void ffrt(dcomplex *ac, dcomplex *ws, double *ffte, double *ffts, CIL *cil) {
}
} // im60 loop
} // l60 loop
- ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
+ ffte[0] = real(ctqce[0] - ctqce[2]) * sq2i;
ffte[1] = real(sq2iti * (ctqce[0] + ctqce[2]));
ffte[2] = real(ctqce[1]);
ffts[0] = -sq2i * real(ctqcs[0] - ctqcs[2]);
@@ -364,7 +361,7 @@ void wamff(
dcomplex **w, *ylm;
const dcomplex cc0 = 0.0 + 0.0 * I;
const dcomplex uim = 0.0 + 1.0 * I;
- complex<double> cfam, cf1, cf2;
+ dcomplex cfam, cf1, cf2;
double rho, cph, sph, cth, sth, r;
double ftc1, ftc2;
double *up = new double[3];
@@ -469,7 +466,7 @@ void wamff(
skip62 = false;
}
} else { // label 57
- double fam = tra * cexp(-apfa * apfa) / sqrt(cth);
+ double fam = tra * exp(-apfa * apfa) / sqrt(cth);
if (lmode == 0) {
double f1 = cph * fam;
double f2 = -sph * fam;
diff --git a/src/sphere/np_sphere.cpp b/src/sphere/np_sphere.cpp
index e3033dae6fa36f373e17c1fda92f802537ce3b53..7f381e41ee28b645b53376828ddf33c4c68b7a88 100644
--- a/src/sphere/np_sphere.cpp
+++ b/src/sphere/np_sphere.cpp
@@ -1,3 +1,5 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file np_sphere.cpp
*
* \brief Single sphere scattering problem handler.
@@ -15,10 +17,13 @@
* while the advanced mode implements the possibility to change input and
* output options, without having to modify the code.
*/
-#include <complex>
#include <cstdio>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 7016d8361eb8d61fef8437f504516352f0444948..b4bf9652d3e0a9401e076f750b92615073750b97 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -1,13 +1,18 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file sphere.cpp
*
* \brief Implementation of the single sphere calculation.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
@@ -37,7 +42,7 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void sphere(string config_file, string data_file, string output_path) {
- complex<double> arg, s0, tfsas;
+ dcomplex arg, s0, tfsas;
double th, ph;
printf("INFO: making legacy configuration...\n");
ScattererConfiguration *sconf = NULL;
@@ -57,7 +62,7 @@ void sphere(string config_file, string data_file, string output_path) {
} catch(const OpenConfigurationFileException &ex) {
printf("\nERROR: failed to open geometry configuration file.\n");
printf("FILE: %s\n", ex.what());
- if (sconf) delete sconf;
+ if (sconf != NULL) delete sconf;
exit(1);
}
if (sconf->number_of_spheres == gconf->number_of_spheres) {
@@ -83,21 +88,21 @@ void sphere(string config_file, string data_file, string output_path) {
cmul[i] = new double[4];
cmullr[i] = new double[4];
}
- complex<double> **tqspe, **tqsps;
+ dcomplex **tqspe, **tqsps;
double **tqse, **tqss;
tqse = new double*[2];
tqss = new double*[2];
- tqspe = new std::complex<double>*[2];
- tqsps = new std::complex<double>*[2];
+ tqspe = new dcomplex*[2];
+ tqsps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[2]();
tqss[ti] = new double[2]();
- tqspe[ti] = new std::complex<double>[2]();
- tqsps[ti] = new std::complex<double>[2]();
+ tqspe[ti] = new dcomplex[2]();
+ tqsps[ti] = new dcomplex[2]();
}
double frx = 0.0, fry = 0.0, frz = 0.0;
double cfmp, cfsp, sfmp, sfsp;
- complex<double> *vint = new complex<double>[16];
+ dcomplex *vint = new dcomplex[16];
int jw;
int nsph = gconf->number_of_spheres;
C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
@@ -244,8 +249,8 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " \n");
}
for (int jxi488 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
- printf("INFO: running scale iteration...\n");
int jxi = jxi488 + 1;
+ printf("INFO: running scale iteration %d of %d...\n", jxi, sconf->number_of_scales);
fprintf(output, "========== JXI =%3d ====================\n", jxi);
double xi = sconf->scale_vec[jxi488];
if (sconf->idfc >= 0) {
@@ -289,7 +294,7 @@ void sphere(string config_file, string data_file, string output_path) {
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
- fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, arg.real(), arg.imag());
+ fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, real(arg), imag(arg));
tppoan.close();
fclose(output);
return;
@@ -323,8 +328,8 @@ void sphere(string config_file, string data_file, string output_path) {
output,
" SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
c2->vsz[i170],
- c2->vkt[i170].real(),
- c2->vkt[i170].imag()
+ real(c2->vkt[i170]),
+ imag(c2->vkt[i170])
);
}
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
@@ -341,12 +346,12 @@ void sphere(string config_file, string data_file, string output_path) {
c1->sqscs[i170], c1->sqabs[i170],
c1->sqexs[i170]
);
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i170]), imag(c1->fsas[i170]));
double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
s0 = c1->fsas[i170] * exri;
- double qschu = csch * s0.imag();
- double pschu = csch * s0.real();
- double s0mag = cs0 * abs(s0);
+ double qschu = csch * imag(s0);
+ double pschu = csch * real(s0);
+ double s0mag = cs0 * cabs(s0);
fprintf(
output,
" QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
@@ -359,33 +364,33 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
- double val = tqspe[0][i170].real();
+ double val = real(tqspe[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[0][i170].imag();
+ val = imag(tqspe[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].real();
+ val = real(tqsps[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].imag();
+ val = imag(tqsps[0][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
- val = tqspe[1][i170].real();
+ val = real(tqspe[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[1][i170].imag();
+ val = imag(tqspe[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].real();
+ val = real(tqsps[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].imag();
+ val = imag(tqsps[1][i170]);
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
} // End if iog[i170] >= i
} // i170 loop
if (nsph != 1) {
- fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
+ fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", real(tfsas), imag(tfsas));
double csch = 2.0 * vk * sqsfi / gcs;
s0 = tfsas * exri;
- double qschu = csch * s0.imag();
- double pschu = csch * s0.real();
- double s0mag = cs0 * abs(s0);
+ double qschu = csch * imag(s0);
+ double pschu = csch * real(s0);
+ double s0mag = cs0 * cabs(s0);
fprintf(
output,
" QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
@@ -491,13 +496,13 @@ void sphere(string config_file, string data_file, string output_path) {
fprintf(output, " SPHERE %2d\n", ns);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
- c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
+ real(c1->sas[ns226][0][0]), imag(c1->sas[ns226][0][0]),
+ real(c1->sas[ns226][1][0]), imag(c1->sas[ns226][1][0])
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
- c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
+ real(c1->sas[ns226][0][1]), imag(c1->sas[ns226][0][1]),
+ real(c1->sas[ns226][1][1]), imag(c1->sas[ns226][1][1])
);
if (jths == 1 && jphs == 1)
fprintf(
@@ -523,9 +528,9 @@ void sphere(string config_file, string data_file, string output_path) {
else fprintf(output, "\n");
}
for (int vi = 0; vi < 16; vi++) {
- double value = vint[vi].real();
+ double value = real(vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = vint[vi].imag();
+ value = imag(vint[vi]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int imul = 0; imul < 4; imul++) {
diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp
index 98a681d3895a63659978cd130352083a7028c6e7..b790afd907df1dbb229e0f0ac4b214fd42ea7d5b 100644
--- a/src/testing/test_TEDF.cpp
+++ b/src/testing/test_TEDF.cpp
@@ -1,11 +1,16 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
//! \file test_TEDF.cpp
-#include <complex>
#include <cstdio>
#include <exception>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
diff --git a/src/testing/test_TTMS.cpp b/src/testing/test_TTMS.cpp
index 7c55d16a7939a67ba211f90644b23ba2d329852e..77b8d7c22b1fdc59a73a7165ae2369c2e3c368ff 100644
--- a/src/testing/test_TTMS.cpp
+++ b/src/testing/test_TTMS.cpp
@@ -1,11 +1,16 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
//! \file test_TTMS.cpp
-#include <complex>
#include <cstdio>
#include <exception>
#include <hdf5.h>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index c4866e1f65073bcc2f61ec54d7a23a125ac315ca..2cd6ab9f33a227aff22d4255d4cdf48d01ccdb08 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -1,14 +1,19 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file cfrfme.cpp
*
* C++ implementation of FRFME functions.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif
@@ -47,10 +52,10 @@ void frfme(string data_file, string output_path) {
Swap2 *tt2 = NULL;
char namef[7];
char more;
- complex<double> **w = NULL;
- complex<double> *wk = NULL;
- const complex<double> cc0(0.0, 0.0);
- const complex<double> uim(0.0, 1.0);
+ dcomplex **w = NULL;
+ dcomplex *wk = NULL;
+ const dcomplex cc0 = 0.0 + 0.0 * I;
+ const dcomplex uim = 0.0 + 1.0 * I;
int line_count = 0, last_read_line = 0;
regex re = regex("-?[0-9]+");
string *file_lines = load_file(data_file, &line_count);
@@ -255,7 +260,7 @@ void frfme(string data_file, string output_path) {
double *_xv = tfrfme->get_x();
double *_yv = tfrfme->get_y();
double *_zv = tfrfme->get_z();
- complex<double> **_wsum = tfrfme->get_matrix();
+ dcomplex **_wsum = tfrfme->get_matrix();
for (int i24 = nxshpo; i24 <= nxs; i24++) {
_xv[i24] = _xv[i24 - 1] + delxyz;
_xv[nxv - i24 - 1] = -_xv[i24];
@@ -307,15 +312,15 @@ void frfme(string data_file, string output_path) {
tt2->set_param("nrvc", 1.0 * nrvc);
tt2->write_binary(temp_name2, "HDF5");
for (int j80 = jlmf; j80 <= jlml; j80++) {
- complex<double> *tt1_wk = tt1->get_vector();
+ dcomplex *tt1_wk = tt1->get_vector();
int wk_index = 0;
// w matrix
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
- w = new complex<double>*[nkv];
- for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
+ w = new dcomplex*[nkv];
+ for (int wi = 0; wi < nkv; wi++) w[wi] = new dcomplex[nkv]();
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
@@ -331,19 +336,19 @@ void frfme(string data_file, string output_path) {
for (int ix65 = 0; ix65 < nxv; ix65++) {
double x = _xv[ix65];
ixyz++;
- complex<double> sumy = cc0;
+ dcomplex sumy = cc0;
for (int jy60 = 0; jy60 < nkv; jy60++) {
double vky = vkv[jy60];
double vkx = vkv[nkv - 1];
double vkzf = vkzm[0][jy60];
- complex<double> phasf = exp(uim * (-vkx * x + vky * y + vkzf * z));
+ dcomplex phasf = cexp(uim * (-vkx * x + vky * y + vkzf * z));
double vkzl = vkzm[nkv - 1][jy60];
- complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));
- complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
+ dcomplex phasl = cexp(uim * (vkx * x + vky * y + vkzl * z));
+ dcomplex sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
for (int jx55 = 2; jx55 <= nks; jx55++) {
vkx = vkv[jx55 - 1];
double vkz = vkzm[jx55 - 1][jy60];
- complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
+ dcomplex phas = cexp(uim * (vkx * x + vky * y + vkz * z));
sumx += (w[jx55 - 1][jy60] * phas);
} // jx55 loop
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index 7e3619a252ada2d7102a4d88b3b2a2ba97e60f19..a8033d8098c08694464fe1ab1ac206c84fd18f04 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -1,14 +1,19 @@
+/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
+
/*! \file clffft.cpp
*
* \brief C++ implementation of LFFFT functions.
*/
-#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif
@@ -37,18 +42,18 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void lffft(string data_file, string output_path) {
- const complex<double> uim(0.0, 1.0);
+ const dcomplex uim = 0.0 + 1.0 * I;
const double sq2i = 1.0 / sqrt(2.0);
- const complex<double> sq2iti = sq2i * uim;
+ const dcomplex sq2iti = sq2i * uim;
fstream tlfff, tlfft;
double ****zpv = NULL;
- complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;
- complex<double> **am0m = NULL;
- complex<double> **amd = NULL;
+ dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;
+ dcomplex **am0m = NULL;
+ dcomplex **amd = NULL;
int **indam = NULL;
- complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
- complex<double> **_wsum = NULL;
+ dcomplex *tmsm = NULL, *tmse = NULL, **tms = NULL;
+ dcomplex **_wsum = NULL;
int jft, jss, jtw;
int is, le, nvam = 0;
double vks, exris;
@@ -79,55 +84,55 @@ void lffft(string data_file, string output_path) {
cil->nlem = le * (le + 2);
cil->nlemt = cil->nlem + cil->nlem;
if (is >= 2222) { // label 120
- tms = new complex<double>*[le];
- for (int ti = 0; ti < le; ti++) tms[ti] = new complex<double>[3]();
+ tms = new dcomplex*[le];
+ for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[3]();
// QUESTION|WARNING: original code uses LM without defining it. Where does it come from?
int lm = le;
for (int i = 0; i < lm; i++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][0] = complex<double>(vreal, vimag);
+ tms[i][0] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][1] = complex<double>(vreal, vimag);
+ tms[i][1] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tms[i][2] = complex<double>(vreal, vimag);
+ tms[i][2] = vreal + vimag * I;
} // i loop
} else if (is >= 1111) { // label 125
- tmsm = new complex<double>[le]();
- tmse = new complex<double>[le]();
+ tmsm = new dcomplex[le]();
+ tmse = new dcomplex[le]();
for (int i = 0; i < le; i++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tmsm[i] = complex<double>(vreal, vimag);
+ tmsm[i] = vreal + vimag * I;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- tmse[i] = complex<double>(vreal, vimag);
+ tmse[i] = vreal + vimag * I;
} // i loop
} else if (is >= 0) { // label 135
- am0m = new complex<double>*[cil->nlemt];
- for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new complex<double>[cil->nlemt]();
+ am0m = new dcomplex*[cil->nlemt];
+ for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new dcomplex[cil->nlemt]();
for (int i = 0; i < cil->nlemt; i++) {
for (int j = 0; j < cil->nlemt; j++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- am0m[i][j] = complex<double>(vreal, vimag);
+ am0m[i][j] = vreal + vimag * I;
} // j loop
} // i loop
} else if (is < 0) {
nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;
- amd = new complex<double>*[nvam];
- for (int ai = 0; ai < nvam; ai++) amd[ai] = new complex<double>[4]();
+ amd = new dcomplex*[nvam];
+ for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[4]();
for (int i = 0; i < nvam; i++) {
for (int j = 0; j < 4; j++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- amd[i][j] = complex<double>(vreal, vimag);
+ amd[i][j] = vreal + vimag * I;
} // j loop
} // i loop
indam = new int*[le];
@@ -265,8 +270,8 @@ void lffft(string data_file, string output_path) {
// label 160
const int nlmm = lm * (lm + 2);
const int nlmmt = nlmm + nlmm;
- ws = new complex<double>[nlmmt]();
- if (lm > le) wsl = new complex<double>[nlmmt]();
+ ws = new dcomplex[nlmmt]();
+ if (lm > le) wsl = new dcomplex[nlmmt]();
// FORTRAN writes two output formatted files without opening them
// explicitly. It is assumed thay can be opened here.
string out66_name = output_path + "/c_out66.txt";
@@ -282,12 +287,10 @@ void lffft(string data_file, string output_path) {
//binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
int row = i;
int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
- complex<double> value = _wsum[row][col];
+ dcomplex value = _wsum[row][col];
if (lm <= le) {
- //ws[i] = complex<double>(vreal, vimag);
ws[i] = value;
} else { // label 170
- //wsl[i] = complex<double>(vreal, vimag);
wsl[i] = value;
for (int i175 = 0; i175 < cil->nlem; i175++) {
int ie = i175 + cil->nlem;
@@ -301,21 +304,21 @@ void lffft(string data_file, string output_path) {
if (is != 2222) {
if (is != 1111) {
if (is > 0) { // Goes to 305
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
camp(ac, am0m, ws, cil);
// Goes to 445
} else if (is < 0) { // Goes to 405
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
czamp(ac, amd, indam, ws, cil);
// Goes to 445
}
} else {
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
samp(ac, tmsm, tmse, ws, cil);
// Goes to 445
}
} else {
- ac = new complex<double>[cil->nlemt]();
+ ac = new dcomplex[cil->nlemt]();
sampoa(ac, tms, ws, cil);
// Goes to 445
}