Implement LFFFT formatted output

src/trapping/clffft.cpp

@@ -251,6 +251,12 @@ void lffft(string data_file, string output_path) {
 	  const int nlmmt = nlmm + nlmm;
 	  ws = new complex<double>[nlmmt]();
 	  if (lm > le) wsl = new complex<double>[nlmmt]();
+	  // FORTRAN writes two output formatted files without opening them
+	  // explicitly. It is assumed thay can be opened here.
+	  string out66_name = output_path + "/c_out66.txt";
+	  string out67_name = output_path + "/c_out67.txt";
+	  FILE *output66 = fopen(out66_name.c_str(), "w");
+	  FILE *output67 = fopen(out67_name.c_str(), "w");
 	  for (int iz475 = 0; iz475 < nzv; iz475++) {
 	    for (int iy475 = 0; iy475 < nyv; iy475++) {
 	      for (int ix475 = 0; ix475 < nxv; ix475++) {
@@ -299,6 +305,18 @@ void lffft(string data_file, string output_path) {
 		    ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
 		    if (jss == 1) {
 		      // Writes to 66
+		      fprintf(
+			      output66, " %18.16lE%18.16lE%18.16lE\n",
+			      fffe[0], fffs[0], fffe[0] - fffs[0]
+			      );
+		      fprintf(
+			      output66, " %18.16lE%18.16lE%18.16lE\n",
+			      fffe[1], fffs[1], fffe[1] - fffs[1]
+			      );
+		      fprintf(
+			      output66, " %18.16lE%18.16lE%18.16lE\n",
+			      fffe[2], fffs[2], fffe[2] - fffs[2]
+			      );
 		    } else { // label 450
 		      for (int i = 0; i < 3; i++) {
 			double value = fffe[i] - fffs[i];
@@ -306,6 +324,11 @@ void lffft(string data_file, string output_path) {
 		      }
 		      if (jtw == 1) {
 			// Writes to 66
+			fprintf(
+				output66, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
+				ix475 + 1, iy475 + 1, iz475 + 1,
+				fffe[0] - fffs[0], fffe[1] - fffs[1], fffe[2] - fffs[2]
+				);
 		      }
 		    }
 		    if (jft < 0) goto475 = true;
@@ -319,6 +342,18 @@ void lffft(string data_file, string output_path) {
 		    ffrt(ac, ws, ffte, ffts, cil);
 		    if (jss == 1) {
 		      // Writes to 67
+		      fprintf(
+			      output67, " %18.16lE%18.16lE%18.16lE\n",
+			      ffte[0], ffts[0], ffte[0] - ffts[0]
+			      );
+		      fprintf(
+			      output67, " %18.16lE%18.16lE%18.16lE\n",
+			      ffte[1], ffts[1], ffte[1] - ffts[1]
+			      );
+		      fprintf(
+			      output67, " %18.16lE%18.16lE%18.16lE\n",
+			      ffte[2], ffts[2], ffte[2] - ffts[2]
+			      );
 		    } else { // label 470
 		      for (int i = 0; i < 3; i++) {
 			double value = ffte[i] - ffts[i];
@@ -326,6 +361,11 @@ void lffft(string data_file, string output_path) {
 		      }
 		      if (jtw == 1) {
 			// Writes to 67
+			fprintf(
+				output67, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
+				ix475 + 1, iy475 + 1, iz475 + 1,
+				ffte[0] - ffts[0], ffte[1] - ffts[1], ffte[2] - ffts[2]
+				);
 		      }
 		    }
 		    delete[] ffte;
@@ -339,6 +379,8 @@ void lffft(string data_file, string output_path) {
 	    if (jft <= 0) tlfff.close();
 	    if (jft >= 0) tlfft.close();
 	  }
+	  fclose(output66);
+	  fclose(output67);
 	}
       }
       binary_input.close();