From 18a207ce8138ea6422b26cf921d08a5ba851144b Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Fri, 5 Jan 2024 17:48:21 +0100
Subject: [PATCH] Finish migration of trapping to C++
---
src/include/tra_subs.h | 76 +++++++++++++++++++++++++++++--
src/trapping/lffft.cpp | 100 ++++++++++++++++++++++-------------------
2 files changed, 126 insertions(+), 50 deletions(-)
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index b100fcd6..5eb99136 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -301,7 +301,6 @@ void camp(
std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
CIL *cil
) {
- ac = new std::complex<double>[cil->nlemt]();
for (int j = 0; j < cil->nlemt; j++) {
for (int i = 0; i < cil->nlemt; i++) {
ac[j] += (am0m[j][i] * ws[i]);
@@ -326,7 +325,6 @@ void czamp(
const std::complex<double> cc0(0.0, 0.0);
const std::complex<double> uim(0.0, 1.0);
std::complex<double> summ, sume;
- ac = new std::complex<double>[cil->nlemt]();
for (int im20 = 1; im20 <= cil->mxim; im20++) {
int m = im20 - cil->mxmpo;
int abs_m = (m < 0) ? -m : m;
@@ -477,6 +475,78 @@ void ffrf(
delete[] gap;
}
+/*! 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.
+ * \param ccr: `CCR *` Pointer to a CCR structure.
+ */
+void ffrt(
+ std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
+ CIL *cil
+) {
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ const double sq2i = 1.0 / sqrt(2.0);
+ const std::complex<double> sq2iti = uim * sq2i;
+ std::complex<double> aca, acw;
+ std::complex<double> *ctqce, *ctqcs;
+
+ ctqce = new std::complex<double>[3]();
+ ctqcs = new std::complex<double>[3]();
+ for (int l60 = 1; l60 < cil->le; l60++) {
+ int lpo = l60 + 1;
+ int il = l60 * lpo;
+ int ltpo = l60 + lpo;
+ for (int im60 = 1; im60 <= ltpo; im60++) {
+ double rmu;
+ int m = im60 - lpo;
+ int i = m + il;
+ int ie = i + cil->nlem;
+ int mmmu = m + 1;
+ int mmmmu = (mmmu < 0) ? -mmmu: mmmu;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + cil->nlem;
+ rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
+ acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
+ aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
+ ctqce[0] += (acw * rmu);
+ ctqcs[0] += (aca * rmu);
+ }
+ // label 30
+ rmu = -1.0 * m;
+ acw = dconjg(ac[i - 1]) * ws[i - 1] + dconjg(ac[ie - 1]) * ws[ie - 1];
+ aca = dconjg(ac[i - 1]) * ac[i - 1] + dconjg(ac[ie - 1]) * ac[ie - 1];
+ ctqce[1] += (acw * rmu);
+ ctqcs[1] += (aca * rmu);
+ mmmu = m - 1;
+ mmmmu = (mmmu < 0) ? -mmmu: mmmu;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + cil->nlem;
+ rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
+ acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
+ aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
+ ctqce[2] += (acw * rmu);
+ ctqcs[2] += (aca * rmu);
+ }
+ } // 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();
+
+ delete[] ctqce;
+ delete[] ctqcs;
+}
+
/*! C++ porting of SAMP
*
* \param ac: Vector of complex. QUESTION: definition?
@@ -492,7 +562,6 @@ void samp(
std::complex<double> *ws, CIL *cil
) {
int i = 0;
- ac = new std::complex<double>[cil->nlemt]();
for (int l20 = 0; l20 < cil->le; l20++) {
int l = l20 + 1;
int ltpo = l + l + 1;
@@ -522,7 +591,6 @@ void sampoa(
tm[0] = new std::complex<double>[2]();
tm[1] = new std::complex<double>[2]();
int i = 0;
- ac = new std::complex<double>[cil->nlemt]();
for (int l20 = 0; l20 < cil->le; l20++) {
tm[0][0] = tms[l20][0];
tm[0][1] = tms[l20][1];
diff --git a/src/trapping/lffft.cpp b/src/trapping/lffft.cpp
index 4b8fb7ab..20bb1aeb 100644
--- a/src/trapping/lffft.cpp
+++ b/src/trapping/lffft.cpp
@@ -262,61 +262,69 @@ int main() {
ws[ie] = wsl[iel];
} // i175 loop
}
- } // i loop
- // label 180
- bool goto475 = false;
- if (is != 2222) {
- if (is != 1111) {
- if (is > 0) { // Goes to 305
- camp(ac, am0m, ws, cil);
- } else if (is < 0) { // Goes to 405
- czamp(ac, amd, indam, ws, cil);
+ // label 180
+ if (is != 2222) {
+ if (is != 1111) {
+ if (is > 0) { // Goes to 305
+ ac = new complex<double>[cil->nlemt]();
+ camp(ac, am0m, ws, cil);
+ // Goes to 445
+ } else if (is < 0) { // Goes to 405
+ ac = new complex<double>[cil->nlemt]();
+ czamp(ac, amd, indam, ws, cil);
+ // Goes to 445
+ }
+ } else {
+ ac = new complex<double>[cil->nlemt]();
+ samp(ac, tmsm, tmse, ws, cil);
+ // Goes to 445
}
} else {
- samp(ac, tmsm, tmse, ws, cil);
+ ac = new complex<double>[cil->nlemt]();
+ sampoa(ac, tms, ws, cil);
+ // Goes to 445
}
- } else {
- sampoa(ac, tms, ws, cil);
- }
- // label 445
- if (jft <= 0) {
- double *fffe = new double[3]();
- double *fffs = new double[3]();
- ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
- if (jss == 1) {
- // Writes to 66
- } else { // label 450
- for (int i = 0; i < 3; i++) {
- double value = fffe[i] - fffs[i];
- tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));
- }
- if (jtw == 1) {
+ bool goto475 = false;
+ // label 445
+ if (jft <= 0) {
+ double *fffe = new double[3]();
+ double *fffs = new double[3]();
+ ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
+ if (jss == 1) {
// Writes to 66
+ } else { // label 450
+ for (int i = 0; i < 3; i++) {
+ double value = fffe[i] - fffs[i];
+ tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ if (jtw == 1) {
+ // Writes to 66
+ }
}
+ if (jft < 0) goto475 = true;
+ delete[] fffe;
+ delete[] fffs;
}
- if (jft < 0) goto475 = true;
- delete[] fffe;
- delete[] fffs;
- }
- // label 460
- if (!goto475) {
- double *ffte = new double[3]();
- double *ffts = new double[3]();
- // Would call FFRT(AC,WS,FFTE,FFTS)
- if (jss == 1) {
- // Writes to 67
- } else { // label 470
- for (int i = 0; i < 3; i++) {
- double value = ffte[i] - ffts[i];
- tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));
- }
- if (jtw == 1) {
+ // label 460
+ if (!goto475) {
+ double *ffte = new double[3]();
+ double *ffts = new double[3]();
+ ffrt(ac, ws, ffte, ffts, cil);
+ if (jss == 1) {
// Writes to 67
+ } else { // label 470
+ for (int i = 0; i < 3; i++) {
+ double value = ffte[i] - ffts[i];
+ tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ if (jtw == 1) {
+ // Writes to 67
+ }
}
+ delete[] ffte;
+ delete[] ffts;
}
- delete[] ffte;
- delete[] ffts;
- }
+ } // i loop
} // ix475 loop
} // iy475 loop
} // iz475 loop
--
GitLab