diff --git a/.gitignore b/.gitignore
index 0f5a2cfd44297542c8b1bac2d1e9b277f4753bba..67cb955181b0aa93c80f1590c7d8ed0d73ee27a2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,7 @@
build/include
build/cluster/*
build/sphere/*
+build/testing/*
build/trapping/*
doc/build/*
src/objects/*
\ No newline at end of file
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index dbe3ed281eeaac3547edda212a9679e827d5a055..d52d23a354a847da5a227e182a99e26655038ba2 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -30,6 +30,7 @@ building_stage:
paths:
- build/cluster/*
- build/sphere/*
+ - build/testing/*
- build/trapping/*
- doc/build/*
exclude:
@@ -45,8 +46,20 @@ building_stage:
- echo "Getting system info ..."
- cat /etc/os-release
- cd src
- - echo "Running make ..."
- - make -j
+ - echo "Running make with gnu compilers version 11..."
+ - make wipe && CXX=g++-11 FC=gfortran-11 make -j
+ - echo "Running make with gnu compilers version 12..."
+ - make wipe && CXX=g++-12 FC=gfortran-12 make -j
+ - echo "Running make with gnu fortran version 12 and clang version 13..."
+ - make wipe && CXX="clang++-13 -stdlib=libc++" FC=gfortran-12 make -j
+ - echo "Running make with gnu fortran version 12 and clang version 14..."
+ - make wipe && CXX="clang++-14 -stdlib=libc++" FC=gfortran-12 make -j
+ - echo "Running make with gnu fortran version 12 and clang version 15..."
+ - make wipe && CXX="clang++-15 -stdlib=libc++" FC=gfortran-12 make -j
+ - echo "Running make with gnu fortran version 12 and clang version 16..."
+ - make wipe && CXX="clang++-16 -stdlib=libc++" FC=gfortran-12 make -j
+ - echo "Finally running make with default compilers..."
+ - make wipe && make -j
- make docs -j && make -C ../doc/build/latex -j
running_stage:
diff --git a/.gitlab/merge_request_templates/.gitkeep b/.gitlab/merge_request_templates/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.gitlab/merge_request_templates/default_mr.md b/.gitlab/merge_request_templates/default_mr.md
new file mode 100644
index 0000000000000000000000000000000000000000..b3bdef741d1c4f6fe9b0f23d83b33429159376aa
--- /dev/null
+++ b/.gitlab/merge_request_templates/default_mr.md
@@ -0,0 +1,22 @@
+<!-- title of the merge request -->
+
+# MR Title
+
+<!-- General merge request description.
+
+Write a small paragraph describing the scope of the merge request. Cite relevant issues, if any. Keep in mind that the text might be used to write a release report, so adopt a descriptive approach.
+
+-->
+
+Text to describe the MR.
+
+## Change log
+
+<!-- List of relevant changes.
+
+Provide a list of the main changes introduced by this MR, possibly in the format of a bullet or numbered item list.
+
+-->
+
+1. Change nr. 1
+2. Change nr. 2 ...
\ No newline at end of file
diff --git a/build/README.md b/build/README.md
index 01d8f651732368753e65a697064164fd789ae3ac..c2634408003e16b6280cc37ce3b2668377894d65 100644
--- a/build/README.md
+++ b/build/README.md
@@ -62,10 +62,10 @@ The default behaviour of `np_sphere` is to take the same input files as `edfb_sp
### trapping
-The execution of trapping programs requires at least one of the previous programs to have produced a complete output set. A light-weight trapping calculation has been configured with input and legacy output files stored in the `../../test_data/trapping/` folder. Since the FORTRAN code assumes the input and output to be defined within the program, it is not yet possible to run the *FORTRAN* version on this case, unless the source code is modified accordingly. Conversely, the *C++* version can be executed without the need to modfy and re-compile the code. The work-flow to test trapping is described below.
+The execution of trapping programs requires at least one of the previous programs to have produced a complete output set. A light-weight trapping calculation has been configured with input and legacy output files stored in the `../../test_data/trapping/` folder. Since the FORTRAN code assumes the input and output to be defined within the program, it is not yet possible to run the *FORTRAN* version on this case, unless the source code is modified accordingly. Conversely, the *C++* version can be executed without the need to modify and re-compile the code. The work-flow to test trapping is described below.
1. cd to the `build/sphere` folder.
-2. run `np_sphere` with arguments to take input from the trapping tast and write output in the trapping folder:
+2. run `np_sphere` with arguments to take input from the trapping test and write output in the trapping build folder:
> ./np_sphere ../../test_data/trapping/DEDFB ../../test_data/trapping/DSPH ../trapping
@@ -74,4 +74,4 @@ The execution of trapping programs requires at least one of the previous program
> ./np_trapping ../../test_data/trapping/DFRFME ../../test_data/trapping/DLFFT .
-5. Check the consistency between `np_trapping` output files (`c_out66.txt` and `c_out67.txt`) and the legacy *FORTRAN* output for this case (the files named, respectively, `fort.66` and `fort.67` in the `test_data/trapping/` folder).
+5. Check the consistency between `np_trapping` output files (`c_out66.txt` and `c_out67.txt`) and the legacy *FORTRAN* output for this case (the files named, respectively, `fort.66` and `fort.67` in the `test_data/trapping/` folder). Consider that some of the output values will be affected by numeric noise and take substantially different values. However, this is expected for results whose order of magnitude is clearly below the precision level of the calculation, as they represent results appraching zero that were just approximated with different precision.
diff --git a/build/testing/.gitkeep b/build/testing/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/containers/docker/Dockerfile b/containers/docker/Dockerfile
index 013c483ca45052c15956e6e50409a698bd260c0a..b72c8b113ecc29212abc06ce6ac68a3cda7f7363 100644
--- a/containers/docker/Dockerfile
+++ b/containers/docker/Dockerfile
@@ -6,7 +6,7 @@ WORKDIR /root
# make sure the debian bullseye us up to date, install needed packages
RUN DEBIAN_FRONTEND=noninteractive apt update && DEBIAN_FRONTEND=noninteractive apt -y upgrade
# install packages needed to build binaries
-RUN DEBIAN_FRONTEND=noninteractive apt -y install g++ gfortran make gcc-offload-nvptx libhdf5-dev
+RUN DEBIAN_FRONTEND=noninteractive apt -y install g++ gfortran gcc-offload-nvptx g++-11 gfortran-11 gcc-11-offload-nvptx g++-12 gfortran-12 gcc-12-offload-nvptx clang libc++-dev libc++abi-dev clang-13 clang-14 clang-15 clang-16 libhdf5-dev make
# install packages needed to run python scripts for checks
RUN DEBIAN_FRONTEND=noninteractive apt -y install python3 python-is-python3 python3-regex
# install packages needed to run doxygen to create html docs
diff --git a/src/Makefile b/src/Makefile
index d0d09a0684cbec1fd17cd5421c09762f3ce779c1..73483ad0d279c2a6028f331ab596afec6c53d591 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,4 +1,4 @@
-SUBDIRS := cluster sphere trapping
+SUBDIRS := cluster sphere testing trapping
SRCDIR=$(PWD)
OBJDIR=$(SRCDIR)/objects
ifndef BUILDDIR
diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index d7a38585693a236fb7a84669f120a8d110de90e8..a8da6ac5b9f241d109c7500c44419399f0e28e9a 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -13,7 +13,7 @@ include ../make.inc
F_CLU_OBJS=$(OBJDIR)/clu.o $(OBJDIR)/edfb_clu.o
-CXX_CLU_OBJS=$(OBJDIR)/np_cluster.o $(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/cluster.o
+CXX_CLU_OBJS=$(OBJDIR)/np_cluster.o $(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/cluster.o $(OBJDIR)/TransitionMatrix.o
all: $(BUILDDIR_CLU)/clu $(BUILDDIR_CLU)/edfb_clu $(BUILDDIR_CLU)/np_cluster
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 92646598166084a8288f6343e7d76f94774a19bb..94bec3c35c601f5433b063b2350673b482e33625 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -1,7 +1,9 @@
/*! \file cluster.cpp
+ *
+ * \brief Implementation of the calculation for a cluster of spheres.
*/
-#include <cstdio>
#include <complex>
+#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
@@ -22,6 +24,10 @@
#include "../include/clu_subs.h"
#endif
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#include "../include/TransitionMatrix.h"
+#endif
+
using namespace std;
/*! \brief C++ implementation of CLU
@@ -108,13 +114,11 @@ void cluster(string config_file, string data_file, string output_path) {
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0, lcalc = 0;
complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);
- int max_ici = 0;
- for (int insh = 0; insh < nsph; insh++) {
- int nsh = sconf->nshl_vec[insh];
- int ici = (nsh + 1) / 2;
- if (ici > max_ici) max_ici = ici;
+ int configurations = 0;
+ for (int ci = 1; ci <= nsph; ci++) {
+ if (sconf->iog_vec[ci -1] >= ci) configurations++;
}
- C2 *c2 = new C2(nsph, max_ici, npnt, npntts);
+ C2 *c2 = new C2(nsph, configurations, npnt, npntts);
complex<double> **am = new complex<double>*[mxndm];
for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm]();
const int ndi = c4->nsph * c4->nlim;
@@ -272,7 +276,7 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " \n");
}
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
- printf("INFO: running scale iteration %d...", jxi488);
+ printf("INFO: running scale iteration %d of %d...", jxi488, nxi);
int jaw = 1;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->scale_vec[jxi488 - 1];
@@ -329,30 +333,12 @@ void cluster(string config_file, string data_file, string output_path) {
ztm(am, c1, c1ao, c4, c6, c9);
if (sconf->idfc >= 0) {
if (jxi488 == gconf->jwtm) {
- int is = 1;
- fstream ttms_file;
- string ttms_name = output_path + "/c_TTMS";
- ttms_file.open(ttms_name.c_str(), ios::out | ios::binary);
- if (ttms_file.is_open()) {
- ttms_file.write(reinterpret_cast<char *>(&is), sizeof(int));
- ttms_file.write(reinterpret_cast<char *>(&lm), sizeof(int));
- ttms_file.write(reinterpret_cast<char *>(&vk), sizeof(double));
- ttms_file.write(reinterpret_cast<char *>(&exri), sizeof(double));
- int nlemt = 2 * c4->nlem;
- for (int ami = 0; ami < nlemt; ami++) {
- for (int amj = 0; amj < nlemt; amj++) {
- complex<double> value = c1ao->am0m[ami][amj];
- double rval = value.real();
- double ival = value.imag();
- ttms_file.write(reinterpret_cast<char *>(&rval), sizeof(double));
- ttms_file.write(reinterpret_cast<char *>(&ival), sizeof(double));
- }
- }
- ttms_file.close();
- } else { // Could not open TM file. Should never occur.
- printf("\nERROR: failed to open TTMS file.\n");
- break;
- }
+ int nlemt = 2 * c4->nlem;
+ string ttms_name = output_path + "/c_TTMS.hd5";
+ TransitionMatrix ttms(nlemt, lm, vk, exri, c1ao->am0m);
+ ttms.write_binary(ttms_name, "HDF5");
+ ttms_name = output_path + "/c_TTMS";
+ ttms.write_binary(ttms_name);
}
}
// label 156: continue from here
@@ -515,7 +501,7 @@ void cluster(string config_file, string data_file, string output_path) {
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
- for (int i = 0; i < 3; i++) {
+ for (int i = 0; i < 2; i++) {
double value = c1ao->scscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->scscpm[i].real();
diff --git a/src/cluster/np_cluster.cpp b/src/cluster/np_cluster.cpp
index 0af2b80f8bd8c3514084f30dc77b6fa38c49e909..73caa518f9d330b8ad4114a30f123ba320035756 100644
--- a/src/cluster/np_cluster.cpp
+++ b/src/cluster/np_cluster.cpp
@@ -1,4 +1,19 @@
-/*! \file np_sphere.cpp
+/*! \file np_cluster.cpp
+ *
+ * \brief Cluster of spheres scattering problem handler.
+ *
+ * This program emulates the execution work-flow originally handled by the
+ * FORTRAN EDFB and CLU codes, which undertook the task of solving the
+ * scattering calculation for the case of a cluster of spheres, with one or
+ * more materials. The program is designed to work in two modes: emulation and
+ * enhanced. The emulation mode is activated by invoking the program without
+ * arguments. In this case, the code looks for hard-coded default input and
+ * writes output in the execution folder, replicating the behaviour of the
+ * original FORTRAN code. Advanced mode, instead, is activated by passing
+ * command line arguments that locate the desired input files and a valid
+ * folder to write the output into. The emulation mode is useful for testing,
+ * while the advanced mode implements the possibility to change input and
+ * output options, without having to modify the code.
*/
#include <cstdio>
@@ -20,6 +35,10 @@ extern void cluster(string config_file, string data_file, string output_path);
* or some otherwise formatted configuration data set. The code can be
* linked to a luncher script or to a GUI oriented application that performs
* the configuration and runs the main program.
+ *
+ * \param argc: `int` The number of arguments given in command-line.
+ * \param argv: `char **` The vector of command-line arguments.
+ * \return result: `int` An exit code passed to the OS (0 for succesful execution).
*/
int main(int argc, char **argv) {
string config_file = "../../test_data/cluster/DEDFB";
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 1050a61acaed2d6bf4754df2eeada5fe0ba68a36..37338bdf290661633c6e33d688dcb13008090f7c 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -1,4 +1,28 @@
/*! \file Configuration.h
+ *
+ * \brief Configuration data structures.
+ *
+ * To handle the calculation of a scattering problem, the code needs a set
+ * of configuration parameters that define the properties of the scattering
+ * particle, of the incident radiation field and of the geometry of the
+ * problem. The necessary information is managed through the use of two
+ * classes: `ScattererConfiguration` and `GeometryConfiguration`. The first
+ * class, `ScattererConfiguration`, undertakes the role of describing the
+ * scattering particle properties, by defining its structure in terms of a
+ * distribution of spherical sub-particles, and the optical properties of
+ * the constituting materials. It also defines the scaling vector, which
+ * tracks the ratio of particle size and radiation wavelength, eventually
+ * allowing for the iteration of the calculation on different wavelengths.
+ * Multiple materials and layered structures are allowed, while sphere
+ * compenetration is not handled, due to the implied discontinuity on the
+ * spherical symmetry of the elementary sub-particles. The second class,
+ * `GeometryConfiguration`, on the contrary, describes the properties of
+ * the radiation field impinging on the particle. It defines the incoming
+ * radiation polarization state, the incident direction and the scattering
+ * directions that need to be accounted for. Both classes expose methods to
+ * read the required configuration data from input files formatted according
+ * to the same rules expected by the original code. In addition, they offer
+ * perform I/O operations towards proprietary and standard binary formats.
*/
#ifndef INCLUDE_CONFIGURATION_H_
@@ -189,7 +213,7 @@ class ScattererConfiguration {
friend void cluster(std::string, std::string, std::string);
friend void sphere(std::string, std::string, std::string);
protected:
- //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x LAYERS].
+ //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
std::complex<double> ***dc0_matrix;
//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
double *radii_of_spheres;
@@ -282,7 +306,7 @@ public:
* for dimensions).
* \param dc_matrix: Matrix of reference dielectric constants.
* \param has_external: `bool` Flag to set whether to add an external spherical layer.
- * \param exdc: `double` EXDC
+ * \param exdc: `double` External medium dielectric constant.
* \param wp: `double` wp
* \param xip: `double` xip
*/
@@ -301,7 +325,7 @@ public:
double exdc,
double wp,
double xip
- );
+ );
/*! \brief Destroy a scatterer configuration instance.
*/
@@ -370,6 +394,19 @@ public:
* \param file_name: `string` Name of the file to be written.
*/
void write_formatted(std::string file_name);
+
+ /*! \brief Test whether two instances of ScattererConfiguration are equal.
+ *
+ * ScattererConfiguration objects can be obtained in a variety of manner. They can
+ * be constructed manually, through the class constructor, they can be read from
+ * formatted configuration files, or they can be loaded from binary files. The `==`
+ * operator tests for the equivalence of two configuration instances, returning `true`
+ * if they are equivalnet, `false` otherwise.
+ *
+ * \param other: `ScattererConfiguration &` Reference to the instance to be compared with.
+ * \return result: `bool` True, if the two instances are equal, false otherwise.
+ */
+ bool operator ==(ScattererConfiguration &other);
};
#endif
diff --git a/src/include/List.h b/src/include/List.h
index 0f9ebeab42e3e2e52d434e854420d061cd0c1d3a..8c51bd1ce7bb5816a48bd958c2c6656216e5a0fa 100644
--- a/src/include/List.h
+++ b/src/include/List.h
@@ -1,4 +1,6 @@
/*! \file List.h
+ *
+ * \brief A library of classes used to manage dynamic lists.
*/
#ifndef INCLUDE_LIST_H_
diff --git a/src/include/Parsers.h b/src/include/Parsers.h
index 94dd68e891b175c7263ca45585625e96122e0bf1..5b11d9f6a383f5c941be1912d20c0fa55860ae90 100644
--- a/src/include/Parsers.h
+++ b/src/include/Parsers.h
@@ -1,4 +1,7 @@
/*! \file Parsers.h
+ *
+ * \brief A library of functions designed to parse formatted input
+ * into memory.
*/
#ifndef FILE_NOT_FOUND_ERROR
diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..9180fabaea8c937227d5574ea0a821ee9a82f2aa
--- /dev/null
+++ b/src/include/TransitionMatrix.h
@@ -0,0 +1,177 @@
+/*! \file TransitionMatrix.h
+ *
+ * \brief Representation of the Transition Matrix.
+ */
+
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#define INCLUDE_TRANSITIONMATRIX_H_
+
+/**
+ * \brief Exception for unrecognized file formats.
+ */
+class UnrecognizedFormatException: public std::exception {
+protected:
+ //! Description of the problem.
+ std::string message;
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param problem: `string` Description of the problem that occurred.
+ */
+ UnrecognizedFormatException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/*! \brief Class to represent the Transition Matrix.
+ */
+class TransitionMatrix {
+ protected:
+ //! Matrix type identifier.
+ int is;
+ //! Maximum field expansion order.
+ int l_max;
+ //! Wave number in scale units.
+ double vk;
+ //! External medium refractive index.
+ double exri;
+ //! Vectorized matrix elements.
+ std::complex<double> *elements;
+ //! Sphere radius.
+ double sphere_radius;
+ //! Matrix shape
+ int *shape;
+
+ /*! \brief Build transition matrix from a HDF5 binary input file.
+ *
+ * This function takes care of the specific task of building a transition
+ * matrix memory data structure from a HDF5 binary input file.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \return config: `TransitionMatrix *` Pointer to object containing the
+ * transition matrix data.
+ */
+ static TransitionMatrix *from_hdf5(std::string file_name);
+
+ /*! \brief Build transition matrix from a legacy binary input file.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \return config: `TransitionMatrix *` Pointer to object containing the
+ * transition matrix data.
+ */
+ static TransitionMatrix *from_legacy(std::string file_name);
+
+ /*! \brief Write the Transition Matrix to HDF5 binary output.
+ *
+ * This function takes care of the specific task of building a transition
+ * matrix memory data structure from a binary input file formatted according
+ * to the structure used by the original FORTRAN code.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Write the Transition Matrix to legacy binary output.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_legacy(std::string file_name);
+ public:
+ /*! \brief Default Transition Matrix instance constructor.
+ *
+ * \param _is: `int` Matrix type identifier
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _vk: `double` Wave number in scale units.
+ * \param _exri: `double` External medium refractive index.
+ * \param _elements: `complex<double> *` Vectorized elements of the matrix.
+ * \param _radius: `double` Radius for the single sphere case (defaults to 0.0).
+ */
+ TransitionMatrix(
+ int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ double _radius=0.0
+ );
+
+ /*! \brief Transition Matrix instance constructor for single sphere.
+ *
+ * This constructor allocates the memory structure needed to represent the transition
+ * matrix for the case of a single sphere.
+ *
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _vk: `double` Wave number in scale units.
+ * \param _exri: `double` External medium refractive index.
+ * \param _rmi: Matrix of complex.
+ * \param _rei: Matrix of complex.
+ * \param _sphere_radius: `double` Radius of the sphere.
+ */
+ TransitionMatrix(
+ int _lm, double _vk, double _exri, std::complex<double> **_rmi,
+ std::complex<double> **_rei, double _sphere_radius
+ );
+
+ /*! \brief Transition Matrix instance constructor for a cluster of spheres.
+ *
+ * This constructor allocates the memory structure needed to represent the transition
+ * matrix for the case of a cluster of spheres.
+ *
+ * \param _nlemt: `int` Size of the matrix (2 * LE * (LE + 2)).
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _vk: `double` Wave number in scale units.
+ * \param _exri: `double` External medium refractive index.
+ * \param _am0m: Matrix of complex.
+ */
+ TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, std::complex<double> **_am0m);
+
+ /*! \brief Transition Matrix instance destroyer.
+ */
+ ~TransitionMatrix();
+
+ /*! \brief Build transition matrix from binary input file.
+ *
+ * In some cases, it is necessary to perform calculations starting from a pre-computed
+ * transition matrix. If this matrix is not available in memory (e.g. because it was
+ * calculated by a different process), but it was saved in a binary file, this function
+ * can be used to load it back in memory. The function operates in two modes: "LEGACY",
+ * which reads the matrix data from a proprietary binary file, having the same structure
+ * as the one used by the original FORTRAN code, and "HDF5", which, instead, reads the
+ * data from an input file complying with the HDF5 format.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"]. Optional
+ * (default is "LEGACY").
+ * \return config: `TransitionMatrix *` Pointer to object containing the transition
+ * matrix data.
+ */
+ static TransitionMatrix* from_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Write the Transition Matrix to a binary file.
+ *
+ * This function writes a hard-copy of the transition matrix to an output file, making
+ * it available for subsequent processes to reload. The function operates in two modes:
+ * "LEGACY", which writes a proprietary binary file, using the same structure of the
+ * original FORTRAN code, and "HDF5", which, instead, writes the output to a file using
+ * the HDF5 format, thus leaving it available for inspection with external tools.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
+ * (default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Test whether two instances of TransitionMatrix are equal.
+ *
+ * Transition matrixes can be the result of a calculation or of a file input operation,
+ * reading from a previously computed object. The `==` operator tests for the equivalence
+ * of two transition matrixes, returning `true` if they are equivalnet, `false` otherwise.
+ *
+ * \param other: `TransitionMatrix &` Reference to the instance to be compared with.
+ * \return result: `bool` True, if the two instances are equal, false otherwise.
+ */
+ bool operator ==(TransitionMatrix &other);
+};
+
+#endif
diff --git a/src/include/file_io.h b/src/include/file_io.h
index 8c228098e872d5941bb7b9b9f51ca3576ffa9b71..a39401dd1fcaea4473543adbc166ad278e260471 100644
--- a/src/include/file_io.h
+++ b/src/include/file_io.h
@@ -29,6 +29,7 @@ class FileSchema {
*
* \param num_rec: `int` Number of records in the file.
* \param rec_types: `string *` Description of the records in the file.
+ * \param rec_names: `string *` Names of the records in the file.
*/
FileSchema(int num_rec, std::string *rec_types, std::string *rec_names=NULL);
@@ -82,8 +83,8 @@ class HDFFile {
* \param fapl_id: `hid_t` File access property list identifier (default is `H5P_DEFAULT`).
*/
HDFFile(
- std::string name, unsigned int flags = H5F_ACC_EXCL,
- hid_t fcpl_id = H5P_DEFAULT, hid_t fapl_id = H5P_DEFAULT
+ std::string name, unsigned int flags=H5F_ACC_EXCL,
+ hid_t fcpl_id=H5P_DEFAULT, hid_t fapl_id=H5P_DEFAULT
);
/*! \brief HDFFile instance destroyer.
@@ -104,8 +105,8 @@ class HDFFile {
* \return hdf_file: `HDFFile *` Pointer to a new, open HDF5 file.
*/
static HDFFile* from_schema(
- FileSchema &schema, std::string name, unsigned int flags = H5F_ACC_EXCL,
- hid_t fcpl_id = H5P_DEFAULT, hid_t fapl_id = H5P_DEFAULT
+ FileSchema &schema, std::string name, unsigned int flags=H5F_ACC_EXCL,
+ hid_t fcpl_id=H5P_DEFAULT, hid_t fapl_id=H5P_DEFAULT
);
/*! \brief Get current status.
@@ -120,6 +121,23 @@ class HDFFile {
*/
bool is_open() { return file_open_flag; }
+ /*! \brief Read data from attached file.
+ *
+ * \param dataset_name: `string` Name of the dataset to read from.
+ * \param data_type: `string` Memory data type identifier.
+ * \param buffer: `hid_t` Starting address of the memory sector to store the data.
+ * \param mem_space_id: `hid_t` Memory data space identifier (defaults to `H5S_ALL`).
+ * \param file_space_id: `hid_t` File space identifier (defaults to `H5S_ALL`).
+ * \param dapl_id: `hid_t` Data access property list identifier (defaults to `H5P_DEFAULT`).
+ * \param dxpl_id: `hid_t` Data transfer property list identifier (defaults to `H5P_DEFAULT`).
+ * \return status: `herr_t` Exit status of the operation.
+ */
+ herr_t read(
+ std::string dataset_name, std::string data_type, void *buffer,
+ hid_t mem_space_id=H5S_ALL, hid_t file_space_id=H5S_ALL,
+ hid_t dapl_id=H5P_DEFAULT, hid_t dxpl_id=H5P_DEFAULT
+ );
+
/*! \brief Write data to attached file.
*
* \param dataset_name: `string` Name of the dataset to write to.
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index f0afd2dff43f125b98dbaf4da5126da62991d37f..9526b34341f85fc4a8be09a3ba259b655a05982a 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -35,7 +35,7 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
* backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
*
* \param n: `int` Order of the function.
- * \param z: `complex<double>` Argumento of the function.
+ * \param z: `complex<double>` Argument of the function.
* \param nm: `int &` Highest computed order.
* \param csj: Vector of complex. The desired function \f$j\f$.
*/
@@ -140,7 +140,7 @@ int msta1(double x, int mp);
* This function determines the starting point for backward recurrence such that
* all \f$J\f$ and \f$j\f$ functions have `mp` significant digits.
*
- * \param x: `double` Absolute value of the argumetn to \f$J\f$ or \$j\f$.
+ * \param x: `double` Absolute value of the argumetn to \f$J\f$ or \f$j\f$.
* \param n: `int` Order of the function.
* \param mp: `int` Requested number of significant digits.
* \return result: `int` The necessary starting point.
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 33d9781bed4c081c8ee34e2e7d322034f6148350..82594bd0744878f8109bf16003686c2f992a7341 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -97,7 +97,6 @@ void ffrf(
* \param ffte: `double *`. QUESTION: definition?
* \param ffts: `double *`. QUESTION: definition?
* \param cil: `CIL *` Pointer to a CIL structure.
- * \param ccr: `CCR *` Pointer to a CCR structure.
*/
void ffrt(
std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 8031eaad68f3127546b222eb3d0e62bff6f5e0be..81cb853bfcb1f757cab084a1f587a620c9b14f0d 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1,14 +1,14 @@
/*! \file Commons.cpp
*
- * DEVELOPMENT NOTE:
- * The construction of common blocks requires some information
- * that is stored in configuration objects and is needed to compute
- * the allocation size of vectors and matrices. Currently, this
- * information is passed as arguments to the constructors of the
- * common blocks. A simpler and more logical way to operate is
- * to design the constructors to take as arguments only pointers
- * to the configuration objects. These, on their turn, need to
- * expose methods to access the relevant data in read-only mode.
+ * DEVELOPMENT NOTE:
+ * The construction of common blocks requires some information
+ * that is stored in configuration objects and is needed to compute
+ * the allocation size of vectors and matrices. Currently, this
+ * information is passed as arguments to the constructors of the
+ * common blocks. A simpler and more logical way to operate would
+ * be to design the constructors to take as arguments only pointers
+ * to the configuration objects. These, on their turn, need to
+ * expose methods to access the relevant data in read-only mode.
*/
#include <complex>
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 5d1c05d898179cca789efc59061ee3600a8a899c..7eb58c4c0bd5738d7632beb650fa943e13f0efad 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -1,4 +1,6 @@
/*! \file Configuration.cpp
+ *
+ * \brief Implementation of the configuration classes.
*/
#include <cmath>
@@ -6,9 +8,9 @@
#include <cstdio>
#include <exception>
#include <fstream>
+#include <hdf5.h>
#include <regex>
#include <string>
-#include <hdf5.h>
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
@@ -37,7 +39,7 @@ GeometryConfiguration::GeometryConfiguration(
double in_ph_start, double in_ph_step, double in_ph_end,
double sc_ph_start, double sc_ph_step, double sc_ph_end,
int _jwtm
- ) {
+) {
number_of_spheres = _nsph;
l_max = _lm;
in_pol = _in_pol;
@@ -182,7 +184,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
in_ph_start, in_ph_step, in_ph_end,
sc_ph_start, sc_ph_step, sc_ph_end,
_jwtm
- );
+ );
return conf;
}
@@ -201,7 +203,7 @@ ScattererConfiguration::ScattererConfiguration(
double ex,
double w,
double x
- ) {
+) {
number_of_spheres = nsph;
number_of_scales = nxi;
reference_variable_name = variable_name;
@@ -232,21 +234,21 @@ ScattererConfiguration::ScattererConfiguration(
}
ScattererConfiguration::~ScattererConfiguration() {
- int max_ici = 0;
+ int configurations = 0;
for (int i = 1; i <= number_of_spheres; i++) {
if (iog_vec[i - 1] < i) continue;
- int nsh = nshl_vec[i - 1];
- if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+ configurations++;
}
- for (int i = 0; i < max_ici; i++) {
+ for (int i = 0; i < configurations; i++) {
for (int j = 0; j < number_of_spheres; j++) {
delete[] dc0_matrix[i][j];
}
}
delete[] dc0_matrix;
- for (int i = 0; i < number_of_spheres; i++) {
+ for (int i = 0; i < configurations; i++) {
delete[] rcf[i];
}
+ delete[] rcf;
delete[] nshl_vec;
delete[] radii_of_spheres;
delete[] iog_vec;
@@ -266,112 +268,6 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
return conf;
}
-ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
- ScattererConfiguration *conf = NULL;
- return conf;
-}
-
-ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
- int nsph;
- int *iog;
- double _exdc, _wp, _xip;
- int _idfc, nxi;
- int *nshl_vector;
- double *xi_vec;
- double *ros_vector;
- double **rcf_vector;
- complex<double> ***dc0m;
-
- int max_ici = 0;
- fstream input;
- input.open(file_name.c_str(), ios::in | ios::binary);
- if (input.is_open()) {
- input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
- iog = new int[nsph]();
- for (int i = 0; i < nsph; i++) {
- input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
- }
- input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
- input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
- input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
- input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
- input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
- try {
- xi_vec = new double[nxi]();
- } catch (const bad_alloc &ex) {
- throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + to_string(nxi));
- }
- for (int i = 0; i < nxi; i++) {
- input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
- }
- nshl_vector = new int[nsph]();
- ros_vector = new double[nsph]();
- rcf_vector = new double*[nsph];
- for (int i115 = 1; i115 <= nsph; i115++) {
- if (iog[i115 - 1] < i115) {
- rcf_vector[i115 - 1] = new double[1]();
- continue;
- }
- input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int));
- input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double));
- int nsh = nshl_vector[i115 - 1];
- if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
- try {
- rcf_vector[i115 - 1] = new double[nsh]();
- } catch (const bad_alloc &ex) {
- throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + to_string(nsh));
- }
- for (int nsi = 0; nsi < nsh; nsi++) {
- input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
- }
- }
- dc0m = new complex<double>**[max_ici];
- int dim3 = (_idfc == 0) ? nxi : 1;
- for (int dim1 = 0; dim1 < max_ici; dim1++) {
- dc0m[dim1] = new complex<double>*[nsph];
- for (int dim2 = 0; dim2 < nsph; dim2++) {
- dc0m[dim1][dim2] = new complex<double>[dim3]();
- }
- }
- 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;
- int nsh = nshl_vector[i162 - 1];
- int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
- for (int i157 = 0; i157 < ici; i157++) {
- double dc0_real, dc0_img;
- input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
- input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
- dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
- }
- }
- }
- input.close();
- } else { // Opening of the input file did not succeed.
- OpenConfigurationFileException ex(file_name);
- throw ex;
- }
-
- ScattererConfiguration *conf = new ScattererConfiguration(
- nsph,
- xi_vec,
- nxi,
- "XIV",
- iog,
- ros_vector,
- nshl_vector,
- rcf_vector,
- _idfc,
- dc0m,
- false,
- _exdc,
- _wp,
- _xip
- );
- return conf;
-}
-
ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_name) {
int num_lines = 0;
int last_read_line = 0;
@@ -385,7 +281,6 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
throw ex;
}
int nsph, ies;
- int max_ici = 0;
re = regex("[0-9]+");
string str_target = file_lines[last_read_line];
for (int ri = 0; ri < 2; ri++) {
@@ -534,88 +429,83 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
}
}
int *iog_vector = new int[nsph]();
- double *ros_vector = new double[nsph]();
- double **rcf_vector = new double*[nsph];
- int *nshl_vector = new int[nsph]();
- int filled_iogs = 0;
+ str_target = file_lines[++last_read_line];
re = regex("[0-9]+");
- while (filled_iogs < nsph) {
- string str_target = file_lines[++last_read_line];
- while(regex_search(str_target, m, re)) {
- iog_vector[filled_iogs++] = stoi(m.str());
+ int configurations = 0;
+ for (int i = 1; i <= nsph; i++) {
+ bool success = regex_search(str_target, m, re);
+ if (success) {
+ iog_vector[i - 1] = stoi(m.str());
str_target = m.suffix().str();
- if (filled_iogs == nsph) break;
+ if (iog_vector[i - 1] >= i) configurations++;
+ } else {
+ str_target = file_lines[++last_read_line];
+ i--;
}
}
+ int index = 0;
+ double *ros_vector = new double[configurations]();
+ int *nshl_vector = new int[configurations]();
+ double **rcf_vector = new double*[configurations];
for (int i113 = 1; i113 <= nsph; i113++) {
- int i_val, nsh;
- double ros_val;
- if (iog_vector[i113 - 1] < i113) {
- rcf_vector[i113 - 1] = new double[1]();
- continue;
- }
- re = regex("[0-9]+");
+ if (iog_vector[i113 - 1] < i113) continue;
str_target = file_lines[++last_read_line];
+ re = regex("[0-9]+");
regex_search(str_target, m, re);
- i_val = stoi(m.str());
- str_target = m.suffix();
+ nshl_vector[i113 - 1] = stoi(m.str());
+ str_target = m.suffix().str();
re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
regex_search(str_target, m, re);
string str_number = m.str();
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
- ros_val = stod(str_number);
- nshl_vector[i113 - 1] = i_val;
- if (max_ici < (i_val + 1) / 2) max_ici = (i_val + 1) / 2;
- ros_vector[i113 - 1] = ros_val;
- nsh = nshl_vector[i113 - 1];
+ ros_vector[i113 - 1] = stod(str_number);
+ int nsh = nshl_vector[i113 - 1];
if (i113 == 1) nsh += ies;
rcf_vector[i113 - 1] = new double[nsh]();
- for (int ns = 0; ns < nsh; ns++) {
- double ns_rcf;
+ for (int ns112 = 0; ns112 < nsh; ns112++) {
str_target = file_lines[++last_read_line];
regex_search(str_target, m, re);
str_number = m.str();
str_number = regex_replace(str_number, regex("D"), "e");
str_number = regex_replace(str_number, regex("d"), "e");
- ns_rcf = stod(str_number);
- rcf_vector[i113 - 1][ns] = ns_rcf;
- }
- }
- complex<double> ***dc0m = new complex<double>**[max_ici];
- int dim3 = (_idfc == 0) ? nxi : 1;
- for (int dim1 = 0; dim1 < max_ici; dim1++) {
- dc0m[dim1] = new complex<double>*[nsph];
- for (int dim2 = 0; dim2 < nsph; dim2++) {
- dc0m[dim1][dim2] = new complex<double>[dim3]();
- }
- }
- re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
+ rcf_vector[i113 - 1][ns112] = stod(str_number);
+ } // ns112 loop
+ } // i113 loop
+
+ complex<double> ***dc0m = new complex<double>**[configurations];
+ complex<double> *dc0 = new complex<double>[configurations]();
+ for (int di = 0; di < configurations; di++) {
+ dc0m[di] = new complex<double>*[nsph];
+ for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[nxi]();
+ } // di loop
for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
- if (_idfc != 0 && jxi468 > 1) continue;
+ if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop
for (int i162 = 1; i162 <= nsph; i162++) {
- if (iog_vector[i162 - 1] < i162) continue;
- int nsh = nshl_vector[i162 - 1];
- int 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;
- str_target = file_lines[++last_read_line];
- regex_search(str_target, m, re);
- string str_number = m.str();
- str_number = regex_replace(str_number, regex("D"), "e");
- str_number = regex_replace(str_number, regex("d"), "e");
- dc0_real = stod(str_number);
- str_target = m.suffix().str();
- regex_search(str_target, m, re);
- str_number = m.str();
- str_number = regex_replace(str_number, regex("D"), "e");
- str_number = regex_replace(str_number, regex("d"), "e");
- dc0_img = stod(str_number);
- dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(dc0_real, dc0_img);
+ if (iog_vector[i162 - 1] >= i162) {
+ int nsh = nshl_vector[i162 - 1];
+ int ici = (nsh + 1) / 2;
+ if (i162 == 1) ici += ies;
+ for (int ic157 = 0; ic157 < ici; ic157++) {
+ str_target = file_lines[++last_read_line];
+ regex_search(str_target, m, re);
+ string str_number = m.str();
+ str_number = regex_replace(str_number, regex("D"), "e");
+ str_number = regex_replace(str_number, regex("d"), "e");
+ double rval = stod(str_number);
+ str_target = m.suffix().str();
+ regex_search(str_target, m, re);
+ str_number = m.str();
+ str_number = regex_replace(str_number, regex("D"), "e");
+ str_number = regex_replace(str_number, regex("d"), "e");
+ double ival = stod(str_number);
+ dc0[ic157] = complex<double>(rval, ival);
+ dc0m[ic157][i162 - 1][jxi468 - 1] = dc0[ic157];
+ } // ic157 loop
}
- }
- }
+ } // i162 loop
+ } // jxi468 loop
+ delete[] dc0;
ScattererConfiguration *config = new ScattererConfiguration(
nsph,
@@ -632,30 +522,210 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
_exdc,
_wp,
_xip
- );
+ );
delete[] file_lines;
return config;
}
+ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
+ ScattererConfiguration *conf = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ string str_name, str_type;
+ if (status == 0) {
+ int nsph;
+ int *iog;
+ double _exdc, _wp, _xip;
+ int _idfc, nxi;
+ int *nshl_vector;
+ double *xi_vec;
+ double *ros_vector;
+ double **rcf_vector;
+ complex<double> ***dc0m;
+ status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
+ status = hdf_file->read("EXDC", "FLOAT64_(1)", &_exdc);
+ status = hdf_file->read("WP", "FLOAT64_(1)", &_wp);
+ status = hdf_file->read("XIP", "FLOAT64_(1)", &_xip);
+ status = hdf_file->read("IDFC", "INT32_(1)", &_idfc);
+ status = hdf_file->read("NXI", "INT32_(1)", &nxi);
+ iog = new int[nsph];
+ string str_type = "INT32_(" + to_string(nsph) + ")";
+ status = hdf_file->read("IOGVEC", str_type, iog);
+ int configuration_count = 0;
+ for (int ci = 0; ci < nsph; ci++) {
+ if (iog[ci] < ci + 1) continue;
+ configuration_count++;
+ }
+ nshl_vector = new int[configuration_count]();
+ ros_vector = new double[configuration_count]();
+ rcf_vector = new double*[configuration_count];
+ for (int i115 = 1; i115 <= nsph; i115++) {
+ if (iog[i115 - 1] < i115) continue;
+ str_name = "NSHL_" + to_string(iog[i115 - 1]);
+ str_type = "INT32_(1)";
+ status = hdf_file->read(str_name, str_type, (nshl_vector + i115 - 1));
+ str_name = "ROS_" + to_string(iog[i115 - 1]);
+ str_type = "FLOAT64_(1)";
+ status = hdf_file->read(str_name, str_type, (ros_vector + i115 - 1));
+ int nsh = nshl_vector[i115 - 1];
+ rcf_vector[i115 - 1] = new double[nsh]();
+ str_name = "RCF_" + to_string(iog[i115 - 1]);
+ str_type = "FLOAT64_(" + to_string(nsh) + ")";
+ status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1]));
+ }
+ xi_vec = new double[nxi];
+ str_name = "XIVEC";
+ str_type = "FLOAT64_(" + to_string(nxi) + ")";
+ status = hdf_file->read(str_name, str_type, xi_vec);
+
+ int dim3 = (_idfc == 0) ? nxi : 1;
+ int element_size = 2 * dim3 * nsph * configuration_count;
+ double *elements = new double[element_size]();
+ str_name = "DC0M";
+ str_type = "FLOAT64_(" + to_string(element_size) + ")";
+ status = hdf_file->read(str_name, str_type, elements);
+ dc0m = new complex<double>**[configuration_count];
+ for (int di = 0; di < configuration_count; di++) {
+ dc0m[di] = new complex<double>*[nsph];
+ for (int dj = 0; dj < nsph; dj++) {
+ dc0m[di][dj] = new complex<double>[dim3]();
+ for (int dk = 0; dk < dim3; dk++) {
+ double rval = elements[(di * nsph) + (dj * dim3) + 2 * dk];
+ double ival = elements[(di * nsph) + (dj * dim3) + 2 * dk + 1];
+ dc0m[di][dj][dk] = complex<double>(rval, ival);
+ }
+ } // dj loop
+ } // di loop
+ delete[] elements;
+ status = hdf_file->close();
+ conf = new ScattererConfiguration(
+ nsph,
+ xi_vec,
+ nxi,
+ "XIV",
+ iog,
+ ros_vector,
+ nshl_vector,
+ rcf_vector,
+ _idfc,
+ dc0m,
+ false,
+ _exdc,
+ _wp,
+ _xip
+ );
+ }
+
+ return conf;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
+ int _nsph;
+ int *_iog_vec;
+ int _ies;
+ double _exdc, _wp, _xip;
+ int _idfc, _nxi;
+ int *_nshl_vec;
+ double *_xi_vec;
+ double *_ros_vec;
+ double **_rcf_vec;
+ complex<double> ***_dc0m;
+
+ fstream input;
+ input.open(file_name.c_str(), ios::in | ios::binary);
+ input.read(reinterpret_cast<char *>(&_nsph), sizeof(int));
+ input.read(reinterpret_cast<char *>(&_ies), sizeof(int));
+ _iog_vec = new int[_nsph]();
+ for (int i = 0; i < _nsph; i++)
+ input.read(reinterpret_cast<char *>(&(_iog_vec[i])), sizeof(int));
+ input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
+ input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
+ input.read(reinterpret_cast<char *>(&_nxi), sizeof(int));
+ _xi_vec = new double[_nxi]();
+ for (int i = 0; i < _nxi; i++)
+ input.read(reinterpret_cast<char *>(&(_xi_vec[i])), sizeof(double));
+ int configurations = 0;
+ for (int i = 1; i <= _nsph; i++) {
+ if (_iog_vec[i - 1] >= i) configurations++;
+ }
+ _nshl_vec = new int[configurations]();
+ _ros_vec = new double[configurations]();
+ _rcf_vec = new double*[configurations];
+ for (int i115 = 1; i115 <= _nsph; i115++) {
+ if (_iog_vec[i115 - 1] < i115) continue;
+ input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));
+ input.read(reinterpret_cast<char *>(&(_ros_vec[i115 - 1])), sizeof(double));
+ int nsh = _nshl_vec[i115 - 1];
+ if (i115 == 1) nsh += _ies;
+ _rcf_vec[i115 - 1] = new double[nsh]();
+ for (int nsi = 0; nsi < nsh; nsi++)
+ input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));
+ }
+ _dc0m = new complex<double>**[configurations];
+ for (int di = 0; di < configurations; di++) {
+ _dc0m[di] = new complex<double>*[_nsph];
+ for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[_nxi]();
+ } // di loop
+ for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
+ if (_idfc != 0 && jxi468 > 1) continue;
+ for (int i162 = 1; i162 <= _nsph; i162++) {
+ if (_iog_vec[i162 - 1] < i162) continue;
+ int nsh = _nshl_vec[i162 - 1];
+ int ici = (nsh + 1) / 2;
+ if (i162 == 1) ici = ici + _ies;
+ for (int i157 = 0; i157 < ici; i157++) {
+ double rval, ival;
+ input.read(reinterpret_cast<char *>(&rval), sizeof(double));
+ input.read(reinterpret_cast<char *>(&ival), sizeof(double));
+ _dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(rval, ival);
+ }
+ }
+ }
+ input.close();
+
+ ScattererConfiguration *conf = new ScattererConfiguration(
+ _nsph,
+ _xi_vec,
+ _nxi,
+ "XIV",
+ _iog_vec,
+ _ros_vec,
+ _nshl_vec,
+ _rcf_vec,
+ _idfc,
+ _dc0m,
+ (_ies == 1),
+ _exdc,
+ _wp,
+ _xip
+ );
+ return conf;
+}
+
void ScattererConfiguration::print() {
int ies = (use_external_sphere)? 1 : 0;
- int max_ici = 0;
+ int configurations = 0;
+ for (int ci = 1; ci <= number_of_spheres; ci++) {
+ if (iog_vec[ci - 1] >= ci) configurations++;
+ }
printf("### CONFIGURATION DATA ###\n");
printf("NSPH = %d\n", number_of_spheres);
printf("ROS = [");
- for (int i = 0; i < number_of_spheres; i++) printf("\t%lg", radii_of_spheres[i]);
+ for (int i = 0; i < configurations; i++) printf("\t%lg", radii_of_spheres[i]);
printf("\t]\n");
printf("IOG = [");
for (int i = 0; i < number_of_spheres; i++) printf("\t%d", iog_vec[i]);
printf("\t]\n");
printf("NSHL = [");
- for (int i = 0; i < number_of_spheres; i++) printf("\t%d", nshl_vec[i]);
+ for (int i = 0; i < configurations; i++) printf("\t%d", nshl_vec[i]);
printf("\t]\n");
printf("RCF = [\n");
- for (int i = 1; i <= number_of_spheres; i++) {
+ for (int i = 1; i <= configurations; i++) {
int nsh = nshl_vec[i - 1];
if (i == 1) nsh += ies;
- if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
printf(" [");
for (int ns = 0; ns < nsh; ns++) {
printf("\t%lg", rcf[i - 1][ns]);
@@ -669,7 +739,7 @@ void ScattererConfiguration::print() {
for (int i = 0; i < number_of_scales; i++) printf("\t%lg", scale_vec[i]);
printf("\t]\n");
printf("DC0M = [\n");
- for (int i = 0; i < max_ici; i++) {
+ for (int i = 0; i < configurations; i++) {
printf(" [\n");
for (int j = 0; j < number_of_spheres; j++) {
printf(" [");
@@ -696,14 +766,15 @@ void ScattererConfiguration::write_binary(string file_name, string mode) {
}
void ScattererConfiguration::write_hdf5(string file_name) {
- const double two_pi = acos(0.0) * 4.0;
- const double evc = 6.5821188e-16;
int ies = (use_external_sphere)? 1 : 0;
- int max_ici = 0;
List<string> rec_name_list(1);
List<string> rec_type_list(1);
List<void *> rec_ptr_list(1);
string str_type, str_name;
+ int configurations = 0;
+ for (int ci = 1; ci <= number_of_spheres; ci++) {
+ if(iog_vec[ci - 1] >= ci) configurations++;
+ }
rec_name_list.set(0, "NSPH");
rec_type_list.set(0, "INT32_(1)");
rec_ptr_list.set(0, &number_of_spheres);
@@ -735,30 +806,29 @@ void ScattererConfiguration::write_hdf5(string file_name) {
str_name = "NSHL_" + to_string(i115);
rec_name_list.append(str_name);
rec_type_list.append("INT32_(1)");
- rec_ptr_list.append(&(nshl_vec[i115 - 1]));
+ rec_ptr_list.append(&(nshl_vec[i115 - 1])); // was not from IOG
str_name = "ROS_" + to_string(i115);
rec_name_list.append(str_name);
rec_type_list.append("FLOAT64_(1)");
- rec_ptr_list.append(&(radii_of_spheres[i115 - 1]));
- int nsh = nshl_vec[i115 - 1];
+ rec_ptr_list.append(&(radii_of_spheres[i115 - 1])); // was not from IOG
+ int nsh = nshl_vec[i115 - 1]; // was not from IOG
if (i115 == 1) nsh += ies;
- if (max_ici < (nsh + 1) / 2) max_ici = nsh + 1 / 2;
- str_name = "RCF_" + to_string(i115);
+ str_name = "RCF_" + to_string(i115); // was not from IOG
str_type = "FLOAT64_(" + to_string(nsh) + ")";
rec_name_list.append(str_name);
rec_type_list.append(str_type);
- rec_ptr_list.append(&(rcf[i115 - 1][0]));
+ rec_ptr_list.append(&(rcf[i115 - 1][0])); // was not from IOG
}
int dim3 = (idfc == 0) ? number_of_scales : 1;
- int dc0m_size = 2 * dim3 * number_of_spheres * max_ici;
+ int dc0m_size = 2 * dim3 * number_of_spheres * configurations;
double *dc0m = new double[dc0m_size];
int dc0_index = 0;
for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
if (idfc != 0 && jxi468 > 1) continue;
for (int i162 = 1; i162 <= number_of_spheres; i162++) {
if (iog_vec[i162 - 1] < i162) continue;
- int nsh = nshl_vec[i162 - 1];
+ int nsh = nshl_vec[i162 - 1]; // was not from IOG
int ici = (nsh + 1) / 2;
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
@@ -780,7 +850,7 @@ void ScattererConfiguration::write_hdf5(string file_name) {
void **rec_pointers = rec_ptr_list.to_array();
const int rec_num = rec_name_list.length();
FileSchema schema(rec_num, rec_types, rec_names);
- HDFFile *hdf_file = HDFFile::from_schema(schema, "c_TEDF.hd5", H5F_ACC_TRUNC);
+ HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
for (int ri = 0; ri < rec_num; ri++)
hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
hdf_file->close();
@@ -794,12 +864,11 @@ void ScattererConfiguration::write_hdf5(string file_name) {
}
void ScattererConfiguration::write_legacy(string file_name) {
- const double two_pi = acos(0.0) * 4.0;
- const double evc = 6.5821188e-16;
fstream output;
int ies = (use_external_sphere)? 1 : 0;
output.open(file_name.c_str(), ios::out | ios::binary);
output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
+ output.write(reinterpret_cast<char *>(&ies), sizeof(int));
for (int i = 0; i < number_of_spheres; i++)
output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
@@ -811,18 +880,18 @@ void ScattererConfiguration::write_legacy(string file_name) {
output.write(reinterpret_cast<char *>(&(scale_vec[i])), sizeof(double));
for (int i115 = 1; i115 <= number_of_spheres; i115++) {
if (iog_vec[i115 - 1] < i115) continue;
- output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int));
- output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double));
- int nsh = nshl_vec[i115 - 1];
+ output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int)); // was not from IOG
+ output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double)); // was not from IOG
+ int nsh = nshl_vec[i115 - 1]; // was not from IOG
if (i115 == 1) nsh += ies;
for (int nsi = 0; nsi < nsh; nsi++)
- output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
+ output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double)); // was not from IOG
}
for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
if (idfc != 0 && jxi468 > 1) continue;
for (int i162 = 1; i162 <= number_of_spheres; i162++) {
if (iog_vec[i162 - 1] < i162) continue;
- int nsh = nshl_vec[i162 - 1];
+ int nsh = nshl_vec[i162 - 1]; // was not from IOG
int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
if (i162 == 1) ici = ici + ies;
for (int i157 = 0; i157 < ici; i157++) {
@@ -878,7 +947,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
pu_vec[i],
ev_vec[i]
- );
+ );
}
break;
case 1:
@@ -898,7 +967,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
pu_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 2:
@@ -918,7 +987,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
pu_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 3:
@@ -938,7 +1007,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
ev_vec[i],
xi_vec[i]
- );
+ );
}
break;
case 4:
@@ -958,7 +1027,7 @@ void ScattererConfiguration::write_formatted(string file_name) {
wl_vec[i],
pu_vec[i],
xi_vec[i]
- );
+ );
}
break;
default:
@@ -1014,3 +1083,60 @@ void ScattererConfiguration::write_formatted(string file_name) {
}
fclose(output);
}
+
+bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
+ if (number_of_spheres != other.number_of_spheres) {
+ return false;
+ }
+ if (number_of_scales != other.number_of_scales) {
+ return false;
+ }
+ if (idfc != other.idfc) {
+ return false;
+ }
+ if (exdc != other.exdc) {
+ return false;
+ }
+ if (wp != other.wp) {
+ return false;
+ }
+ if (xip != other.xip) {
+ return false;
+ }
+ //if (use_external_sphere != other.use_external_sphere) return false;
+ for (int svi = 0; svi < number_of_scales; svi++)
+ if (scale_vec[svi] != other.scale_vec[svi]) {
+ return false;
+ }
+ int dj_index = 0;
+ for (int dj = 0; dj < number_of_spheres; dj++) {
+ bool check_matrixes = false;
+ if (iog_vec[dj] >= dj + 1) {
+ dj_index = iog_vec[dj] - 1;
+ check_matrixes = true;
+ }
+ if (iog_vec[dj] != other.iog_vec[dj]) {
+ return false;
+ }
+ if (radii_of_spheres[dj_index] != other.radii_of_spheres[dj_index]) {
+ return false;
+ }
+ int layers = nshl_vec[dj_index];
+ if (layers != other.nshl_vec[dj_index]) {
+ return false;
+ }
+ if (check_matrixes) {
+ for (int di = 0; di < layers; di++) {
+ if (rcf[dj_index][di] != other.rcf[dj_index][di]) {
+ return false;
+ }
+ for (int dk = 0; dk < number_of_scales; dk++) {
+ if (dc0_matrix[di][dj_index][dk] != other.dc0_matrix[di][dj_index][dk]) {
+ return false;
+ }
+ } // dk loop
+ } // di loop
+ }
+ } // dj loop
+ return true;
+}
diff --git a/src/libnptm/Parsers.cpp b/src/libnptm/Parsers.cpp
index a327df810c95b58e3c0e94cfe207f104d482a5bd..3f8c9ff5116126b909baa9b7755f9bad09659642 100644
--- a/src/libnptm/Parsers.cpp
+++ b/src/libnptm/Parsers.cpp
@@ -1,4 +1,6 @@
/*! \file Parsers.cpp
+ *
+ * \brief Implementation of the parsing functions.
*/
#include <fstream>
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..688607811c261e46142354e5a3f0f9f719529c35
--- /dev/null
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -0,0 +1,330 @@
+/*! \file TransitionMatrix.cpp
+ *
+ * \brief Implementation of the Transition Matrix structure.
+ */
+#include <complex>
+#include <exception>
+#include <fstream>
+#include <hdf5.h>
+
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#include "../include/TransitionMatrix.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
+using namespace std;
+
+TransitionMatrix::~TransitionMatrix() {
+ if (elements != NULL) delete[] elements;
+ if (shape != NULL) delete[] shape;
+}
+
+TransitionMatrix::TransitionMatrix(
+ int _is, int _lm, double _vk, double _exri, std::complex<double> *_elements,
+ double _radius
+) {
+ is = _is;
+ l_max = _lm;
+ vk = _vk;
+ exri = _exri;
+ elements = _elements;
+ sphere_radius = _radius;
+ shape = new int[2]();
+ if (is == 1111) {
+ shape[0] = l_max;
+ shape[1] = 2;
+ } else if (is == 1) {
+ const int nlemt = 2 * l_max * (l_max + 2);
+ shape[0] = nlemt;
+ shape[1] = nlemt;
+ }
+}
+
+TransitionMatrix::TransitionMatrix(
+ int _lm, double _vk, double _exri, std::complex<double> **_rmi,
+ std::complex<double> **_rei, double _sphere_radius
+) {
+ is = 1111;
+ shape = new int[2];
+ shape[0] = _lm;
+ shape[1] = 2;
+ l_max = _lm;
+ vk = _vk;
+ exri = _exri;
+ sphere_radius = _sphere_radius;
+ elements = new complex<double>[2 * l_max]();
+ for (int ei = 0; ei < l_max; ei++) {
+ elements[2 * ei] = -1.0 / _rmi[ei][0];
+ elements[2 * ei + 1] = -1.0 / _rei[ei][0];
+ }
+}
+
+TransitionMatrix::TransitionMatrix(
+ int _nlemt, int _lm, double _vk, double _exri,
+ std::complex<double> **_am0m
+) {
+ is = 1;
+ shape = new int[2];
+ shape[0] = _nlemt;
+ shape[1] = _nlemt;
+ l_max = _lm;
+ vk = _vk;
+ exri = _exri;
+ sphere_radius = 0.0;
+ elements = new complex<double>[_nlemt * _nlemt]();
+ for (int ei = 0; ei < _nlemt; ei++) {
+ for (int ej = 0; ej < _nlemt; ej++) elements[_nlemt * ei + ej] = _am0m[ei][ej];
+ }
+}
+
+TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
+ TransitionMatrix *tm = NULL;
+ if (mode.compare("LEGACY") == 0) {
+ tm = TransitionMatrix::from_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ tm = TransitionMatrix::from_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+ return tm;
+}
+
+TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
+ TransitionMatrix *tm = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ if (status == 0) {
+ int _is;
+ int _lm;
+ double _vk;
+ double _exri;
+ // This vector will be passed to the new object. DO NOT DELETE HERE!
+ complex<double> *_elements;
+ double _radius = 0.0;
+ status = hdf_file->read("IS", "INT32", &_is);
+ status = hdf_file->read("L_MAX", "INT32", &_lm);
+ status = hdf_file->read("VK", "FLOAT64", &_vk);
+ status = hdf_file->read("EXRI", "FLOAT64", &_exri);
+ if (_is == 1111) {
+ int num_elements = 2 * _lm;
+ double *file_vector = new double[2 * num_elements]();
+ hid_t file_id = hdf_file->get_file_id();
+ hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
+ status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
+ _elements = new complex<double>[num_elements]();
+ for (int ei = 0; ei < num_elements; ei++) {
+ _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ }
+ status = H5Dclose(dset_id);
+ status = hdf_file->read("RADIUS", "FLOAT64", &_radius);
+ tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);
+ delete[] file_vector;
+ } else if (_is == 1) {
+ int nlemt = 2 * _lm * (_lm + 2);
+ int num_elements = nlemt * nlemt;
+ double *file_vector = new double[2 * num_elements]();
+ hid_t file_id = hdf_file->get_file_id();
+ hid_t dset_id = H5Dopen2(file_id, "ELEMENTS", H5P_DEFAULT);
+ status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, file_vector);
+ _elements = new complex<double>[num_elements]();
+ for (int ei = 0; ei < num_elements; ei++) {
+ _elements[ei] = complex<double>(file_vector[2 * ei], file_vector[2 * ei + 1]);
+ }
+ status = H5Dclose(dset_id);
+ tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);
+ delete[] file_vector;
+ }
+ status = hdf_file->close();
+ } else {
+ printf("ERROR: could not open file \"%s\"\n", file_name.c_str());
+ }
+ return tm;
+}
+
+TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
+ fstream ttms;
+ TransitionMatrix *tm = NULL;
+ ttms.open(file_name, ios::binary | ios::in);
+ if (ttms.is_open()) {
+ int num_elements = 0;
+ int _is;
+ int _lm;
+ double _vk;
+ double _exri;
+ // This vector will be passed to the new object. DO NOT DELETE HERE!
+ complex<double> *_elements;
+ double _radius = 0.0;
+ ttms.read(reinterpret_cast<char *>(&_is), sizeof(int));
+ ttms.read(reinterpret_cast<char *>(&_lm), sizeof(int));
+ ttms.read(reinterpret_cast<char *>(&_vk), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&_exri), sizeof(double));
+ if (_is == 1111) {
+ num_elements = _lm * 2;
+ _elements = new complex<double>[num_elements]();
+ for (int ei = 0; ei < num_elements; ei++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ _elements[ei] = complex<double>(vreal, vimag);
+ }
+ double _radius;
+ ttms.read(reinterpret_cast<char *>(&_radius), sizeof(double));
+ tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements, _radius);
+ } else if (_is == 1) {
+ int nlemt = 2 * _lm * (_lm + 2);
+ num_elements = nlemt * nlemt;
+ _elements = new complex<double>[num_elements]();
+ for (int ei = 0; ei < num_elements; ei++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ _elements[ei] = complex<double>(vreal, vimag);
+ }
+ tm = new TransitionMatrix(_is, _lm, _vk, _exri, _elements);
+ }
+ } else {
+ printf("ERROR: could not open file \"%s\"\n", file_name.c_str());
+ }
+ return tm;
+}
+
+void TransitionMatrix::write_binary(string file_name, string mode) {
+ if (mode.compare("LEGACY") == 0) {
+ write_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ write_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+}
+
+void TransitionMatrix::write_hdf5(string file_name) {
+ if (is == 1 || is == 1111) {
+ List<string> rec_name_list(1);
+ List<string> rec_type_list(1);
+ List<void *> rec_ptr_list(1);
+ string str_type, str_name;
+ rec_name_list.set(0, "IS");
+ rec_type_list.set(0, "INT32_(1)");
+ rec_ptr_list.set(0, &is);
+ rec_name_list.append("L_MAX");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&l_max);
+ rec_name_list.append("VK");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&vk);
+ rec_name_list.append("EXRI");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&exri);
+ rec_name_list.append("ELEMENTS");
+ str_type = "FLOAT64_(" + to_string(shape[0]) + "," + to_string(2 * shape[1]) + ")";
+ rec_type_list.append(str_type);
+ // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,
+ // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)
+ int num_elements = 2 * shape[0] * shape[1];
+ double *ptr_elements = new double[num_elements]();
+ for (int ei = 0; ei < num_elements / 2; ei++) {
+ ptr_elements[2 * ei] = elements[ei].real();
+ ptr_elements[2 * ei + 1] = elements[ei].imag();
+ }
+ rec_ptr_list.append(ptr_elements);
+ if (is == 1111) {
+ rec_name_list.append("RADIUS");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&sphere_radius);
+ }
+
+ string *rec_names = rec_name_list.to_array();
+ string *rec_types = rec_type_list.to_array();
+ void **rec_pointers = rec_ptr_list.to_array();
+ const int rec_num = rec_name_list.length();
+ FileSchema schema(rec_num, rec_types, rec_names);
+ HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
+ for (int ri = 0; ri < rec_num; ri++)
+ hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ hdf_file->close();
+
+ delete[] ptr_elements;
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete hdf_file;
+ } else {
+ string message = "Unrecognized matrix data.";
+ throw UnrecognizedFormatException(message);
+ }
+}
+
+void TransitionMatrix::write_legacy(string file_name) {
+ fstream ttms;
+ if (is == 1111 || is == 1) {
+ ttms.open(file_name, ios::binary | ios::out);
+ if (ttms.is_open()) {
+ ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&l_max), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ }
+ } else {
+ string message = "Unrecognized matrix data.";
+ throw UnrecognizedFormatException(message);
+ }
+ if (ttms.is_open()) {
+ int num_elements = shape[0] * shape[1];
+ for (int ei = 0; ei < num_elements; ei++) {
+ complex<double> element1 = elements[ei];
+ double vreal, vimag;
+ vreal = element1.real();
+ vimag = element1.imag();
+ ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+ }
+ if (is == 1111) {
+ ttms.write(reinterpret_cast<char *>(&sphere_radius), sizeof(double));
+ }
+ ttms.close();
+ } else { // Failed to open output file. Should never happen.
+ printf("ERROR: could not open Transition Matrix file for writing.\n");
+ }
+}
+
+bool TransitionMatrix::operator ==(TransitionMatrix &other) {
+ if (is != other.is) {
+ return false;
+ }
+ if (l_max != other.l_max) {
+ return false;
+ }
+ if (vk != other.vk) {
+ return false;
+ }
+ if (exri != other.exri) {
+ return false;
+ }
+ if (sphere_radius != other.sphere_radius) {
+ return false;
+ }
+ if (shape[0] != other.shape[0]) {
+ return false;
+ }
+ if (shape[1] != other.shape[1]) {
+ return false;
+ }
+ int num_elements = shape[0] * shape[1];
+ for (int ei = 0; ei < num_elements; ei++) {
+ if (elements[ei] != other.elements[ei]) {
+ return false;
+ }
+ }
+ return true;
+}
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index b516e576ce250e6259d33e7b8141c1efb19270dd..22287bb7c734fcfc067c4e185b4b4f7932c03918 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -19,8 +19,8 @@
using namespace std;
void apc(
- double ****zpv, int le, complex<double> **am0m, complex<double> **w,
- double sqk, double **gapr, complex<double> **gapp
+ double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
+ double sqk, double **gapr, std::complex<double> **gapp
) {
complex<double> **ac, **gap;
const complex<double> cc0(0.0, 0.0);
@@ -137,8 +137,8 @@ void apc(
}
void apcra(
- double ****zpv, const int le, complex<double> **am0m, int inpol, double sqk,
- double **gaprm, complex<double> **gappm
+ double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
+ double **gaprm, std::complex<double> **gappm
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
@@ -358,7 +358,7 @@ void apcra(
}
complex<double> cdtp(
- complex<double> z, complex<double> **am, int i, int jf,
+ std::complex<double> z, std::complex<double> **am, int i, int jf,
int k, int nj
) {
/* NOTE: the original FORTRAN code treats the AM matrix as a
@@ -409,7 +409,7 @@ double cgev(int ipamo, int mu, int l, int m) {
return result;
}
-void cms(complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
+void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
complex<double> dm, de, cgh, cgk;
const complex<double> cc0(0.0, 0.0);
int ndi = c4->nsph * c4->nlim;
@@ -626,7 +626,7 @@ complex<double> ghit(
if (ihi == 2) {
if (c1->rxx[nbl - 1] == 0.0 && c1->ryy[nbl - 1] == 0.0 && c1->rzz[nbl - 1] == 0.0) {
if (ipamo == 0) {
- if (l1 == l2 && m1 == m2) result = complex(1.0, 0.0);
+ if (l1 == l2 && m1 == m2) result = complex<double>(1.0, 0.0);
}
return result;
}
@@ -823,7 +823,7 @@ complex<double> ghit(
}
void hjv(
- double exri, double vk, int &jer, int &lcalc, complex<double> &arg,
+ double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
) {
int nsphmo = c4->nsph - 1;
@@ -892,7 +892,7 @@ void hjv(
delete[] rfn;
}
-void lucin(complex<double> **am, const int nddmst, int n, int &ier) {
+void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
/* NDDMST FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
* STATEMENT.
* N NUMBER OF ROWS IN AM.
@@ -997,7 +997,7 @@ void lucin(complex<double> **am, const int nddmst, int n, int &ier) {
delete[] v;
}
-void mextc(double vk, double exri, complex<double> **fsac, double **cextlr, double **cext) {
+void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext) {
double fa11r = fsac[0][0].real();
double fa11i = fsac[0][0].imag();
double fa21r = fsac[1][0].real();
@@ -1538,8 +1538,8 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
}
void raba(
- int le, complex<double> **am0m, complex<double> **w, double **tqce,
- complex<double> **tqcpe, double **tqcs, complex<double> **tqcps
+ int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
+ std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
) {
complex<double> **a, **ctqce, **ctqcs;
complex<double> acw, acwp, aca, acap, c1, c2, c3;
@@ -1901,7 +1901,7 @@ void tqr(
tsk = u[0] * tqsv[0] + u[1] * tqsv[1] + u[2] * tqsv[2];
}
-void ztm(complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
+void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
const complex<double> cc0(0.0, 0.0);
int ndi = c4->nsph * c4->nlim;
diff --git a/src/libnptm/file_io.cpp b/src/libnptm/file_io.cpp
index a5df8a6ab6ae42f7d758830dfad5954518991fca..0c0b0c170009bb40abf290564ee838b68f1260a9 100644
--- a/src/libnptm/file_io.cpp
+++ b/src/libnptm/file_io.cpp
@@ -47,7 +47,10 @@ string* FileSchema::get_record_types() {
HDFFile::HDFFile(string name, unsigned int flags, hid_t fcpl_id, hid_t fapl_id) {
file_name = name;
- file_id = H5Fcreate(name.c_str(), flags, fcpl_id, fapl_id);
+ if (flags == H5F_ACC_EXCL || flags == H5F_ACC_TRUNC)
+ file_id = H5Fcreate(name.c_str(), flags, fcpl_id, fapl_id);
+ else if (flags == H5F_ACC_RDONLY || flags == H5F_ACC_RDWR)
+ file_id = H5Fopen(name.c_str(), flags, fapl_id);
id_list = new List<hid_t>(1);
id_list->set(0, file_id);
if (file_id != H5I_INVALID_HID) file_open_flag = true;
@@ -132,6 +135,45 @@ HDFFile* HDFFile::from_schema(
return hdf_file;
}
+herr_t HDFFile::read(
+ string dataset_name, string data_type, void *buffer,
+ hid_t mem_space_id, hid_t file_space_id, hid_t dapl_id,
+ hid_t dxpl_id
+) {
+ string known_types[] = {"INT32", "FLOAT64"};
+ regex re;
+ smatch m;
+ bool found_type = false;
+ int type_index = 0;
+ while (!found_type) {
+ re = regex(known_types[type_index++]);
+ found_type = regex_search(data_type, m, re);
+ if (type_index == 2) break;
+ }
+ if (found_type) {
+ hid_t dataset_id = H5Dopen2(file_id, dataset_name.c_str(), dapl_id);
+ hid_t mem_type_id;
+ switch (type_index) {
+ case 1:
+ mem_type_id = H5T_NATIVE_INT; break;
+ case 2:
+ mem_type_id = H5T_NATIVE_DOUBLE; break;
+ default:
+ throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ }
+ if (dataset_id != H5I_INVALID_HID) {
+ status = H5Dread(dataset_id, mem_type_id, mem_space_id, file_space_id, dxpl_id, buffer);
+ if (status == 0) status = H5Dclose(dataset_id);
+ else status = (herr_t)-2;
+ } else {
+ status = (herr_t)-1;
+ }
+ } else {
+ throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ }
+ return status;
+}
+
herr_t HDFFile::write(
string dataset_name, string data_type, const void *buffer,
hid_t mem_space_id, hid_t file_space_id, hid_t dapl_id,
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 99aa4d7f25981421648ca34267e322db332c59aa..30e31bf9ec78e83a8e545a408d93c22fb6a267af 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -64,7 +64,7 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
}
}
-void cbf(int n, complex<double> z, int &nm, complex<double> csj[]) {
+void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
/*
* FROM CSPHJY OF LIBRARY specfun
*
@@ -177,7 +177,7 @@ double cg1(int lmpml, int mu, int l, int m) {
return result;
}
-complex<double> dconjg(complex<double> z) {
+complex<double> dconjg(std::complex<double> z) {
double zreal = z.real();
double zimag = z.imag();
return complex<double>(zreal, -zimag);
@@ -204,7 +204,7 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- C1 *c1, C2 *c2, int &jer, int &lcalc, complex<double> &arg
+ C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
) {
const int lipo = li + 1;
const int lipt = li + 2;
@@ -353,7 +353,7 @@ double envj(int n, double x) {
return result;
}
-void mmulc(complex<double> *vint, double **cmullr, double **cmul) {
+void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
double sm2 = vint[0].real();
double s24 = vint[1].real();
double d24 = vint[1].imag();
@@ -481,7 +481,7 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
}
void pwma(
- double *up, double *un, complex<double> *ylm, int inpol, int lw,
+ double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
int isq, C1 *c1
) {
const double four_pi = 8.0 * acos(0.0);
@@ -560,8 +560,8 @@ void pwma(
}
void rabas(
- int inpol, int li, int nsph, C1 *c1, double **tqse, complex<double> **tqspe,
- double **tqss, complex<double> **tqsps
+ int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
+ double **tqss, std::complex<double> **tqsps
) {
complex<double> cc0 = complex<double>(0.0, 0.0);
complex<double> uim = complex<double>(0.0, 1.0);
@@ -672,9 +672,9 @@ void rbf(int n, double x, int &nm, double sj[]) {
}
void rkc(
- int npntmo, double step, complex<double> dcc, double &x, int lpo,
- complex<double> &y1, complex<double> &y2, complex<double> &dy1,
- complex<double> &dy2
+ int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
+ std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
+ std::complex<double> &dy2
) {
complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
complex<double> cy4, yy, c14, c21, c22, c23, c24;
@@ -710,8 +710,8 @@ void rkc(
}
void rkt(
- int npntmo, double step, double &x, int lpo, complex<double> &y1,
- complex<double> &y2, complex<double> &dy1, complex<double> &dy2,
+ int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
+ std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
C2 *c2
) {
complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
@@ -804,7 +804,7 @@ void rnf(int n, double x, int &nm, double sy[]) {
void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, complex<double> *ylm
+ int ll, std::complex<double> *ylm
) {
const int rmp_size = ll;
const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
@@ -881,7 +881,7 @@ void sphar(
goto label40;
}
-void sscr0(complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
+void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
complex<double> sum21, rm, re, csam;
const complex<double> cc0 = complex<double>(0.0, 0.0);
const double exdc = exri * exri;
@@ -923,7 +923,7 @@ void sscr0(complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1
}
void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
- complex<double> s11, s21, s12, s22, rm, re, csam, cc;
+ std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
const complex<double> cc0(0.0, 0.0);
double ccs = 1.0 / (vk * vk);
csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index 792ed3feaaa38f22b25c727d3f52def65a7da2ca..aef2de470a5c3fe0b68aa95f05c0dbdac9d63615 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -21,7 +21,7 @@
using namespace std;
void camp(
- complex<double> *ac, complex<double> **am0m, complex<double> *ws,
+ std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
CIL *cil
) {
for (int j = 0; j < cil->nlemt; j++) {
@@ -32,8 +32,8 @@ void camp(
}
void czamp(
- complex<double> *ac, complex<double> **amd, int **indam,
- complex<double> *ws, CIL *cil
+ std::complex<double> *ac, std::complex<double> **amd, int **indam,
+ std::complex<double> *ws, CIL *cil
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
@@ -61,7 +61,7 @@ void czamp(
}
void ffrf(
- double ****zpv, complex<double> *ac, complex<double> *ws, double *fffe,
+ double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
) {
const complex<double> cc0(0.0, 0.0);
@@ -178,7 +178,7 @@ void ffrf(
}
void ffrt(
- complex<double> *ac, complex<double> *ws, double *ffte, double *ffts,
+ std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
CIL *cil
) {
const complex<double> cc0(0.0, 0.0);
@@ -282,7 +282,7 @@ void frfmer(
delete[] wk;
}
-void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, int lw) {
+void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
complex<double> cp1, cm1, cp2, cm2, cl;
const complex<double> uim(0.0, 1.0);
const double four_pi = acos(0.0) * 8.0;
@@ -320,8 +320,8 @@ void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, i
}
void samp(
- complex<double> *ac, complex<double> *tmsm, complex<double> *tmse,
- complex<double> *ws, CIL *cil
+ std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
+ std::complex<double> *ws, CIL *cil
) {
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
@@ -337,7 +337,7 @@ void samp(
}
void sampoa(
- complex<double> *ac, complex<double> **tms, complex<double> *ws,
+ std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
CIL *cil
) {
complex<double> **tm = new complex<double>*[2];
@@ -364,7 +364,7 @@ void sampoa(
}
void wamff(
- complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);
diff --git a/src/scripts/pycompare.py b/src/scripts/pycompare.py
index b253da6682e141465dfa56690cf8b08ca97a5071..3720701751ae34a02a30326acf7473cf00b85721 100755
--- a/src/scripts/pycompare.py
+++ b/src/scripts/pycompare.py
@@ -166,6 +166,7 @@ def compare_lines(f_line, c_line, config, line_num=0, num_len=4, log_file=None):
f_line = f_line.replace("D-","E-").replace("D+","E+")
ref_format = " <div><span style=\"font-weight: bold; color: rgb(125,125,125)\"><pre><code>{0:%ds}"%num_len
ref_line = (ref_format + ": {1:s}</code></pre></span></div>\n").format("ORIG", f_line[:-1])
+ log_line = ""
if (f_line == c_line):
if log_file is not None:
if (config['full_log']):
@@ -189,9 +190,10 @@ def compare_lines(f_line, c_line, config, line_num=0, num_len=4, log_file=None):
if (len(f_groups) == len(c_groups)):
severities = mismatch_severities(f_groups, c_groups, config)
if log_file is not None:
- num_format = " <div><pre><code>{0:0%dd}"%num_len
- log_line = (num_format + ": ").format(line_num)
- log_line = log_line + c_line[0:c_starts[0]]
+ if (len(severities) > 0):
+ num_format = " <div><pre><code>{0:0%dd}"%num_len
+ log_line = (num_format + ": ").format(line_num)
+ log_line = log_line + c_line[0:c_starts[0]]
for si in range(len(severities) - 1):
if (severities[si] == 1): noisy += 1
elif (severities[si] == 2): warnings += 1
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index 305ee34a3ba3abceece74146a12cf0f412fbdfdf..d012622e1d346bceca21f3bcf3f819b8614b80b7 100644
--- a/src/sphere/Makefile
+++ b/src/sphere/Makefile
@@ -13,8 +13,7 @@ include ../make.inc
F_SPH_OBJS=$(OBJDIR)/edfb_sph.o $(OBJDIR)/sph.o
-CXX_SPH_OBJS=$(OBJDIR)/np_sphere.o $(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/sphere.o
-
+CXX_SPH_OBJS=$(OBJDIR)/np_sphere.o $(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/sphere.o $(OBJDIR)/TransitionMatrix.o
all: $(BUILDDIR_SPH)/edfb_sph $(BUILDDIR_SPH)/sph $(BUILDDIR_SPH)/np_sphere
diff --git a/src/sphere/file_io.f b/src/sphere/file_io.f
deleted file mode 100644
index ebd708c5e50fb58d0dec7962cd8ff5391c19e0ef..0000000000000000000000000000000000000000
--- a/src/sphere/file_io.f
+++ /dev/null
@@ -1,33 +0,0 @@
- 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
diff --git a/src/sphere/np_sphere.cpp b/src/sphere/np_sphere.cpp
index 1846b36973859d1a175d2ff7f579ed1de064b6fe..e3033dae6fa36f373e17c1fda92f802537ce3b53 100644
--- a/src/sphere/np_sphere.cpp
+++ b/src/sphere/np_sphere.cpp
@@ -1,4 +1,19 @@
/*! \file np_sphere.cpp
+ *
+ * \brief Single sphere scattering problem handler.
+ *
+ * This program emulates the execution work-flow originally handled by the
+ * FORTRAN EDFB and SPH codes, which undertook the task of solving the
+ * scattering calculation for the case of a single, possibly multi-layered
+ * sphere. The program is designed to work in two modes: emulation and
+ * enhanced. The emulation mode is activated by invoking the program without
+ * arguments. In this case, the code looks for hard-coded default input and
+ * writes output in the execution folder, replicating the behaviour of the
+ * original FORTRAN code. Advanced mode, instead, is activated by passing
+ * command line arguments that locate the desired input files and a valid
+ * folder to write the output into. The emulation mode is useful for testing,
+ * while the advanced mode implements the possibility to change input and
+ * output options, without having to modify the code.
*/
#include <complex>
#include <cstdio>
@@ -19,16 +34,20 @@ extern void sphere(string config_file, string data_file, string output_path);
* or some otherwise formatted configuration data set. The code can be
* linked to a luncher script or to a GUI oriented application that performs
* the configuration and runs the main program.
+ *
+ * \param argc: `int` The number of arguments given in command-line.
+ * \param argv: `char **` The vector of command-line arguments.
+ * \return result: `int` An exit code passed to the OS (0 for succesful execution).
*/
int main(int argc, char **argv) {
- string config_file = "../../test_data/sphere/DEDFB";
- string data_file = "../../test_data/sphere/DSPH";
- string output_path = ".";
- if (argc == 4) {
- config_file = string(argv[1]);
- data_file = string(argv[2]);
- output_path = string(argv[3]);
- }
- sphere(config_file, data_file, output_path);
- return 0;
+ string config_file = "../../test_data/sphere/DEDFB";
+ string data_file = "../../test_data/sphere/DSPH";
+ string output_path = ".";
+ if (argc == 4) {
+ config_file = string(argv[1]);
+ data_file = string(argv[2]);
+ output_path = string(argv[3]);
+ }
+ sphere(config_file, data_file, output_path);
+ return 0;
}
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 2e179fc59ecb5aaff9e23f68a6a2494757a632e8..efad7ad1819f879ab0128b2ed08e10ca7fd4b33f 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -1,9 +1,12 @@
/*! \file sphere.cpp
+ *
+ * \brief Implementation of the single sphere calculation.
*/
+#include <complex>
#include <cstdio>
+#include <exception>
#include <fstream>
#include <string>
-#include <complex>
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
@@ -17,6 +20,10 @@
#include "../include/sph_subs.h"
#endif
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#include "../include/TransitionMatrix.h"
+#endif
+
using namespace std;
/*! \brief C++ implementation of SPH
@@ -39,6 +46,7 @@ void sphere(string config_file, string data_file, string output_path) {
}
sconf->write_formatted(output_path + "/c_OEDFB");
sconf->write_binary(output_path + "/c_TEDF");
+ sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");
GeometryConfiguration *gconf = NULL;
try {
gconf = GeometryConfiguration::from_legacy(data_file);
@@ -216,6 +224,7 @@ void sphere(string config_file, string data_file, string output_path) {
tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nsph), sizeof(int));
for (int nsi = 0; nsi < nsph; nsi++)
tppoan.write(reinterpret_cast<char *>(&(sconf->iog_vec[nsi])), sizeof(int));
@@ -284,39 +293,14 @@ void sphere(string config_file, string data_file, string output_path) {
} // i132
if (sconf->idfc >= 0 and nsph == 1 and jxi == gconf->jwtm) {
// This is the condition that writes the transition matrix to output.
- int is = 1111;
- fstream ttms;
- string ttms_name = output_path + "/c_TTMS";
- ttms.open(ttms_name.c_str(), ios::binary | ios::out);
- if (ttms.is_open()) {
- ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
- ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
- ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
- for (int lmi = 0; lmi < gconf->l_max; lmi++) {
- complex<double> element1 = -1.0 / c1->rmi[lmi][0];
- complex<double> element2 = -1.0 / c1->rei[lmi][0];
- double vreal, vimag;
- vreal = element1.real();
- vimag = element1.imag();
- ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- vreal = element2.real();
- vimag = element2.imag();
- ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
- ttms.close();
- } else { // Failed to open output file. Should never happen.
- printf("ERROR: could not open TTMS file.\n");
- tppoan.close();
- }
+ TransitionMatrix ttms(gconf->l_max, vk, exri, c1->rmi, c1->rei, sconf->radii_of_spheres[0]);
+ string ttms_name = output_path + "/c_TTMS.hd5";
+ ttms.write_binary(ttms_name, "HDF5");
+ ttms_name = output_path + "/c_TTMS";
+ ttms.write_binary(ttms_name);
}
double cs0 = 0.25 * vk * vk * vk / half_pi;
- //printf("DEBUG: cs0 = %lE\n", cs0);
sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
- //printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
@@ -333,7 +317,7 @@ void sphere(string config_file, string data_file, string output_path) {
} else {
fprintf(
output,
- " SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
+ " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
c2->vsz[i170],
c2->vkt[i170].real(),
c2->vkt[i170].imag()
diff --git a/src/testing/Makefile b/src/testing/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..de9a26d7f070decadead8901a7b32f0bb06c625b
--- /dev/null
+++ b/src/testing/Makefile
@@ -0,0 +1,27 @@
+ifndef BUILDDIR
+override BUILDDIR=../../build/testing
+endif
+ifndef OBJDIR
+override OBJDIR=../objects
+endif
+
+include ../make.inc
+
+
+CXX_TEDF_OBJS=$(OBJDIR)/test_TEDF.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o
+
+CXX_TTMS_OBJS=$(OBJDIR)/test_TTMS.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/TransitionMatrix.o
+
+all: $(BUILDDIR)/test_TEDF $(BUILDDIR)/test_TTMS
+
+$(BUILDDIR)/test_TEDF: $(CXX_TEDF_OBJS)
+ $(CXX) $(CXXFLAGS) -o $(BUILDDIR)/test_TEDF $(CXX_TEDF_OBJS) $(CXXLDFLAGS)
+
+$(BUILDDIR)/test_TTMS: $(CXX_TTMS_OBJS)
+ $(CXX) $(CXXFLAGS) -o $(BUILDDIR)/test_TTMS $(CXX_TTMS_OBJS) $(CXXLDFLAGS)
+
+clean:
+ rm -f $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS)
+
+wipe:
+ rm -f $(BUILDDIR)/test_TEDF $(BUILDDIR)/test_TTMS $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS)
diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..037a7351505c809ca8d2609e81a68ef9b91baf71
--- /dev/null
+++ b/src/testing/test_TEDF.cpp
@@ -0,0 +1,63 @@
+//! \file test_TEDF.cpp
+
+#include <complex>
+#include <cstdio>
+#include <hdf5.h>
+#include <string>
+#include "../include/List.h"
+#include "../include/file_io.h"
+#include "../include/Configuration.h"
+
+using namespace std;
+
+/*! \brief Main program execution body.
+ *
+ * This program executes a test to compare whether three configuration
+ * instances, loaded respectively from an EDFB configuration file, a
+ * legacy binary, or a HDF5 binary are actually equivalent. The test
+ * writes a result message to `stdout` then it returns 0 (OS flag for
+ * a successful process) or some kind of error code, depending on
+ * whether the test files were found all equal or not. The test accepts
+ * three command line arguments: the name of the EDFB configuration
+ * file, the name of the legacy binary file and the name of the HDF5
+ * binary configuration file.
+ *
+ * \param argc: `int` Number of command line arguments
+ * \param argv: `char **` Array of command argument character strings.
+ * \return result: `int` Can be: 0 (all files equal); 1 (EDFB and
+ * legacy binary are different); 10 (EDFB and HDF5 are different);
+ * 100 (legacy and HDF5 are different). In case more differences are
+ * found, the error codes sum up together (e.g. 111 means all files
+ * are different).
+ */
+int main(int argc, char **argv) {
+ int result = 0;
+ string dedfb_file = "DEDFB";
+ string legacy_file = "c_TEDF";
+ string hdf5_file = "c_TEDF.hd5";
+ if (argc == 4) {
+ dedfb_file = string(argv[1]);
+ legacy_file = string(argv[2]);
+ hdf5_file = string(argv[3]);
+ }
+ ScattererConfiguration *a, *b, *c;
+ a = ScattererConfiguration::from_dedfb(dedfb_file);
+ b = ScattererConfiguration::from_binary(legacy_file);
+ c = ScattererConfiguration::from_binary(hdf5_file, "HDF5");
+ if (*a == *b) printf("Configuration objects a and b are equal.\n");
+ else {
+ printf("Configuration objects a and b are different.\n");
+ result += 1;
+ }
+ if (*a == *c) printf("Configuration objects a and c are equal.\n");
+ else {
+ printf("Configuration objects a and c are different.\n");
+ result += 10;
+ }
+ if (*c == *b) printf("Configuration objects c and b are equal.\n");
+ else {
+ printf("Configuration objects c and b are different.\n");
+ result += 100;
+ }
+ return 0;
+}
diff --git a/src/testing/test_TTMS.cpp b/src/testing/test_TTMS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bdf0daaedc510a210ac1e314f69cf93b91a36f88
--- /dev/null
+++ b/src/testing/test_TTMS.cpp
@@ -0,0 +1,44 @@
+//! \file test_TTMS.cpp
+
+#include <complex>
+#include <cstdio>
+#include <hdf5.h>
+#include <string>
+#include "../include/List.h"
+#include "../include/file_io.h"
+#include "../include/TransitionMatrix.h"
+
+using namespace std;
+
+/*! \brief Main program execution body.
+ *
+ * This program executes a test to compare whether two transition
+ * matrix instances, loaded respectively from a legacy and a HDF5
+ * binary file are actually equivalent. The test writes a result
+ * message to `stdout` then it returns 0 (OS flag for a successful
+ * process) or 1 (OS flag for failing process) depending on whether
+ * the two instances were considered equivalent or not.
+ *
+ * \param argc: `int` Number of command line arguments
+ * \param argv: `char **` Array of command argument character strings.
+ * \return result: `int` Can be 0 (files are equal) or 1 (files are
+ * different).
+ */
+int main(int argc, char **argv) {
+ int result = 0;
+ TransitionMatrix *a, *b;
+ string legacy_file = "c_TTMS";
+ string hdf5_file = "c_TTMS.hd5";
+ if (argc == 3) {
+ legacy_file = string(argv[1]);
+ hdf5_file = string(argv[2]);
+ }
+ a = TransitionMatrix::from_binary(legacy_file);
+ b = TransitionMatrix::from_binary(hdf5_file, "HDF5");
+ if (*a == *b) printf("Transition matrixes a and b are equal.\n");
+ else {
+ printf("Transition matrixes a and b are different.\n");
+ result = 1;
+ }
+ return result;
+}
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index fd7c342083a557c2718ed41c7c3b617593649f75..f5d7af1cfefb7a5a1d41591962acbc9fbf3a0ae4 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -1,4 +1,6 @@
-/*! \file frfme.cpp
+/*! \file cfrfme.cpp
+ *
+ * C++ implementation of FRFME functions.
*/
#include <complex>
#include <cstdio>
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index e9c1f42147d42c4db1b1cac7ee978649bf582337..9cd9b50b409b48383fb02d74983ef761a5c71506 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -1,4 +1,6 @@
-/*! \file lffft.cpp
+/*! \file clffft.cpp
+ *
+ * \brief C++ implementation of LFFFT functions.
*/
#include <complex>
#include <cstdio>
diff --git a/src/trapping/np_trapping.cpp b/src/trapping/np_trapping.cpp
index cb3d9be08dbdbb8d0e374039d6fd7ba90f336b5c..99b1f7898e900cd852113f427b926152cce08c91 100644
--- a/src/trapping/np_trapping.cpp
+++ b/src/trapping/np_trapping.cpp
@@ -1,4 +1,6 @@
-/*! \file trapping.cpp
+/*! \file np_trapping.cpp
+ *
+ * \brief Trapping problem handler.
*/
#include <cstdio>
#include <string>
@@ -8,7 +10,7 @@ using namespace std;
extern void frfme(string data_file, string output_path);
extern void lffft(string data_file, string output_path);
-/* \brief Main program execution body.
+/*! \brief Main program execution body.
*
* \param argc: `int`
* \param argv: `char **`