Migrate FFRF subroutine to C++

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?

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