Complete migration of FRFME to C++

src/include/tra_subs.h

@@ -7,12 +7,58 @@
 #endif
 
 //Externally defined subroutines
+extern std::complex<double> dconjg(std::complex<double> z);
 extern void sphar(
 	   double cosrth, double sinrth, double cosrph, double sinrph,
 	   int ll, std::complex<double> *ylm
 );
 //End of externally defined subroutines
 
+/*! C++ porting of PWMALP
+ *
+ * \param w: Matrix of complex. QUESTION: definition?
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ylm: Vector of complex
+ * \param lw: `int`
+ */
+void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
+  std::complex<double> cp1, cm1, cp2, cm2, cl;
+  const std::complex<double> uim(0.0, 1.0);
+  const double four_pi = acos(0.0) * 8.0;
+  const int nlwm = lw * (lw + 2);
+  cm1 = 0.5 * std::complex<double>(up[0], up[1]);
+  cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
+  double cz1 = up[2];
+  cm2 = 0.5 * std::complex<double>(un[0], un[1]);
+  cp2 = 0.5 * std::complex<double>(un[0], -un[1]);
+  double cz2 =un[2];
+  for (int l20 = 1; l20 <= lw; l20++) {
+    int lf = l20 + 1;
+    int lftl = lf * l20;
+    double x = 1.0 * lftl;
+    std::complex<double> cl = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);
+    int mv = l20 + lf;
+    int m = -lf;
+    for (int mf20 = 1; mf20 <= mv; mf20++) {
+      m++;
+      int k = lftl + m;
+      x = 1.0 * (lftl - m * (m + 1));
+      double cp = sqrt(x);
+      x = 1.0 * (lftl - m * (m - 1));
+      double cm = sqrt(x);
+      double cz = 1.0 * m;
+      w[k - 1][0] = dconjg(cp1 * cp * ylm[k + 1] + cm1 * cm * ylm[k - 1] + cz1 * cz * ylm[k]) * cl;
+      w[k - 1][1] = dconjg(cp2 * cp * ylm[k + 1] + cm2 * cm * ylm[k - 1] + cz2 * cz * ylm[k]) * cl;
+    } // mf20 loop
+  } // l20 loop
+  for (int k30 = 0; k30 < nlwm; k30++) {
+    int i = k30 + nlwm;
+    w[i][0] = uim * w[k30][1];
+    w[i][1] = -uim * w[k30][0];
+  } // k30 loop
+}
+
 /*! C++ porting of WAMFF
  *
  * \param wk: Vector of complex. QUESTION: definition?
@@ -103,6 +149,7 @@ void wamff(
     un[1] = cph;
     un[2] = 0.0;
     // Would call PWMALP(W,UP,UN,YLM,LM)
+    pwmalp(w, up, un, ylm, lm);
     double apfa = sth * apfafa;
     if (spd > 0.0) {
       double sthl = sth * rir;
@@ -206,7 +253,7 @@ void frfmer(
       double vkn = sqrt(sqn);
       if (vkn <= vknmx) {
 	double vkz = sqrt(sq - sqn);
-	// Would call WAMFF(WK,VKX,VKY,VKZ,LE,APFAFA,TRA,SPD,RIR,FTCN,LMODE,PMF)
+	wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
 	for (int j = 0; j < nlemt; j++) {
 	  double vreal = wk[j].real();
 	  double vimag = wk[j].imag();

src/trapping/frfme.cpp

@@ -24,12 +24,15 @@ using namespace std;
 int main() {
   string data_file = "../../test_data/trapping/DFRFME";
   string output_path = ".";
+  string tfrfme_name = output_path + "/c_TFRFME";
+  fstream tfrfme;
   char namef[7];
   char more;
   double *xv = NULL, *yv = NULL, *zv = NULL;
   double *vkv = NULL, **vkzm = NULL;
   complex<double> *wk = NULL, **w = NULL, **wsum = NULL;
   const complex<double> cc0(0.0, 0.0);
+  const complex<double> uim(0.0, 1.0);
   int line_count = 0, last_read_line = 0;
   regex re = regex("-?[0-9]+");
   string *file_lines = load_file(data_file, &line_count);
@@ -53,8 +56,7 @@ int main() {
   // End of vector size variables
   if (jlmf != 1) {
     int nxv, nyv, nzv;
-    fstream tfrfme;
-    tfrfme.open("c_TFRFME", ios::in | ios::binary);
+    tfrfme.open(tfrfme_name, ios::in | ios::binary);
     if (tfrfme.is_open()) {
       tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
       tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
@@ -79,7 +81,7 @@ int main() {
       fstream temptape2;
       temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);
       if (temptape2.is_open()) {
-	vkv = new double[nkv]();
+	//vkv = new double[nkv]();
 	for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
 	vkzm = new double*[nkv];
 	vkzm_size = nkv;
@@ -106,9 +108,6 @@ int main() {
       } else {
 	printf("ERROR: could not open TEMPTAPE2 file.\n");
       }
-      wsum = new complex<double>*[nlmmt];
-      wsum_size = nlmmt;
-      for (int wsi = 0; wsi < jlmf - 1; wsi++) wsum[wsi] = new complex<double>[nrvc]();
       for (int ixyz12 = 0; ixyz12 < nrvc; ixyz12++) {
 	for (int j12 = 0; j12 < jlmf - 1; j12++) {
 	  double vreal, vimag;
@@ -220,8 +219,6 @@ int main() {
 	// label 22
 	nlmmt = lm * (lm + 2) * 2;
 	nks = nksh * 2;
-	//wsum = new complex<double>*[nlmmt];
-	wsum_size = nlmmt;
 	nkv = nks + 1;
 	double vkm = vk * exri;
 	vknmx = vk * an;
@@ -259,14 +256,15 @@ int main() {
 	} // i27 loop
 	int nrvc = nxv * nyv * nzv;
 	int nkshpo = nksh + 1;
+	wsum = new complex<double>*[nlmmt];
+	wsum_size = nlmmt;
+	for (int wsi = 0; wsi < nlmmt; wsi++) wsum[wsi] = new complex<double>[nrvc]();
 	vkv = new double[nkv]();
 	// vkv[nksh] = 0.0;
 	for (int i28 = nkshpo; i28 <= nks; i28++) {
 	  vkv[i28] = vkv[i28 - 1] + delk;
 	  vkv[nkv - i28 - 1] = -vkv[i28];
 	} // i28 loop
-	fstream tfrfme;
-	string tfrfme_name = output_path + "/c_TFRFME";
 	tfrfme.open(tfrfme_name.c_str(), ios::out | ios::binary);
 	if (tfrfme.is_open()) {
 	  tfrfme.write(reinterpret_cast<char *>(&lmode), sizeof(int));
@@ -297,7 +295,6 @@ int main() {
 	  for (int jx = 0; jx < nkv; jx++)
 	    temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
 	  frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, temptape1, temptape2);
-	  temptape1.write(reinterpret_cast<char *>(&apfafa), sizeof(double)); // Place-holder
 	  temptape1.close();
 	  temptape2.write(reinterpret_cast<char *>(&apfafa), sizeof(double));
 	  temptape2.write(reinterpret_cast<char *>(&pmf), sizeof(double));
@@ -313,7 +310,103 @@ int main() {
 	  temptape2.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
 	  temptape2.write(reinterpret_cast<char *>(&nrvc), sizeof(int));
 	  temptape2.close();
+	  temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);
+	  vkv = new double[nkv]();
+	  for (int jx = 0; jx < nkv; jx++)
+	    temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
+	  vkzm = new double*[nkv];
+	  vkzm_size = nkv;
+	  for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();
+	  for (int jy40 = 0; jy40 < nkv; jy40++) {
+	    for (int jx40 = 0; jx40 < nkv; jx40++)
+	      temptape2.read(reinterpret_cast<char *>(&(vkzm[jx40][jy40])), sizeof(double));
+	  } // jy40 loop
+	  temptape2.close();
+	  wk = new complex<double>[nlmmt];
+	  w = new complex<double>*[nkv];
+	  for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
+	  for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+	    temptape1.open("c_TEMPTAPE1", ios::in | ios::binary);
+	    for (int jy50 = 0; jy50 < nkv; jy50++) {
+	      for (int jx50 = 0; jx50 < nkv; jx50++) {
+		for (int i = 0; i < nlmmt; i++) {
+		  double vreal, vimag;
+		  temptape1.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+		  temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+		  wk[i] = complex<double>(vreal, vimag);
+		}
+		w[jx50][jy50] = wk[j80];
+	      } // jx50
+	    } // jy50 loop
+	    temptape1.close();
+	    int ixyz = 0;
+	    for (int iz75 = 0; iz75 < nzv; iz75++) {
+	      double z = zv[iz75]  + frsh;
+	      for (int iy70 = 0; iy70 < nyv; iy70++) {
+		double y = yv[iy70];
+		for (int ix65 = 0; ix65 < nxv; ix65++) {
+		  double x = xv[ix65];
+		  ixyz++;
+		  complex<double> 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));
+		    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);
+		    for (int jx55 = 1; jx55 < nks; jx55++) {
+		      vkx = vkv[jx55];
+		      double vkz = vkzm[jx55][jy60];
+		      complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
+		      sumx += (w[jx55][jy60] * phas);
+		    } // jx55 loop
+		    if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
+		    sumy += sumx;
+		  } // jy60 loop
+		  wsum[j80][ixyz - 1] = sumy * delks;
+		} // ix65 loop
+	      } // iy70 loop
+	    } // iz75 loop
+	  } // j80 loop
+	  if (jlmf != 1) {
+	    tfrfme.open(tfrfme_name, ios::in | ios::out | ios::binary);
+	    tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+	    tfrfme.read(reinterpret_cast<char *>(&vk), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&exri), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&an), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&ff), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&tra), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+	    tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
+	    for (int i = 0; i < nxv; i++) tfrfme.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
+	    for (int i = 0; i < nyv; i++) tfrfme.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
+	    for (int i = 0; i < nzv; i++) tfrfme.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+	  }
+	  // label 88
+	  for (int ixyz = 0; ixyz < nrvc; ixyz++) {
+	    for (int j = 0; j< jlml; j++) {
+	      double vreal = wsum[j][ixyz].real();
+	      double vimag = wsum[j][ixyz].imag();
+	      tfrfme.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+	      tfrfme.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+	    } // j loop
+	  } // ixyz loop
 	  tfrfme.close();
+	  string output_name = output_path + "/c_OFRFME";
+	  FILE *output = fopen(output_name.c_str(), "w");
+	  fprintf(output, " IF JLML < NLMMT, PRESERVE TEMPTAPE1, TEMPTAPE2, AND TFRFRME,\n");
+	  fprintf(output, " AND RESTART LM RUN WITH JLMF = JLML+1\n");
+	  if (spd > 0.0) fprintf(output, "  FRSHMX =%15.7lE\n", fshmx);
+	  fprintf(output, "  FRSH =%15.7lE\n", frsh);
+	  fclose(output);
 	} else { // Should never happen.
 	  printf("ERROR: could not open TFRFME file for output.\n");
 	}
@@ -328,7 +421,7 @@ int main() {
     }
   }
   // label 45
-  
+  if (tfrfme.is_open()) tfrfme.close();
   delete[] file_lines;
   if (xv != NULL) delete[] xv;
   if (yv != NULL) delete[] yv;
@@ -342,6 +435,11 @@ int main() {
     for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
     delete[] wsum;
   }
+  if (wk != NULL) delete[] wk;
+  if (w != NULL) {
+    for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
+    delete[] w;
+  }
   printf("Done.\n");
   return 0;
 }