Start porting edfb.f to edfb.cpp

src/sphere/edfb.cpp

@@ -5,12 +5,133 @@
 
 #include <cstdio>
 #include <cmath>
-#include <string>
+#include <complex>
+#include <cstring>
+#include <iostream>
+#include <fstream>
 #include "List.h"
-#include "edfb.h"
+
+using namespace std;
+
+string *load_file(string file_name, int *count);
 
 int main(int argc, char **argv) {
-  DEDFB input_data = DEDFB(std::string("../../test_data/sphere/DEDFB"));
-  input_data.print();
+  // Common variables set
+  complex<double> *dc0, *dc0m;
+  double *ros, **rcf;
+  int *iog, *nshl;
+  double *xiv, *wns, *wls, *pus, *evs, *vss;
+  string vns[5];
+
+  // Input file reading section
+  int num_lines = 0;
+  int last_read_line; //!< Keep track of where INXI left the input stream
+  string *file_lines = load_file("../../test_data/sphere/DEDFB", &num_lines);
+
+  // Configuration code
+  int nsph, ies;
+  sscanf(file_lines[0].c_str(), " %d %d", &nsph, &ies);
+  if (ies != 0) ies = 1;
+  double exdc, wp, xip;
+  int exdc_exp, wp_exp, xip_exp;
+  int idfc, nxi, instpc, insn;
+  sscanf(
+    file_lines[1].c_str(),
+    " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d",
+    &exdc, &exdc_exp,
+    &wp, &wp_exp,
+    &xip, &xip_exp,
+    &idfc, &nxi, &instpc, &insn
+  );
+  exdc *= pow(10.0, exdc_exp);
+  wp *= pow(10.0, wp_exp);
+  xip *= pow(10.0, xip_exp);
+
+  FILE *output = fopen("c_OEDFB", "w");
+  // FORTRAN starts subroutine INXI at this point
+  const double pigt = acos(0.0) * 4.0;
+  const double evc = 6.5821188e-16;
+  if (idfc >= 0) {
+    printf("Not walked by input data.\n");
+  } else {
+    if (instpc < 1) {
+      // In this case the XI vector is explicitly defined.
+      // Test input comes this way.
+      vns[insn] = "XIV";
+      List<double> xi_vector;
+      double xi;
+      int xi_exp;
+      sscanf(file_lines[2].c_str(), " %9lE D%d", &xi, &xi_exp);
+      xi *= pow(10.0, xi_exp);
+      xi_vector.set(0, xi);
+      for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
+	sscanf(file_lines[2 + jxi310].c_str(), " %9lE D%d", &xi, &xi_exp);
+	xi *= pow(10.0, xi_exp);
+	xi_vector.append(xi);
+	last_read_line = 2 + jxi310;
+      }
+      vss = xi_vector.to_array();
+      xiv = xi_vector.to_array();
+      double pu = xip + wp;
+      double wn = pu / 3.0e8;
+      fprintf(output, "          XIP          WN           WL           PU           EV\n");
+      fprintf(output, "     %13.4lE", xip);
+      fprintf(output, "%13.4lE", wn);
+      fprintf(output, "%13.4lE", pigt / wn);
+      fprintf(output, "%13.4lE", pu);
+      fprintf(output, "%13.4lE\n", pu * evc);
+      fprintf(output, "  SCALE FACTORS XI\n", pu * evc);
+      for (int jxi6612 = 1; jxi6612 <= nxi; jxi6612++)
+	fprintf(output, "%5d%13.4lE\n", jxi6612, xiv[jxi6612 - 1]);
+      //INXI branch ends here.
+    }
+  }
+  last_read_line++;
+  iog = new int[nsph];
+  for (int i = 0; i < nsph; i++) {
+    sscanf(file_lines[last_read_line].c_str(), " %d", (iog + i));
+  }
+  nshl = new int[nsph];
+  ros = new double[nsph];
+  rcf = new double*[nsph];
+  for (int i113 = 1; i113 <= nsph; i113++) {
+    int i_val;
+    double ros_val;
+    int ros_val_exp;
+    if (iog[i113 - 1] < i113) continue;
+    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];
+    if (i113 == 1) nsh += ies;
+    rcf[i113 - 1] = new double[nsh];
+    for (int ns = 0; ns < nsh; ns++) {
+      double ns_rcf;
+      int ns_rcf_exp;
+      sscanf(file_lines[++last_read_line].c_str(), " %8lf D%d", &ns_rcf, &ns_rcf_exp);
+      rcf[i113 -1][ns] = ns_rcf * pow(10.0, ns_rcf_exp);
+    }
+  }
+  if (idfc < 0) {
+    fprintf(output, "  DIELECTRIC CONSTANTS\n");
+  }
+  fclose(output);
   return 0;
 }
+
+string *load_file(string file_name, int *count) {
+  fstream input_file(file_name.c_str(), ios::in);
+  List<string> file_lines = List<string>();
+  string line;
+  if (input_file.is_open()) {
+    getline(input_file, line);
+    file_lines.set(0, line);
+    while (getline(input_file, line)) {
+      file_lines.append(line);
+    }
+    input_file.close();
+  }
+  string *array_lines = file_lines.to_array();
+  *count = file_lines.length();
+  return array_lines;
+}

src/sphere/edfb.h

-/**
- * \brief A class to represent EDFB configuration data.
- *
- * This class replicates the structure of formatted DEDFB input data
- * files for the EDFB program. The class is built on the SPHERE/DEDFB
- * file template, but it is meant to grant compliance with the CLUSTER
- * case too.
- *
- * Compatibility between the SPHERE and the CLUSTER cases is planned
- * through the introduction of dynamic pointers that can fit both cases,
- * provided that the DEDFB formatting has been properly understood. In
- * case not, the suggested work-around is to develop specialized classes
- * for the SPHERE case and the CLUSTER case. In that case, the current
- * class would already represent a valid implementation for the single
- * sphere problem.
- *
- * It should be noted that the variable names have been left intentionally
- * equal to the FORTRAN equivalents. In addition they are not documented
- * yet. Documentation and naming may be revised after the code has been
- * tested to succesfully reproduce the original work-flow.
- */
-class DEDFB {
- protected:
-  int nsph, ies;
-  double exdc, wp, xip;
-  int idfc, nxi, instpc, insn;
-  double *xi_vector;
-  int *iog;
-  int *nshl;
-  double *ros;
-  double *rcf;
-  double *dc0;
-
-public:
-  /*! \fn DEDFB(std::string)
-   * \brief Data structure constructor.
-   *
-   * The default procedure to import a configuration file is to build a
-   * configuration object, according to one of the following flavours:
-   *
-   * DEDFB edfb_cfg = DEDFB("FILE_NAME"); // configuration object
-   *
-   * or
-   *
-   * DEDFB *edfb_cfg_ptr = new DEDFB("FILE_NAME"); // pointer version
-   */
-  DEDFB(std::string file_name);
-  
-  /*! \fn print()
-   * \brief Print the data structure to console.
-   *
-   * This is a simple function to print the data read from a configuration
-   * file to the screen. This operation can be useful for debug and testing.
-   */
-  void print();
-};
-
-DEDFB::DEDFB(std::string file_name) {
-  FILE *input_file = fopen(file_name.c_str(), "r");
-  int exp_exdc, exp_wp, exp_xip;
-  List<double> xi_list(1);
-  fscanf(input_file, " %d %d", &nsph, &ies);
-  fscanf(
-	 input_file,
-	 " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d",
-	 &exdc, &exp_exdc,
-	 &wp, &exp_wp,
-	 &xip, &exp_xip,
-	 &idfc, &nxi, &instpc, &insn
-	 );
-  exdc *= pow(10.0, exp_exdc);
-  wp *= pow(10.0, exp_wp);
-  xip *= pow(10.0, exp_xip);
-  for (int i_xi = 0; i_xi < nxi; i_xi++) {
-    double xi;
-    int exp_xi;
-    fscanf(input_file, " %9lf D%d", &xi, &exp_xi);
-    xi *= pow(10.0, exp_xi);
-    if (i_xi == 0) {
-      xi_list.set(0, xi);
-    } else {
-      xi_list.append(xi);
-    }
-  }
-  xi_vector = xi_list.to_array();
-  iog = new int[nsph];
-  for (int i_iog = 0; i_iog < nsph; i_iog++) {
-    fscanf(input_file, " %d", (iog + i_iog));
-  }
-  nshl = new int[nsph];
-  ros = new double[nsph];
-  for (int i_nshl = 0; i_nshl < nsph; i_nshl++) {
-    double i_ros;
-    int i_ros_exp;
-    fscanf(input_file, " %d %9lf D%d", (nshl + i_nshl), &i_ros, &i_ros_exp);
-    i_ros *= pow(10.0, i_ros_exp);
-    ros[i_nshl] = i_ros;
-  }
-  rcf = new double[nsph * nshl[0]];
-  for (int i_rcf = 0; i_rcf < nshl[0]; i_rcf++) {
-    for (int j_rcf = 0; j_rcf < nsph; j_rcf++) {
-      double ij_rcf;
-      int ij_rcf_exp;
-      fscanf(input_file, " %8lf D%d", &ij_rcf, &ij_rcf_exp);
-      ij_rcf *= pow(10.0, ij_rcf_exp);
-      rcf[(nsph * i_rcf) + j_rcf] = ij_rcf;
-    }
-  }
-  dc0 = new double[nsph * nshl[0]];
-  for (int i_dc0 = 0; i_dc0 < nsph; i_dc0++) {
-    double ij_dc0;
-    int ij_dc0_exp;
-    fscanf(input_file, " (");
-    for (int j_dc0 = 0; j_dc0 < nshl[0] - 1; j_dc0++) {
-      fscanf(input_file, " %8lf D%d,", &ij_dc0, &ij_dc0_exp);
-      ij_dc0 *= pow(10.0, ij_dc0_exp);
-      dc0[(nshl[0] * i_dc0) + j_dc0] = ij_dc0;
-    }
-    fscanf(input_file, " %8lf D%d)", &ij_dc0, &ij_dc0_exp);
-    ij_dc0 *= pow(10.0, ij_dc0_exp);
-    dc0[(nshl[0] * i_dc0) + nshl[0] - 1] = ij_dc0;
-  }
-  fclose(input_file);
-}
-
-void DEDFB::print() {
-  printf("### CONFIGURATION DATA ###\n");
-  printf("NSPH   = %d\n", nsph);
-  printf("IES    = %d\n", ies);
-  printf("EXDC   = %E\n", exdc);
-  printf("WP     = %E\n", wp);
-  printf("XIP    = %E\n", xip);
-  printf("IDFC   = %d\n", idfc);
-  printf("NXI    = %d\n", nxi);
-  printf("INSTPC = %d\n", instpc);
-  printf("INSN   = %d\n", insn);
-  printf("XIV    = [ %lE", xi_vector[0]);
-  for (int i_xiv = 1; i_xiv < nxi; i_xiv++) printf(", %lE", xi_vector[i_xiv]);
-  printf(" ]\n");
-  printf("IOG    = [ %d", iog[0]);
-  for (int i_iog = 1; i_iog < nsph; i_iog++) printf(", %d", iog[i_iog]);
-  printf(" ]\n");
-  printf("NSHL   = [ %d", nshl[0]);
-  for (int i_nshl = 1; i_nshl < nsph; i_nshl++) printf(", %d", nshl[i_nshl]);
-  printf(" ]\n");
-  printf("ROS    = [ %lE", ros[0]);
-  for (int i_ros = 1; i_ros < nsph; i_ros++) printf(", %lE", ros[i_ros]);
-  printf(" ]\n");
-  printf("RCF    = [ %lE", rcf[0]);
-  for (int i_rcf = 1; i_rcf < nsph * nshl[0]; i_rcf++) printf(", %lE", rcf[i_rcf]);
-  printf(" ]\n");
-  printf("DC0    = [ %lE", dc0[0]);
-  for (int i_dc0 = 1; i_dc0 < nsph * nshl[0]; i_dc0++) printf(", %lE", dc0[i_dc0]);
-  printf(" ]\n");
-  printf("######   END DATA   ######\n");
-  
-}
-
-void inxi(
-  DEDFB *data
-); //!< \brief Initialization process