From 0618d629d1d591eb41e65128e798863b90ab10b0 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Sun, 24 Sep 2023 20:37:07 +0200
Subject: [PATCH] Replicate edfb.f as edfb.cpp
---
src/sphere/edfb.cpp | 100 +++++++++++++++++++++++++++++++++++++++++---
1 file changed, 95 insertions(+), 5 deletions(-)
diff --git a/src/sphere/edfb.cpp b/src/sphere/edfb.cpp
index 542e1924..1f38c0a1 100644
--- a/src/sphere/edfb.cpp
+++ b/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;
}
--
GitLab