Merge branch 'c_fortran_binary_io' into 'master'

Wrap FORTRAN binary I/O in C++

See merge request giacomo.mulas/np_tmcode!4

.gitignore

+build/include
 build/cluster/*
 build/sphere/*
 build/trapping/*

doc/src/config.dox

@@ -342,7 +342,7 @@ OPTIMIZE_OUTPUT_SLICE  = NO
 #
 # Note see also the list of default file extension mappings.
 
-EXTENSION_MAPPING      =
+EXTENSION_MAPPING      = f=FortranFixed f90=FortranFree
 
 # If the MARKDOWN_SUPPORT tag is enabled then doxygen pre-processes all comments
 # according to the Markdown format, which allows for more readable
@@ -984,6 +984,7 @@ INPUT_FILE_ENCODING    =
 # *.f18, *.f, *.for, *.vhd, *.vhdl, *.ucf, *.qsf and *.ice.
 
 FILE_PATTERNS          = *.cpp \
+		         *.f \
 		         *.h \
 			 *.md
 

src/sphere/List.h → src/include/List.h

@@ -154,8 +154,8 @@ template<class T> class List {
   T* to_array() {
     T *array = new T[size];
     current = last;
-    for (int i = size; i > 0; i--) {
-      array[i - 1] = current->value;
+    for (int i = size - 1; i > -1; i--) {
+      array[i] = current->value;
       current = current->p_prev;
     }
     return array;

src/include/file_io.h

+/*! \file file_io.h
+ *
+ * C++ wrapper of FORTRAN I/O operations with files.
+ */
+
+/*! \brief Open a file for subsequent access.
+ *
+ * \param uid: `int*` Pointer to the unit ID to be associated to the file.
+ * \param name: `char*` C-string for file name (max. length of 63).
+ * \param sta: `char*` C-string for file status (max. length of 7).
+ * \param mode: `char*` C-string for access mode (max. length of 11).
+ */
+extern "C" void open_file_(int* uid, const char* name, const char* sta, const char* mode);
+/*! \brief Close a previously opened file.
+ *
+ * \param uid: `int*` Pointer to the unit ID of the file.
+ */
+extern "C" void close_file_(int* uid);
+/*! \brief Read an integer value from a file.
+ *
+ * \param uid: `int*` Pointer to the unit ID of the file.
+ * \param value: `int*` Pointer of the variable to be updated.
+ */
+extern "C" void read_int_(int* uid, int* value);
+/*! \brief Write a complex value to a file.
+ *
+ * \param uid: `int*` Pointer to the unit ID of the file.
+ * \param real: `double*` Pointer to the real part of the value.
+ * \param imag: `double*` Pointer to the imaginary part of the value.
+ */
+extern "C" void write_complex_(int* uid, double* real, double* imag);
+/*! \brief Write a double precision float value to a file.
+ *
+ * \param uid: `int*` Pointer to the unit ID of the file.
+ * \param value: `double*` Pointer to the variable to be written.
+ */
+extern "C" void write_double_(int* uid, double* value);
+/*! \brief Write an integer value to a file.
+ *
+ * \param uid: `int*` Pointer to the unit ID of the file.
+ * \param value: `int*` Pointer to the variable to be written.
+ */
+extern "C" void write_int_(int* uid, int* value);

src/sphere/edfb.cpp

@@ -7,7 +7,8 @@
 #include <cstring>
 #include <iostream>
 #include <fstream>
-#include "List.h"
+#include "../include/file_io.h"
+#include "../include/List.h"
 
 using namespace std;
 
@@ -45,7 +46,6 @@ int main(int argc, char **argv) {
   double *xiv, *wns, *wls, *pus, *evs, *vss;
   string vns[5];
 
-  int max_nsh = 0; // A helper variable to set the size of dc0m
   int ici;
 
   // Input file reading section
@@ -344,7 +344,7 @@ int main(int argc, char **argv) {
       }
       vss = xi_vector.to_array();
       xiv = xi_vector.to_array();
-      pu = xip + wp;
+      pu = xip * wp;
       wn = pu / 3.0e8;
       fprintf(output, "          XIP          WN           WL           PU           EV\n");
       fprintf(output, "     %13.4lE", xip);
@@ -352,7 +352,7 @@ int main(int argc, char **argv) {
       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);
+      fprintf(output, "  SCALE FACTORS XI\n");
       for (int jxi6612 = 1; jxi6612 <= nxi; jxi6612++)
 	fprintf(output, "%5d%13.4lE\n", jxi6612, xiv[jxi6612 - 1]);
       //INXI branch ends here.
@@ -361,7 +361,10 @@ int main(int argc, char **argv) {
   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));
+    string read_format = "";
+    for (int j = 0; j < i; j++) read_format += " %*d";
+    read_format += " %d";
+    sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog + i));
   }
   nshl = new int[nsph];
   ros = new double[nsph];
@@ -376,7 +379,6 @@ int main(int argc, char **argv) {
     ros[i113 - 1] = ros_val * pow(10.0, ros_val_exp);
     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;
@@ -388,28 +390,33 @@ int main(int argc, char **argv) {
   // 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));
+  int uid = 27;
+  string bin_file_name = "c_TEDF";
+  string status = "UNKNOWN";
+  string mode = "UNFORMATTED";
+  open_file_(&uid, bin_file_name.c_str(), status.c_str(), mode.c_str());
+  write_int_(&uid, &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));
+    write_int_(&uid, (iog + iogi));
+  write_double_(&uid, &exdc);
+  write_double_(&uid, &wp);
+  write_double_(&uid, &xip);
+  write_int_(&uid, &idfc);
+  write_int_(&uid, &nxi);
+  for (int xivi = 0; xivi < nxi; xivi++)
+    write_double_(&uid, (xiv + xivi));
   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]));
+    write_int_(&uid, (nshl + i115 -1));
+    write_double_(&uid, (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]));
+      write_double_(&uid, (rcf[i115 - 1] + ins));
   }
   // Remake the dc0m matrix.
-  dc0m = new complex<double>**[max_nsh];
-  for (int dim1 = 0; dim1 < max_nsh; dim1++) {
+  dc0m = new complex<double>**[nsph];
+  for (int dim1 = 0; dim1 < nsph; dim1++) {
     dc0m[dim1] = new complex<double>*[nsph];
     for (int dim2 = 0; dim2 < nxi; dim2++) {
       dc0m[dim1][dim2] = new complex<double>[nxi];
@@ -432,12 +439,11 @@ int main(int argc, char **argv) {
 	// 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));
+	write_complex_(&uid, &dc0_real, &dc0_img);
       }
     }
   }
-  tedf_file.close();
+  close_file_(&uid);
   if (idfc != 0) {
     fprintf(output, "  DIELECTRIC CONSTANTS\n");
     for (int i473 = 1; i473 <= nsph; i473++) {

src/sphere/file_io.f

+      SUBROUTINE OPEN_FILE(UID,NAME,STA,MODE)
+      INTEGER, INTENT(IN):: UID
+      CHARACTER*64, INTENT(IN):: NAME
+      CHARACTER*8, INTENT(IN):: STA
+      CHARACTER*12, INTENT(IN):: MODE
+      OPEN(UNIT=UID, FILE=NAME, STATUS=STA(1:7), FORM=MODE(1:11))
+      END
+      SUBROUTINE CLOSE_FILE(UID)
+      INTEGER, INTENT(IN):: UID
+      CLOSE(UID)
+      END
+      SUBROUTINE READ_INT(UID, VALUE)
+      INTEGER, INTENT(IN):: UID
+      INTEGER, INTENT(OUT):: VALUE
+      READ(UID)VALUE
+      END
+      SUBROUTINE WRITE_COMPLEX(UID, RVAL, IVAL)
+      INTEGER, INTENT(IN):: UID
+      REAL*8, INTENT(IN):: RVAL, IVAL
+      COMPLEX*16:: VALUE
+      VALUE=COMPLEX(RVAL,IVAL)
+      WRITE(UID)VALUE
+      END
+      SUBROUTINE WRITE_DOUBLE(UID, VALUE)
+      INTEGER, INTENT(IN):: UID
+      REAL*8, INTENT(IN):: VALUE
+      WRITE(UID)VALUE
+      END
+      SUBROUTINE WRITE_INT(UID,VALUE)
+      INTEGER, INTENT(IN):: UID
+      INTEGER, INTENT(IN):: VALUE
+      WRITE(UID)VALUE
+      END