Replicate edfb.f as edfb.cpp

src/sphere/edfb.cpp

@@ -17,12 +17,16 @@ string *load_file(string file_name, int *count);
 
 int main(int argc, char **argv) {
   // Common variables set
-  complex<double> *dc0, *dc0m;
+  complex<double> *dc0, ***dc0m;
   double *ros, **rcf;
   int *iog, *nshl;
   double *xiv, *wns, *wls, *pus, *evs, *vss;
   string vns[5];
 
+  /// A helper variable to set the size of dc0m
+  int max_nsh = 0;
+  int ici;
+
   // Input file reading section
   int num_lines = 0;
   int last_read_line; //!< Keep track of where INXI left the input stream
@@ -35,6 +39,7 @@ int main(int argc, char **argv) {
   double exdc, wp, xip;
   int exdc_exp, wp_exp, xip_exp;
   int idfc, nxi, instpc, insn;
+  int nsh;
   sscanf(
     file_lines[1].c_str(),
     " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d",
@@ -102,8 +107,9 @@ int main(int argc, char **argv) {
     sscanf(file_lines[++last_read_line].c_str(), " %d %9lf D%d", &i_val, &ros_val, &ros_val_exp);
     nshl[i113 - 1] = i_val;
     ros[i113 - 1] = ros_val * pow(10.0, ros_val_exp);
-    int nsh = nshl[i113 -1];
+    nsh = nshl[i113 -1];
     if (i113 == 1) nsh += ies;
+    if ((nsh + 1) / 2 + ies > max_nsh) max_nsh = (nsh + 1) / 2 + ies;
     rcf[i113 - 1] = new double[nsh];
     for (int ns = 0; ns < nsh; ns++) {
       double ns_rcf;
@@ -112,14 +118,98 @@ int main(int argc, char **argv) {
       rcf[i113 -1][ns] = ns_rcf * pow(10.0, ns_rcf_exp);
     }
   }
-  if (idfc < 0) {
+  // The FORTRAN code writes an auxiliary file in binary format. This should
+  // be avoided or possibly replaced with the use of standard file formats for
+  // scientific use (e.g. FITS).
+  ofstream tedf_file;
+  tedf_file.open("c_TEDF", ofstream::binary);
+  tedf_file.write(reinterpret_cast<char *>(&nsph), sizeof(nsph));
+  for (int iogi = 0; iogi < nsph; iogi++)
+    tedf_file.write(reinterpret_cast<char *>(iog + iogi), sizeof(iog[iogi]));
+  tedf_file.write(reinterpret_cast<char *>(&exdc), sizeof(exdc));
+  tedf_file.write(reinterpret_cast<char *>(&wp), sizeof(wp));
+  tedf_file.write(reinterpret_cast<char *>(&xip), sizeof(xip));
+  tedf_file.write(reinterpret_cast<char *>(&idfc), sizeof(idfc));
+  tedf_file.write(reinterpret_cast<char *>(&nxi), sizeof(nxi));
+  for (int i115 = 1; i115 <= nsph; i115++) {
+    if (iog[i115 - 1] < i115) continue;
+    tedf_file.write(reinterpret_cast<char *>(nshl + i115 - 1), sizeof(nshl[i115 - 1]));
+    tedf_file.write(reinterpret_cast<char *>(ros + i115 - 1), sizeof(ros[i115 - 1]));
+    nsh = nshl[i115 - 1];
+    if (i115 == 1) nsh += ies;
+    for (int ins = 0; ins < nsh; ins++)
+      tedf_file.write(reinterpret_cast<char *>(rcf[i115 - 1] + ins), sizeof(rcf[i115 - 1][ins]));
+  }
+  // Remake the dc0m matrix.
+  dc0m = new complex<double>**[max_nsh];
+  for (int dim1 = 0; dim1 < max_nsh; dim1++) {
+    dc0m[dim1] = new complex<double>*[nsph];
+    for (int dim2 = 0; dim2 < nxi; dim2++) {
+      dc0m[dim1][dim2] = new complex<double>[nxi];
+    }
+  }
+  for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
+    if (idfc != 0 && jxi468 > 1) continue;
+    for (int i162 = 1; i162 <= nsph; i162++) {
+      if (iog[i162 - 1] < i162) continue;
+      nsh = nshl[i162 - 1];
+      ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+      if (i162 == 1) ici = ici + ies;
+      for (int i157 = 0; i157 < ici; i157++) {
+	double dc0_real, dc0_img;
+	int dc0_real_exp, dc0_img_exp;
+	sscanf(file_lines[++last_read_line].c_str(), " (%8lf D%d, %8lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
+	dc0_real *= pow(10.0, dc0_real_exp);
+	dc0_img *= pow(10.0, dc0_img_exp);
+	dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
+	// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
+	// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
+	// real part and 8 bytes to represent the imaginary one.
+	tedf_file.write(reinterpret_cast<char *>(&dc0_real), sizeof(dc0_real));
+	tedf_file.write(reinterpret_cast<char *>(&dc0_img), sizeof(dc0_img));
+      }
+    }
+  }
+  tedf_file.close();
+  if (idfc != 0) {
     fprintf(output, "  DIELECTRIC CONSTANTS\n");
+    for (int i473 = 1; i473 <= nsph; i473++) {
+      if (iog[i473 - 1] != i473) continue;
+      ici = (nshl[i473 - 1] + 1) / 2;
+      if (i473 == 1) ici += ies;
+      fprintf(output, " SPHERE N. %d\n", i473);
+      for (int ic472 = 0; ic472 < ici; ic472++) {
+	double dc0_real = dc0m[ic472][i473 - 1][0].real(), dc0_img = dc0m[ic472][i473 - 1][0].imag();
+	fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
+      }
+    }
+  } else {
+    fprintf(output, "  DIELECTRIC FUNCTIONS\n");
   }
   fclose(output);
   return 0;
 }
 
-string *load_file(string file_name, int *count) {
+/*! \fn load_file(string, int*)
+ * \brief Load a text file as a sequence of strings in memory.
+ *
+ * The configuration of the field expansion code in FORTRAN uses
+ * shared memory access and file I/O operations managed by different
+ * functions. Although this approach could be theoretically replicated,
+ * it is more convenient to handle input and output to distinct files
+ * using specific functions. load_file() helps in the task of handling
+ * input such as configuration files or text data structures that need
+ * to be loaded entirely. The function performs a line-byline scan of
+ * the input file and returns an array of strings that can be later
+ * parsed and ingested by the concerned code blocks. An optional pointer
+ * to integer allows the function to keep track of the number of file
+ * lines that were read, if needed.
+ *
+ * \param string file_name: The path of the file to be read.
+ * \param [int *count = NULL]: Pointer to an integer recording the number of lines.
+ * \return string*: An array of strings, one for each input file line.
+ */
+string *load_file(string file_name, int *count = 0) {
   fstream input_file(file_name.c_str(), ios::in);
   List<string> file_lines = List<string>();
   string line;
@@ -132,6 +222,6 @@ string *load_file(string file_name, int *count) {
     input_file.close();
   }
   string *array_lines = file_lines.to_array();
-  *count = file_lines.length();
+  if (count != 0) *count = file_lines.length();
   return array_lines;
 }