Finish migration of trapping to C++

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];

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