diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index bf3b108a9d2ada949446800aa4d09c6721c2cbfe..b100fcd6f08be3a701f7d950486e88656af0e036 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -17,6 +17,7 @@ struct CCR {
//End of TRAPPING structures
//Externally defined subroutines
+extern double cg1(int lmpml, int mu, int l, int m);
extern std::complex<double> dconjg(std::complex<double> z);
extern void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
@@ -287,6 +288,223 @@ void frfmer(
delete[] wk;
}
+/*! C++ porting of CAMP
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param am0m: Matrix of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ *
+ * This function builds the AC vector using AM0M and WS.
+ */
+void camp(
+ std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
+ CIL *cil
+) {
+ 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]);
+ } // i loop
+ } // j loop
+}
+
+/*! C++ porting of CZAMP
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param amd: Matrix of complex. QUESTION: definition?
+ * \param indam: `int **`. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ *
+ * This function builds the AC vector using AMD, INDAM and WS.
+ */
+void czamp(
+ std::complex<double> *ac, std::complex<double> **amd, int **indam,
+ std::complex<double> *ws, CIL *cil
+) {
+ 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;
+ int lmn = (abs_m > 1) ? abs_m : 1;
+ for (int l20 = lmn; l20 <= cil->le; l20++) {
+ int i = m + l20 * (l20 + 1);
+ int ie = i + cil->nlem;
+ summ = cc0;
+ sume = cc0;
+ for (int ls15 = lmn; ls15 <= cil->le; ls15++) {
+ int is = m + ls15 * (ls15 + 15) - 1;
+ int ise = is + cil->nlem;
+ int num = indam[l20 - 1][ls15 - 1] + m - 1;
+ summ += (amd[num][0] * ws[is] + amd[num][1] * ws[ise]);
+ sume += (amd[num][2] * ws[is] + amd[num][3] * ws[ise]);
+ } // ls15 loop
+ ac[i - 1] = summ;
+ ac[ie - 1] = sume;
+ } // l20 loop
+ } // im20 loop
+}
+
+/*! C++ porting of FFRF
+ *
+ * \param zpv: `double ****`. QUESTION: definition?
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param fffe: `double *`. QUESTION: definition?
+ * \param fffs: `double *`. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ * \param ccr: `CCR *` Pointer to a CCR structure.
+ */
+void ffrf(
+ double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
+ double *fffs, CIL *cil, CCR *ccr
+) {
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ //const std::complex<double> sq2iti = uim / sqrt(2.0);
+ std::complex<double> uimmp, summ, sume, suem, suee;
+ std::complex<double> *gap = new std::complex<double>[3]();
+
+ for (int imu50 = 1; imu50 <= 3; imu50++) {
+ int mu = imu50 - 2;
+ gap[imu50 - 1] = cc0;
+ for (int l40 = 1; l40 <= cil->le; l40++) {
+ int lpo = l40 + 1;
+ int ltpo = lpo + l40;
+ int imm = l40 * lpo;
+ for (int ilmp40 = 1; ilmp40 <= 3; ilmp40++) {
+ if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) continue; // ilmp40 loop
+ int lmpml = ilmp40 - 2;
+ int lmp = l40 + lmpml;
+ uimmp = uim * (-1.0 * lmpml);
+ int impmmmp = lmp * (lmp + 1);
+ for (int im30 = 1; im30 <= ltpo; im30++) {
+ int m = im30 - lpo;
+ int mmp = m - mu;
+ int abs_mmp = (mmp < 0) ? -mmp : mmp;
+ if (abs_mmp <= lmp) {
+ int i = imm + m;
+ int ie = i + cil->nlem;
+ int imp = impmmmp + mmp;
+ int impe = imp + cil->nlem;
+ double cgc = cg1(lmpml, mu, l40, m);
+ summ = dconjg(ws[i - 1]) * ac[imp - 1];
+ sume = dconjg(ws[i - 1]) * ac[impe - 1];
+ suem = dconjg(ws[ie - 1]) * ac[imp - 1];
+ suee = dconjg(ws[ie - 1]) * ac[impe - 1];
+ if (lmpml != 0) {
+ summ *= uimmp;
+ sume *= uimmp;
+ suem *= uimmp;
+ suee *= uimmp;
+ }
+ // label 25
+ gap[imu50 - 1] += (cgc * (
+ summ * zpv[l40 - 1][ilmp40 - 1][0][0] +
+ sume * zpv[l40 - 1][ilmp40 - 1][0][1] +
+ suem * zpv[l40 - 1][ilmp40 - 1][1][0] +
+ suee * zpv[l40 - 1][ilmp40 - 1][1][1]
+ )
+ );
+ }
+ } // im30 loop
+ } // ilmp40
+ } // l40 loop
+ } // imu50 loop
+ sume = -gap[0] * ccr->cimu;
+ suee = -gap[1] * ccr->cof;
+ suem = -gap[2] * ccr->cimu;
+ fffe[0] = (sume - suem).real();
+ fffe[1] = ((sume + suem) * uim).real();
+ fffe[2] = suee.real();
+
+ for (int imu90 = 1; imu90 <= 3; imu90++) {
+ int mu = imu90 - 2;
+ gap[imu90 - 1] = cc0;
+ for (int l80 = 1; l80 <= cil->le; l80++) {
+ int lpo = l80 + 1;
+ int ltpo = lpo + l80;
+ int imm = l80 * lpo;
+ for (int ilmp80 = 1; ilmp80 <= 3; ilmp80++) {
+ if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) continue; // ilmp80 loop
+ int lmpml = ilmp80 - 2;
+ int lmp = l80 + lmpml;
+ uimmp = uim * (-1.0 * lmpml);
+ int impmmmp = lmp * (lmp + 1);
+ for (int im70 = 1; im70 <= ltpo; im70++) {
+ int m = im70 - lpo;
+ int mmp = m - mu;
+ int abs_mmp = (mmp < 0) ? -mmp : mmp;
+ if (abs_mmp <= lmp) {
+ int i = imm + m;
+ int ie = i + cil->nlem;
+ int imp = impmmmp + mmp;
+ int impe = imp + cil->nlem;
+ double cgc = cg1(lmpml, mu, l80, m);
+ summ = dconjg(ac[i - 1]) * ac[imp - 1];
+ sume = dconjg(ac[i - 1]) * ac[impe - 1];
+ suem = dconjg(ac[ie - 1]) * ac[imp - 1];
+ suee = dconjg(ac[ie - 1]) * ac[impe - 1];
+ if (lmpml != 0) {
+ summ *= uimmp;
+ sume *= uimmp;
+ suem *= uimmp;
+ suee *= uimmp;
+ }
+ // label 65
+ gap[imu90 - 1] += (cgc * (
+ summ * zpv[l80 - 1][ilmp80 - 1][0][0] +
+ sume * zpv[l80 - 1][ilmp80 - 1][0][1] +
+ suem * zpv[l80 - 1][ilmp80 - 1][1][0] +
+ suee * zpv[l80 - 1][ilmp80 - 1][1][1]
+ )
+ );
+ }
+ } // im70 loop
+ } // ilmp80 loop
+ } // l80 loop
+ } // imu90 loop
+ sume = gap[0] * ccr->cimu;
+ suee = gap[1] * ccr->cof;
+ suem = gap[2] * ccr->cimu;
+ fffs[0] = (sume - suem).real();
+ fffs[1] = ((sume + suem) * uim).real();
+ fffs[2] = suee.real();
+ delete[] gap;
+}
+
+/*! C++ porting of SAMP
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param tmsm: Vector of complex. QUESTION: definition?
+ * \param tmse: Vector of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ *
+ * This function builds the AC vector using TMSM, TMSE and WS.
+ */
+void samp(
+ std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
+ std::complex<double> *ws, CIL *cil
+) {
+ 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;
+ for (int im20 = 0; im20 < ltpo; im20++) {
+ int ie = i + cil->nlem;
+ ac[i] = tmsm[l20] * ws[i];
+ ac[ie] = tmse[l20] * ws[ie];
+ i++;
+ } // im20 loop
+ } // l20 loop
+}
+
/*! C++ porting of SAMPOA
*
* \param ac: Vector of complex. QUESTION: definition?
diff --git a/src/trapping/lffft.cpp b/src/trapping/lffft.cpp
index 33d516a5ae5929d5154f4d0821072751986f31de..4b8fb7ab94c63eb0d662233c78bdfb428c0a2018 100644
--- a/src/trapping/lffft.cpp
+++ b/src/trapping/lffft.cpp
@@ -268,12 +268,12 @@ int main() {
if (is != 2222) {
if (is != 1111) {
if (is > 0) { // Goes to 305
- // Would call CAMP(AC,AM0M,WS)
+ camp(ac, am0m, ws, cil);
} else if (is < 0) { // Goes to 405
- // Would call CZAMP(AC,AMD,INDAM,WS)
+ czamp(ac, amd, indam, ws, cil);
}
} else {
- // Would call SAMP(AC,TMSM,TMSE,WS)
+ samp(ac, tmsm, tmse, ws, cil);
}
} else {
sampoa(ac, tms, ws, cil);
@@ -282,7 +282,7 @@ int main() {
if (jft <= 0) {
double *fffe = new double[3]();
double *fffs = new double[3]();
- // Would call FFRF(ZPV,AC,WS,FFFE,FFFS)
+ ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
if (jss == 1) {
// Writes to 66
} else { // label 450