diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index e38967644a90d8377c0eb28895aca283feec760c..22bc637a66fc67ed87f40b4355315a0e34ba87c8 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -24,6 +24,8 @@ 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 $(OBJDIR)/TransitionMatrix.o
CXX_CLU_OBJS=$(OBJDIR)/np_cluster.o $(OBJDIR)/cluster.o $(LIBNPTM)
+CXX_CLU_DEBUG=$(OBJDIR)/np_cluster.g* $(OBJDIR)/cluster.g*
+
all: $(LIBNPTM) $(BUILDDIR_CLU)/clu $(BUILDDIR_CLU)/edfb_clu $(BUILDDIR_CLU)/np_cluster
@@ -47,8 +49,8 @@ $(BUILDDIR_CLU)/np_cluster: $(OBJDIR) $(CXX_CLU_OBJS) $(BUILDDIR_CLU) $(LIBNPTM)
$(CXX) $(CXXFLAGS) -o $(BUILDDIR_CLU)/np_cluster $(CXX_CLU_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
clean:
- rm -f $(F_CLU_OBJS) $(CXX_CLU_OBJS)
+ rm -f $(F_CLU_OBJS) $(CXX_CLU_OBJS) $(CXX_CLU_DEBUG)
wipe:
- rm -f $(BUILDDIR_CLU)/clu $(BUILDDIR_CLU)/edfb_clu $(BUILDDIR_CLU)/np_cluster $(F_CLU_OBJS) $(CXX_CLU_OBJS)
+ rm -f $(BUILDDIR_CLU)/clu $(BUILDDIR_CLU)/edfb_clu $(BUILDDIR_CLU)/np_cluster $(F_CLU_OBJS) $(CXX_CLU_OBJS) $(CXX_CLU_DEBUG)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 94bec3c35c601f5433b063b2350673b482e33625..4591b28aa6bda3f1483f0e0190af7bea94d934f2 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -8,6 +8,10 @@
#include <fstream>
#include <string>
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 37338bdf290661633c6e33d688dcb13008090f7c..f3c4df91bcc64e6673865b739b5f7dd9a06af5a3 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -28,52 +28,6 @@
#ifndef INCLUDE_CONFIGURATION_H_
#define INCLUDE_CONFIGURATION_H_
-/**
- * \brief Exception for open file error handlers.
- */
-class OpenConfigurationFileException: public std::exception {
-protected:
- //! \brief Name of the file that was accessed.
- std::string file_name;
-
-public:
- /**
- * \brief Exception instance constructor.
- *
- * \param name: `string` Name of the file that was accessed.
- */
- OpenConfigurationFileException(std::string name) { file_name = name; }
-
- /**
- * \brief Exception message.
- */
- virtual const char* what() const throw() {
- return file_name.c_str();
- }
-};
-
-/**
- * \brief Exception for unrecognized configuration data sets.
- */
-class UnrecognizedConfigurationException: 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.
- */
- UnrecognizedConfigurationException(std::string problem) { message = problem; }
- /**
- * \brief Exception message.
- */
- virtual const char* what() const throw() {
- return message.c_str();
- }
-};
-
/**
* \brief A class to represent the configuration of the scattering geometry.
*
@@ -361,6 +315,44 @@ public:
*/
static ScattererConfiguration* from_dedfb(std::string file_name);
+ /*! \brief Get the ID of a sphere by its index.
+ *
+ * The proper way to access read-only parameters from outside a class is to define
+ * public methods that return their values. For arrays, particularly those that
+ * are accessed multiple times, it is convenient to have specialized methods that
+ * return the required values based on their index in the array.
+ *
+ * \param index: `int` Index of the ID to be retrieved.
+ * \return id: `int` The desired identifier.
+ */
+ int get_iog(int index) { return iog_vec[index]; }
+
+ /*! \brief Get the value of a parameter by name.
+ *
+ * The proper way to access read-only parameters from outside a class is to define
+ * public methods that return their values. For configuration operations, whose
+ * optimization is not critical, it is possible to define a single function that
+ * returns simple scalar values called by name. Access to more complicated data
+ * structures, on the other hand, require specialized methods which avoid the
+ * burden of searching the necessary value across the whole arrya every time.
+ *
+ * \param param_name: `string` Name of the parameter to be retrieved.
+ * \return value: `double` Value of the requested parameter.
+ */
+ double get_param(std::string param_name);
+
+ /*! \brief Get the value of a scale by its index.
+ *
+ * The proper way to access read-only parameters from outside a class is to define
+ * public methods that return their values. For arrays, particularly those that
+ * are accessed multiple times, it is convenient to have specialized methods that
+ * return the required values based on their index in the array.
+ *
+ * \param index: `int` Index of the scale to be retrieved.
+ * \return scale: `double` The desired scale.
+ */
+ double get_scale(int index) { return scale_vec[index]; }
+
/*! \brief Print the contents of the configuration object to terminal.
*
* In case of quick debug testing, `ScattererConfiguration.print()` allows printing
diff --git a/src/include/List.h b/src/include/List.h
index 8c51bd1ce7bb5816a48bd958c2c6656216e5a0fa..f97ebc6004c0867ce552480c69e16a11d671bc07 100644
--- a/src/include/List.h
+++ b/src/include/List.h
@@ -6,38 +6,7 @@
#ifndef INCLUDE_LIST_H_
#define INCLUDE_LIST_H_
-/**
- * \brief Exception for out of bounds List requests.
- */
-class ListOutOfBoundsException: public std::exception {
-protected:
- //! Description of the problem.
- std::string message;
-
-public:
- /**
- * \brief Exception instance constructor.
- *
- * \param requested: `int` The index that was requested.
- * \param min: `int` The minimum index allowed by the list.
- * \param max: `int` The maximum index allowed by the list.
- */
- ListOutOfBoundsException(int requested, int min, int max) {
- message = "Error: requested index " + std::to_string(requested)
- + " out of list allowed bounds [" + std::to_string(min) + ", "
- + std::to_string(max - 1) + "]";
- }
-
- /**
- * \brief Exception message.
- */
- virtual const char* what() const throw() {
- return message.c_str();
- }
-};
-
-/**
- * \brief A class to represent dynamic lists.
+/*! \brief A class to represent dynamic lists.
*
* This class helps in the creation and management of dynamic lists of
* objects, whose size is not known in advance. List offers the advantage
diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
index 25e5cf86b0b5b3a801bca5db26e274c510f7e748..447457f43eb009556fcc422b75e1fd171f636ed0 100644
--- a/src/include/TransitionMatrix.h
+++ b/src/include/TransitionMatrix.h
@@ -6,28 +6,6 @@
#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 {
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 9bc2f72a9c50808d78d9f035fac89cf22446f217..49b89e20875919f916cea21a269fa7b210da3a9f 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -15,7 +15,11 @@
#ifndef INCLUDE_CLU_SUBS_H_
#define INCLUDE_CLU_SUBS_H_
-/*! \brief C++ porting of APC
+/*! \brief Compute the asymmetry-corrected scattering cross-section of a cluster.
+ *
+ * This function computes the product between the geometrical asymmetry parameter and
+ * the scattering cross-section, like `aps()`, but for a cluster of spheres. See Eq. (3.16)
+ * in Borghese, Denti & Saija (2007).
*
* \param zpv: `double ****`
* \param le: `int`
@@ -30,7 +34,11 @@ void apc(
double sqk, double **gapr, std::complex<double> **gapp
);
-/*! \brief C++ porting of APCRA
+/*! \brief Compute the asymmetry-corrected scattering cross-section under random average
+ * conditions.
+ *
+ * This function computes the product between the geometrical asymmetry parameter and
+ * the scattering cross-section of a cluster using the random average directions.
*
* \param zpv: `double ****`
* \param le: `int`
@@ -45,7 +53,9 @@ void apcra(
double **gaprm, std::complex<double> **gappm
);
-/*! \brief C++ porting of CDTP
+/*! \brief Complex inner product.
+ *
+ * This function performs the complex inner product. It is used by `lucin()`.
*
* \param z: `complex<double>`
* \param am: Matrix of complex.
@@ -59,7 +69,7 @@ std::complex<double> cdtp(
int k, int nj
);
-/*! \brief C++ porting of CGEV
+/*! \brief C++ porting of CGEV. QUESTION: description?
*
* \param ipamo: `int`
* \param mu: `int`
@@ -69,7 +79,10 @@ std::complex<double> cdtp(
*/
double cgev(int ipamo, int mu, int l, int m);
-/*! \brief C++ porting of CMS
+/*! \brief Build the multi-centered M-matrix of the cluster.
+ *
+ * This function constructs the multi-centered M-matrix of the cluster, according
+ * to Eq. (5.28) of Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param c1: `C1 *`
@@ -79,7 +92,11 @@ double cgev(int ipamo, int mu, int l, int m);
*/
void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
-/*! \brief C++ porting of CRSM1
+/*! \brief Compute orientation-averaged scattered field intensity.
+ *
+ * This function computes the intensity of the scattered field for the cluster,
+ * averaged on the orientations. It is invoked for IAVM=1 (geometry referred to
+ * the meridional plane). QUESTION: correct?
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -90,7 +107,10 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
*/
void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
-/*! \brief C++ porting of GHIT
+/*! \brief Compute the transfer vector from N2 to N1.
+ *
+ * This function computes the transfer vector going from N2 to N1, using either
+ * Hankel, Bessel or Bessel from origin functions.
*
* \param ihi: `int`
* \param ipamo: `int`
@@ -109,7 +129,10 @@ std::complex<double> ghit(
C1_AddOns *c1ao, C4 *c4, C6 *c6
);
-/*! \brief C++ porting of HJV
+/*! \brief Compute Hankel funtion and Bessel functions.
+ *
+ * This function constructs the Hankel function and the Bessel functions vectors. See
+ * page 331 in Borghese, Denti & Saija (2007).
*
* \param exri: `double` External medium refractive index.
* \param vk: `double` Wave number.
@@ -125,7 +148,10 @@ void hjv(
C1 *c1, C1_AddOns *c1ao, C4 *c4
);
-/*! \brief C++ porting of LUCIN
+/*! \brief Invert the multi-centered M-matrix.
+ *
+ * This function performs the inversion of the multi-centered M-matrix through
+ * LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param nddmst: `const int`
@@ -134,7 +160,11 @@ void hjv(
*/
void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
-/*! \brief C++ porting of MEXTC
+/*! \brief Compute the average extinction cross-section.
+ *
+ * This funbction computes the average extinction cross-section starting
+ * from the definition of the scattering amplitude. See Sec. 3.2.1 of Borghese,
+ * Denti & Saija (2007).
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -144,7 +174,10 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
*/
void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext);
-/*! \brief C++ porting of PCROS
+/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.
+ *
+ * This function computes the scattering, absorption and extinction cross-sections
+ * of the cluster, together with the Forward Scattering Amplitude.
*
* This function is intended to evaluate the particle cross-section. QUESTIUON: correct?
* \param vk: `double` Wave number.
@@ -155,7 +188,10 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
*/
void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4);
-/*! \brief C++ porting of PCRSM0
+/*! \brief Compute orientation-averaged cross-sections and forward scattering amplitude.
+ *
+ * This function computes the orientation-averaged scattering, absorption and extinction
+ * cross-sections of the cluster, together with the averaged Forward Scattering Amplitude.
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -166,7 +202,10 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4);
*/
void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4);
-/*! \brief C++ porting of POLAR
+/*! \brief Transform Cartesian quantities to spherical coordinates ones.
+ *
+ * This function performs a conversion from the Cartesian coordinates system to the
+ * spherical one. It is used by `sphar()`.
*
* \param x: `double` X-axis Cartesian coordinate.
* \param y: `double` Y-axis Cartesian coordinate.
@@ -182,7 +221,10 @@ void polar(
double &cph, double &sph
);
-/*! \brief C++ porting of R3J000
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients towards J=0.
+ *
+ * This function calculates the 3j(J,J2,J3;0,0,0) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j2: `int`
* \param j3: `int`
@@ -190,7 +232,10 @@ void polar(
*/
void r3j000(int j2, int j3, C6 *c6);
-/*! \brief C++ porting of R3JJR
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients for JJ transitions.
+ *
+ * This function calculates the 3j(J,J2,J3;-M2-M3,M2,M3) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j2: `int`
* \param j3: `int`
@@ -200,7 +245,10 @@ void r3j000(int j2, int j3, C6 *c6);
*/
void r3jjr(int j2, int j3, int m2, int m3, C6 *c6);
-/*! \brief C++ porting of R3JMR
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients for JM transitions.
+ *
+ * This function calculates the 3j(J,J2,J3;M1,M,-M1-M) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j1: `int`
* \param j2: `int`
@@ -210,7 +258,12 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6);
*/
void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
-/*! \brief C++ porting of RABA
+/*! \brief Compute radiation torques on a particle in Cartesian coordinates.
+ *
+ * This function computes radiation torque on on a cluster of spheres as the
+ * result of the difference between the extinction and the scattering
+ * contributions for a Cartesian coordinate system, as `rabas()`. See Sec. 4.9
+ * in Borghese, Denti & Saija (2007).
*
* \param le: `int`
* \param am0m: Matrix of complex.
@@ -225,7 +278,10 @@ void raba(
std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
);
-/*! \brief C++ porting of RFTR
+/*! \brief Compute the radiation force Cartesian components.
+ *
+ * This function computes the Cartesian components of the radiation force
+ * exerted on a particle. See Sec. 3.2.1 in Borghese, Denti & Saija (2007).
*
* \param u: `double *`
* \param up: `double *`
@@ -248,7 +304,12 @@ void rftr(
double &fy, double &fz
);
-/*! \brief C++ porting of SCR0
+/*! \brief Compute Mie cross-sections for the sphere units in the cluster.
+ *
+ * This function computes the scattering, absorption and extinction cross-sections
+ * for the spheres composing the cluster, in terms of Mie coefficients, together
+ * with the Forward Scattering Amplitude. See Sec. 4.2.1 in Borghese, Denti & Saija
+ * (2007).
*
* \param vk: `double` Wave number
* \param exri: `double` External medium refractive index.
@@ -259,7 +320,10 @@ void rftr(
*/
void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4);
-/*! \brief C++ porting of SCR2
+/*! \brief Compute the scattering amplitude for a single sphere in an aggregate.
+ *
+ * This function computes the scattering amplitude for single spheres constituting
+ * an aggregate. See Sec. 4.2.1 in Borghese, Denti & Saija (2007).
*
* \param vk: `double` Wave number.
* \param vkarg: `double` QUESTION: definition?
@@ -275,7 +339,11 @@ void scr2(
C3 *c3, C4 *c4
);
-/*! \brief C++ porting of STR
+/*! \brief Transform sphere Cartesian coordinates to spherical coordinates.
+ *
+ * This function transforms the Cartesian coordinates of the spheres in an aggregate
+ * to radial coordinates, then it calls `sphar()` to calculate the vector of spherical
+ * harmonics of the incident field.
*
* \param rcf: `double **` Matrix of double.
* \param c1: `C1 *` Pointer to a C1 instance.
@@ -286,7 +354,12 @@ void scr2(
*/
void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
-/*! \brief C++ porting of TQR
+/*! \brief Compute radiation torques on particles on a k-vector oriented system.
+ *
+ * This function computes the radiation torques resulting from the difference
+ * between absorption and scattering contributions, like `rabas()`, but for a
+ * coordinate system oriented along the wave vector and its orthogonal axis. See
+ * Sec. 4.9 in Borghese, Denti & Saija (2007).
*
* \param u: `double *`
* \param up: `double *`
@@ -305,7 +378,10 @@ void tqr(
double &ten, double &tek, double &tsp, double &tsn, double &tsk
);
-/*! \brief C++ porting of ZTM
+/*! \brief Calculate the single-centered inversion of the M-matrix.
+ *
+ * This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
+ * of Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param c1: `C1 *` Pointer to a C1 instance.
diff --git a/src/include/errors.h b/src/include/errors.h
new file mode 100644
index 0000000000000000000000000000000000000000..8730b36b70c251e93d4d2707b8900bdcd13870e3
--- /dev/null
+++ b/src/include/errors.h
@@ -0,0 +1,161 @@
+/*! \file errors.h
+ *
+ * \brief Collection of proprietary code exceptions.
+ *
+ * There are many circumstances that can prevent the correct execution
+ * of a code. These range from user mistakes, to improper configuration,
+ * to unsupported hardware and all the way up to various system failures.
+ * Although it is not possible to grant proper execution in all cases,
+ * it is often possible to design a code in such a way that the program
+ * detects unexpected conditions, informs the user and takes the proper
+ * actions, eventually stopping without crash, if no other options are
+ * available. C++ handles such unexpected circumstances by means of
+ * `exceptions`. These are special procedures that can be launched
+ * whenever an unexpected situation occurs and they allow to restore the
+ * code work-flow and attempt recovery. Exceptions can be divided in
+ * different cathegories, which respond to various types of problems.
+ * This library contains a set of exceptions designed to the most common
+ * problems that may occur while executing an application of the `NP_TMcode`
+ * suite.
+ */
+
+#ifndef INCLUDE_ERRORS_H_
+#define INCLUDE_ERRORS_H_
+
+/*! \brief Exception for out of bounds List requests.
+ */
+class ListOutOfBoundsException: public std::exception {
+protected:
+ //! Description of the problem.
+ std::string message;
+
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param requested: `int` The index that was requested.
+ * \param min: `int` The minimum index allowed by the list.
+ * \param max: `int` The maximum index allowed by the list.
+ */
+ ListOutOfBoundsException(int requested, int min, int max) {
+ message = "Error: requested index " + std::to_string(requested)
+ + " out of list allowed bounds [" + std::to_string(min) + ", "
+ + std::to_string(max - 1) + "]";
+ }
+
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/*! \brief Exception for open file error handlers.
+ */
+class OpenConfigurationFileException: public std::exception {
+protected:
+ //! \brief Name of the file that was accessed.
+ std::string file_name;
+
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param name: `string` Name of the file that was accessed.
+ */
+ OpenConfigurationFileException(std::string name) { file_name = name; }
+
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return file_name.c_str();
+ }
+};
+
+/*! \brief Exception for access requests out of matrix bounds.
+ */
+class MatrixOutOfBoundsException: 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.
+ */
+ MatrixOutOfBoundsException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/*! \brief Exception for unrecognized configuration data sets.
+ */
+class UnrecognizedConfigurationException: 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.
+ */
+ UnrecognizedConfigurationException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/*! \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 Exception for unrecognized parameters.
+ */
+class UnrecognizedParameterException: 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.
+ */
+ UnrecognizedParameterException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+#endif
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 9526b34341f85fc4a8be09a3ba259b655a05982a..2aafef2233caa43a8f76377afc65db380f056c45 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -14,7 +14,7 @@
#ifndef INCLUDE_SPH_SUBS_H_
#define INCLUDE_SPH_SUBS_H_
-/*! \brief Compute the asymmetry-corrected scattering cross-section.
+/*! \brief Compute the asymmetry-corrected scattering cross-section of a sphere.
*
* This function computes the product between the geometrical asymmetry parameter and
* the scattering cross-section. See Sec. 3.2.1 of Borghese, Denti & Saija (2007).
@@ -179,9 +179,9 @@ void pwma(
int isq, C1 *c1
);
-/*! \brief Compute radiation torques on particles.
+/*! \brief Compute radiation torques on a single sphere in Cartesian coordinates.
*
- * This function computes radiation torque on the particle as a result
+ * This function computes radiation torque on a sphere unit as the result
* of the difference between the extinction and the scattering contributions.
* See Sec. 4.9 in Borghese, Denti & Saija (2007).
*
@@ -321,7 +321,12 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1);
*/
void thdps(int lm, double ****zpv);
-/*! \brief C++ porting of UPVMP
+/*! \brief Compute the unitary vectors onb the polarization plane and its orthogonal
+ * direction.
+ *
+ * This function computes the unitary vectors lying on the polarization plane and on
+ * its orthogonal direction, to optimize the identification of the scattering geometry.
+ * See Sec. 2.3 in Borghese, Denti & Saija (2007).
*
* \param thd: `double`
* \param phd: `double`
diff --git a/src/include/tfrfme.h b/src/include/tfrfme.h
new file mode 100644
index 0000000000000000000000000000000000000000..ba2954a742b5458f8cabb1ca64729befef232e20
--- /dev/null
+++ b/src/include/tfrfme.h
@@ -0,0 +1,442 @@
+/*! \file tfrfme.h
+ *
+ * \brief Representation of the trapping calculation objects.
+ */
+
+#ifndef INCLUDE_TFRFME_H_
+#define INCLUDE_TFRFME_H_
+
+/*! \brief Class to represent the first group of trapping swap data.
+ */
+class Swap1 {
+protected:
+ //! Index of the last element to be filled.
+ int last_index;
+ //! Number of vector coordinates. QUESTION: correct?
+ int nkv;
+ //! NLMMT = 2 * LM * (LM + 2)
+ int nlmmt;
+
+ //! QUESTION: definition?
+ std::complex<double> *wk;
+
+ /*! \brief Load a Swap1 instance from a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `Swap1 *` Pointer to a new Swap1 instance.
+ */
+ static Swap1 *from_hdf5(std::string file_name);
+
+ /*! \brief Load a Swap1 instance from a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `Swap1 *` Pointer to a new Swap1 instance.
+ */
+ static Swap1 *from_legacy(std::string file_name);
+
+ /*! \brief Save a Swap1 instance to a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Save a Swap1 instance to a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_legacy(std::string file_name);
+
+public:
+ /*! \brief Swap1 instance constructor.
+ *
+ * \param lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
+ */
+ Swap1(int lm, int _nkv);
+
+ /*! \brief Swap1 instance destroyer.
+ */
+ ~Swap1() { delete[] wk; }
+
+ /*! \brief Append an element at the end of the vector.
+ *
+ * \param value: `complex<double>` The value to be added to the vector.
+ */
+ void append(std::complex<double> value) { wk[last_index++] = value; }
+
+ /*! \brief Load a Swap1 instance from binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ * \return instance: `Swap1 *` Pointer to a newly created Swap1 instance.
+ */
+ static Swap1* from_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Calculate the necessary amount of memory to create a new instance.
+ *
+ * \param lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
+ * \return size: `long` The necessary memory size in bytes.
+ */
+ static long get_memory_requirement(int lm, int _nkv);
+
+ /*! \brief Get the pointer to the WK vector.
+ *
+ * \return value: `complex<double> *` Memory address of the WK vector.
+ */
+ std::complex<double> *get_vector() { return wk; }
+
+ /*! \brief Bring the pointer to the next element at the start of vector.
+ */
+ void reset() { last_index = 0; }
+
+ /*! \brief Write a Swap1 instance to binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Test whether two instances of Swap1 are equal.
+ *
+ * \param other: `Swap1 &` Reference to the instance to be compared
+ * with.
+ * \return result: `bool` True, if the two instances are equal,
+ * false otherwise.
+ */
+ bool operator ==(Swap1 &other);
+};
+
+/*! \brief Class to represent the second group of trapping swap data.
+ */
+class Swap2 {
+protected:
+ //! Index of the last vector element to be filled.
+ int last_vector;
+ //! Index of the last matrix element to be filled.
+ int last_matrix;
+ //! Number of vector coordinates. QUESTION: correct?
+ int nkv;
+ //! QUESTION: definition?
+ double *vkv;
+ //! QUESTION: definition?
+ double **vkzm;
+ //! QUESTION: definition?
+ double apfafa;
+ //! QUESTION: definition?
+ double pmf;
+ //! QUESTION: definition?
+ double spd;
+ //! QUESTION: definition?
+ double rir;
+ //! QUESTION: definition?
+ double ftcn;
+ //! QUESTION: definition?
+ double fshmx;
+ //! QUESTION: definition?
+ double vxyzmx;
+ //! Cartesian displacement. QUESTION: correct?
+ double delxyz;
+ //! QUESTION: definition?
+ double vknmx;
+ //! QUESTION: definition?
+ double delk;
+ //! QUESTION: definition?
+ double delks;
+ //! NLMMT = LM * (LM + 2) * 2
+ int nlmmt;
+ //! Number of radial vector coordinates. QUESTION: correct?
+ int nrvc;
+
+ /*! \brief Load a Swap2 instance from a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `Swap2 *` Pointer to a new Swap2 instance.
+ */
+ static Swap2 *from_hdf5(std::string file_name);
+
+ /*! \brief Load a Swap2 instance from a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `Swap2 *` Pointer to a new Swap2 instance.
+ */
+ static Swap2 *from_legacy(std::string file_name);
+
+ /*! \brief Save a Swap2 instance to a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Save a Swap2 instance to a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_legacy(std::string file_name);
+
+public:
+ /*! \brief Swap2 instance constructor.
+ *
+ * \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
+ */
+ Swap2(int _nkv);
+
+ /*! \brief Swap2 instance destroyer.
+ */
+ ~Swap2();
+
+ /*! \brief Load a Swap2 instance from binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ * \return instance: `Swap2 *` Pointer to a newly created Swap2 instance.
+ */
+ static Swap2* from_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Get the pointer to the VKZM matrix.
+ *
+ * \return value: `double **` Pointer to the VKZM matrix.
+ */
+ double **get_matrix() { return vkzm; }
+
+ /*! \brief Calculate the necessary amount of memory to create a new instance.
+ *
+ * \param _nkv: `int` Number of radial vector coordinates. QUESTION: correct?
+ * \return size: `long` The necessary memory size in bytes.
+ */
+ static long get_memory_requirement(int _nkv);
+
+ /*! \brief Get a parameter by its name.
+ *
+ * \param param_name: `string` Name of the parameter.
+ * \return value: `double` The value of the requested parameter.
+ */
+ double get_param(std::string param_name);
+
+ /*! \brief Get the pointer to the VKV vector.
+ *
+ * \return value: `double *` Pointer to the VKV vector.
+ */
+ double *get_vector() { return vkv; }
+
+ /*! \brief Append an element at the end of the matrix.
+ *
+ * \param value: `double` The value to be pushed in the matrix.
+ */
+ void push_matrix(double value);
+
+ /*! \brief Append an element at the end of the vector.
+ *
+ * \param value: `double` The value to be pushed in the vector.
+ */
+ void push_vector(double value) { vkv[last_vector++] = value; }
+
+ /*! \brief Bring the matrix pointer to the start of the array.
+ */
+ void reset_matrix() { last_matrix = 0; }
+
+ /*! \brief Bring the vector pointer to the start of the array.
+ */
+ void reset_vector() { last_vector = 0; }
+
+ /*! \brief Set a parameter by its name and value.
+ *
+ * \param param_name: `string` Name of the parameter.
+ * \param value: `double` The value of the parameter.
+ */
+ void set_param(std::string param_name, double value);
+
+ /*! \brief Write a Swap2 instance to binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Test whether two instances of Swap2 are equal.
+ *
+ * \param other: `Swap1 &` Reference to the instance to be compared
+ * with.
+ * \return result: `bool` True, if the two instances are equal,
+ * false otherwise.
+ */
+ bool operator ==(Swap2 &other);
+};
+
+/*! \brief Class to represent the trapping configuration.
+ */
+class TFRFME {
+protected:
+ //! NLMMT = 2 * LM * (LM + 2)
+ int nlmmt;
+ //! NRVC = NXV * NYV * NZV
+ int nrvc;
+
+ //! Field expansion mode identifier.
+ int lmode;
+ //! Maximim field expansion order;
+ int lm;
+ //! QUESTION: definition?
+ int nkv;
+ //! Number of computed X coordinates.
+ int nxv;
+ //! Number of computed Y coordinates.
+ int nyv;
+ //! Number of computed Z coordinates.
+ int nzv;
+ //! Wave number in scale units
+ double vk;
+ //! External medium refractive index
+ double exri;
+ //! QUESTION: definition?
+ double an;
+ //! QUESTION: definition?
+ double ff;
+ //! QUESTION: definition?
+ double tra;
+ //! QUESTION: definition?
+ double spd;
+ //! QUESTION: definition?
+ double frsh;
+ //! QUESTION: definition?
+ double exril;
+ //! Vector of computed x positions
+ double *xv;
+ //! Vector of computed y positions
+ double *yv;
+ //! Vector of computed z positions
+ double *zv;
+ //! QUESTION: definition?
+ std::complex<double> **wsum;
+
+ /*! \brief Load a configuration instance from a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `TFRFME *` Pointer to a new trapping configuration
+ * instance.
+ */
+ static TFRFME *from_hdf5(std::string file_name);
+
+ /*! \brief Load a configuration instance from a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be loaded.
+ * \return instance: `TFRFME *` Pointer to a new trapping configuration
+ * instance.
+ */
+ static TFRFME *from_legacy(std::string file_name);
+
+ /*! \brief Save a configuration instance to a HDF5 binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Save a configuration instance to a legacy binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_legacy(std::string file_name);
+
+public:
+ /*! \brief Trapping configuration instance constructor.
+ *
+ * \param _lmode: `int` Order expansion mode flag.
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of wave vector coordinates. QUESTION: correct?
+ * \param _nxv: `int` Number of computed X coordinates.
+ * \param _nyv: `int` Number of computed Y coordinates.
+ * \param _nzv: `int` Number of computed Z coordinates.
+ */
+ TFRFME(
+ int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv
+);
+
+ /*! \brief Trapping configuration instance destroyer.
+ */
+ ~TFRFME();
+
+ /*! \brief Load a trapping configuration instance from binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ * \return instance: `TFRFME *` Pointer to a newly created configuration
+ * instance.
+ */
+ static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Get the pointer to the WSUM matrix.
+ *
+ * \return value: `complex<double> **` Pointer to the WSUM matrix.
+ */
+ std::complex<double> **get_matrix() { return wsum; }
+
+ /*! \brief Calculate the necessary amount of memory to create a new instance.
+ *
+ * \param _lmode: `int` Order expansion mode flag.
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _nkv: `int` Number of radial vector coordinates. QUESTION: correct?
+ * \param _nxv: `int` Number of computed X coordinates.
+ * \param _nyv: `int` Number of computed Y coordinates.
+ * \param _nzv: `int` Number of computed Z coordinates.
+ * \return size: `long` The necessary memory size in bytes.
+ */
+ static long get_memory_requirement(
+ int _lmode, int _lm, int _nkv, int _nxv,
+ int _nyv, int _nzv
+ );
+
+ /*! \brief Get a configuration parameter.
+ *
+ * \param param_name: `string` Name of the parameter.
+ * \return value: `double` The value of the requested parameter.
+ */
+ double get_param(std::string param_name);
+
+ /*! \brief Get the pointer to the X coordinates vector.
+ *
+ * \return x: `double *` Pointer to X coordinates vector.
+ */
+ double *get_x() { return xv; }
+
+ /*! \brief Get the pointer to the Y coordinates vector.
+ *
+ * \return y: `double *` Pointer to Y coordinates vector.
+ */
+ double *get_y() { return yv; }
+
+ /*! \brief Get the pointer to the Z coordinates vector.
+ *
+ * \return z: `double *` Pointer to Z coordinates vector.
+ */
+ double *get_z() { return zv; }
+
+ /*! \brief Set a configuration parameter.
+ *
+ * \param param_name: `string` Name of the parameter.
+ * \param value: `double` Value to be stored as parameter.
+ */
+ void set_param(std::string param_name, double value);
+
+ /*! \brief Write a trapping configuration instance to binary file.
+ *
+ * \param file_name: `string` Name of the file.
+ * \param mode: `string` Format of the file (can be either "HDF5"
+ * or "LGEACY". Default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode="LEGACY");
+
+ /*! \brief Test whether two instances of configuration are equal.
+ *
+ * \param other: `TFRFME &` Reference to the instance to be compared
+ * with.
+ * \return result: `bool` True, if the two instances are equal,
+ * false otherwise.
+ */
+ bool operator ==(TFRFME &other);
+};
+#endif
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 82594bd0744878f8109bf16003686c2f992a7341..c32aa9ae9bffb58387d936afcf6547c12713be21 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -90,24 +90,10 @@ void ffrf(
double *fffs, CIL *cil, CCR *ccr
);
-/*! C++ porting of FFRT
- *
- * \param ac: Vector of complex. QUESTION: definition?
- * \param ws: Vector of complex. QUESTION: definition?
- * \param ffte: `double *`. QUESTION: definition?
- * \param ffts: `double *`. QUESTION: definition?
- * \param cil: `CIL *` Pointer to a CIL structure.
- */
-void ffrt(
- std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
- CIL *cil
-);
-
/*! C++ porting of FRFMER
*
* \param nkv: `int` QUESTION: definition?
* \param vkm: `double` QUESTION: definition?
- * \param vkv: `double *` QUESTION: definition?
* \param vknmx: `double` QUESTION: definition?
* \param apfafa: `double` QUESTION: definition?
* \param tra: `double` QUESTION: definition?
@@ -117,13 +103,26 @@ void ffrt(
* \param le: `int` QUESTION: definition?
* \param lmode: `int` QUESTION: definition?
* \param pmf: `double` QUESTION: definition?
- * \param tt1: `fstream &` Handle to first temporary binary file.
- * \param tt2: `fstream &` Handle to second temporary binary file.
+ * \param tt1: `Swap1 *` Pointer to first swap object.
+ * \param tt2: `Swap2 *` Pointer to second swap object.
*/
-void frfmer(
- int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
+std::complex<double> *frfmer(
+ int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
- std::fstream &tt1, std::fstream &tt2
+ Swap1 *tt1, Swap2 *tt2
+);
+
+/*! C++ porting of FFRT
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param ffte: `double *`. QUESTION: definition?
+ * \param ffts: `double *`. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ */
+void ffrt(
+ std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
+ CIL *cil
);
/*! C++ porting of PWMALP
@@ -167,7 +166,7 @@ void sampoa(
/*! C++ porting of WAMFF
*
- * \param wk: Vector of complex. QUESTION: definition?
+ * \param wk: `complex<double> *` QUESTION: definition?
* \param x: `double`
* \param y: `double`
* \param z: `double`
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 7c0b63d9bc54135e33d127949490382b64c4664d..e7ed703baf80b284252177a8732fe28930616c8c 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -12,6 +12,10 @@
#include <regex>
#include <string>
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
#endif
@@ -535,7 +539,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
herr_t status = hdf_file->get_status();
string str_name, str_type;
if (status == 0) {
- int nsph;
+ int nsph, ies;
int *iog;
double _exdc, _wp, _xip;
int _idfc, nxi;
@@ -545,6 +549,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
double **rcf_vector;
complex<double> ***dc0m;
status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
+ status = hdf_file->read("IES", "INT32_(1)", &ies);
status = hdf_file->read("EXDC", "FLOAT64_(1)", &_exdc);
status = hdf_file->read("WP", "FLOAT64_(1)", &_wp);
status = hdf_file->read("XIP", "FLOAT64_(1)", &_xip);
@@ -613,7 +618,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
rcf_vector,
_idfc,
dc0m,
- false,
+ (ies == 1),
_exdc,
_wp,
_xip
@@ -709,6 +714,23 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
return conf;
}
+double ScattererConfiguration::get_param(string param_name) {
+ double value;
+ if (param_name.compare("number_of_spheres") == 0) value = 1.0 * number_of_spheres;
+ else if (param_name.compare("nsph") == 0) value = 1.0 * number_of_spheres;
+ else if (param_name.compare("number_of_scales") == 0) value = 1.0 * number_of_scales;
+ else if (param_name.compare("nxi") == 0) value = 1.0 * number_of_scales;
+ else if (param_name.compare("idfc") == 0) value = 1.0 * idfc;
+ else if (param_name.compare("exdc") == 0) value = exdc;
+ else if (param_name.compare("wp") == 0) value = wp;
+ else if (param_name.compare("xip") == 0) value = xip;
+ else {
+ // TODO: add exception code for unknown parameter.
+ return 0.0;
+ }
+ return value;
+}
+
void ScattererConfiguration::print() {
int ies = (use_external_sphere)? 1 : 0;
int configurations = 0;
@@ -782,6 +804,9 @@ void ScattererConfiguration::write_hdf5(string file_name) {
rec_name_list.set(0, "NSPH");
rec_type_list.set(0, "INT32_(1)");
rec_ptr_list.set(0, &number_of_spheres);
+ rec_name_list.append("IES");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&ies);
rec_name_list.append("IOGVEC");
str_type = "INT32_(" + to_string(number_of_spheres) + ")";
rec_type_list.append(str_type);
diff --git a/src/libnptm/Makefile b/src/libnptm/Makefile
index 33dfa35e2a709c5b23381ea367b0e074b19d5f0d..b2fbc322c6b1a227ba0e29e80cd2d8e8c2c688bd 100644
--- a/src/libnptm/Makefile
+++ b/src/libnptm/Makefile
@@ -19,10 +19,11 @@ endif
include ../make.inc
-CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o
+CXX_NPTM_OBJS=$(OBJDIR)/Commons.o $(OBJDIR)/Configuration.o $(OBJDIR)/file_io.o $(OBJDIR)/Parsers.o $(OBJDIR)/sph_subs.o $(OBJDIR)/clu_subs.o $(OBJDIR)/tfrfme.o $(OBJDIR)/tra_subs.o $(OBJDIR)/TransitionMatrix.o
-CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o
+CXX_NPTM_DYNOBJS=$(DYNOBJDIR)/Commons.o $(DYNOBJDIR)/Configuration.o $(DYNOBJDIR)/file_io.o $(DYNOBJDIR)/Parsers.o $(DYNOBJDIR)/sph_subs.o $(DYNOBJDIR)/clu_subs.o $(DYNOBJDIR)/tfrfme.o $(DYNOBJDIR)/tra_subs.o $(DYNOBJDIR)/TransitionMatrix.o
+CXX_NPTM_DEBUG=$(OBJDIR)/Commons.g* $(OBJDIR)/Configuration.g* $(OBJDIR)/file_io.g* $(OBJDIR)/Parsers.g* $(OBJDIR)/sph_subs.g* $(OBJDIR)/clu_subs.g* $(OBJDIR)/tfrfme.g* $(OBJDIR)/tra_subs.g* $(OBJDIR)/TransitionMatrix.g*
all: $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so
@@ -42,10 +43,9 @@ $(BUILDDIR_NPTM)/libnptm.so: $(BUILDDIR_NPTM) $(DYNOBJDIR) $(CXX_NPTM_DYNOBJS)
$(CXX) $(CXXFLAGS) $(PICFLAGS) $(SOFLAGS) -o $(BUILDDIR_NPTM)/libnptm.so $(CXX_NPTM_DYNOBJS)
clean:
- rm -f $(CXX_NPTM_OBJS) $(CXX_NPTM_DYNOBJS)
+ rm -f $(CXX_NPTM_OBJS) $(CXX_NPTM_DYNOBJS) $(CXX_NPTM_DEBUG)
wipe:
#echo "BUILDDIR_NPTM in libnptm is $(BUILDDIR_NPTM)"
#echo "LIBNPTM in libnptm is $(LIBNPTM)"
- rm -f $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so $(CXX_NPTM_OBJS) $(CXX_NPTM_DYNOBJS)
-
+ rm -f $(BUILDDIR_NPTM)/libnptm.a $(BUILDDIR_NPTM)/libnptm.so $(CXX_NPTM_OBJS) $(CXX_NPTM_DYNOBJS) $(CXX_NPTM_DEBUG)
diff --git a/src/libnptm/Parsers.cpp b/src/libnptm/Parsers.cpp
index 3f8c9ff5116126b909baa9b7755f9bad09659642..22ff56db05a692980df8f4140f97a1ba6c227c12 100644
--- a/src/libnptm/Parsers.cpp
+++ b/src/libnptm/Parsers.cpp
@@ -3,10 +3,21 @@
* \brief Implementation of the parsing functions.
*/
+#include <exception>
#include <fstream>
#include <string>
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
+#endif
std::string *load_file(std::string file_name, int *count = 0) {
std::fstream input_file(file_name.c_str(), std::ios::in);
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
index 688607811c261e46142354e5a3f0f9f719529c35..30126ff17ad8cf8511ce395ce3fa0df3f69b80e4 100644
--- a/src/libnptm/TransitionMatrix.cpp
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -6,6 +6,11 @@
#include <exception>
#include <fstream>
#include <hdf5.h>
+#include <string>
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
diff --git a/src/libnptm/file_io.cpp b/src/libnptm/file_io.cpp
index 0c0b0c170009bb40abf290564ee838b68f1260a9..f960e26c6bd3317031d63aae225ff4729d19ef52 100644
--- a/src/libnptm/file_io.cpp
+++ b/src/libnptm/file_io.cpp
@@ -2,11 +2,15 @@
*
* \brief Implementation of file I/O operations.
*/
-#include <stdexcept>
+#include <exception>
#include <regex>
#include <string>
#include <hdf5.h>
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
#endif
@@ -126,7 +130,7 @@ HDFFile* HDFFile::from_schema(
delete[] max_dims;
} else {
string message = "unrecognized type \"" + rec_types[ri] + "\"\n";
- throw runtime_error(message);
+ throw UnrecognizedParameterException(message);
}
}
@@ -159,7 +163,7 @@ herr_t HDFFile::read(
case 2:
mem_type_id = H5T_NATIVE_DOUBLE; break;
default:
- throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ throw UnrecognizedParameterException("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);
@@ -169,7 +173,7 @@ herr_t HDFFile::read(
status = (herr_t)-1;
}
} else {
- throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ throw UnrecognizedParameterException("unrecognized data type \"" + data_type + "\"");
}
return status;
}
@@ -198,7 +202,7 @@ herr_t HDFFile::write(
case 2:
mem_type_id = H5T_NATIVE_DOUBLE; break;
default:
- throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ throw UnrecognizedParameterException("unrecognized data type \"" + data_type + "\"");
}
if (dataset_id != H5I_INVALID_HID) {
status = H5Dwrite(dataset_id, mem_type_id, mem_space_id, file_space_id, dxpl_id, buffer);
@@ -208,7 +212,7 @@ herr_t HDFFile::write(
status = (herr_t)-1;
}
} else {
- throw runtime_error("Unrecognized data type \"" + data_type + "\"");
+ throw UnrecognizedParameterException("unrecognized data type \"" + data_type + "\"");
}
return status;
}
diff --git a/src/libnptm/tfrfme.cpp b/src/libnptm/tfrfme.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b3478d3cc34850289b6a583a6db78fbd0ae0f663
--- /dev/null
+++ b/src/libnptm/tfrfme.cpp
@@ -0,0 +1,1008 @@
+/*! \file tfrfme.cpp
+ *
+ * \brief Implementation of the trapping calculation objects.
+ */
+
+#include <complex>
+#include <exception>
+#include <fstream>
+#include <hdf5.h>
+#include <string>
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_TFRFME_H_
+#include "../include/tfrfme.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
+using namespace std;
+
+// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
+Swap1::Swap1(int lm, int _nkv) {
+ nkv = _nkv;
+ nlmmt = 2 * lm * (lm + 2);
+ const int size = nkv * nkv * nlmmt;
+ wk = new complex<double>[size]();
+ last_index = 0;
+}
+
+Swap1* Swap1::from_binary(string file_name, string mode) {
+ Swap1 *instance = NULL;
+ if (mode.compare("LEGACY") == 0) {
+ instance = from_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ instance = from_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+ return instance;
+}
+
+Swap1* Swap1::from_hdf5(string file_name) {
+ Swap1 *instance = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ double *elements;
+ string str_type;
+ int _nlmmt, _nkv, lm, num_elements, index;
+ complex<double> value;
+ complex<double> *_wk = NULL;
+ if (status == 0) {
+ status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
+ status = hdf_file->read("NKV", "INT32", &_nkv);
+ lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
+ num_elements = 2 * _nlmmt * _nkv * _nkv;
+ instance = new Swap1(lm, _nkv);
+ _wk = instance->get_vector();
+ elements = new double[num_elements]();
+ str_type = "FLOAT64_(" + to_string(num_elements) + ")";
+ status = hdf_file->read("WK", str_type, elements);
+ for (int wi = 0; wi < num_elements / 2; wi++) {
+ index = 2 * wi;
+ value = complex<double>(elements[index], elements[index + 1]);
+ _wk[wi] = value;
+ } // wi loop
+ delete[] elements;
+ status = hdf_file->close();
+ delete hdf_file;
+ }
+ return instance;
+}
+
+Swap1* Swap1::from_legacy(string file_name) {
+ fstream input;
+ Swap1 *instance = NULL;
+ complex<double> *_wk = NULL;
+ int _nlmmt, _nkv, lm;
+ double rval, ival;
+ input.open(file_name.c_str(), ios::in | ios::binary);
+ if (input.is_open()) {
+ input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
+ lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
+ input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
+ instance = new Swap1(lm, _nkv);
+ _wk = instance->get_vector();
+ int num_elements = _nlmmt * _nkv * _nkv;
+ for (int j = 0; j < num_elements; j++) {
+ input.read(reinterpret_cast<char *>(&rval), sizeof(double));
+ input.read(reinterpret_cast<char *>(&ival), sizeof(double));
+ _wk[j] = complex<double>(rval, ival);
+ }
+ input.close();
+ } else {
+ printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
+ }
+ return instance;
+}
+
+long Swap1::get_memory_requirement(int lm, int _nkv) {
+ long size = (long)(3 * sizeof(int));
+ size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2) * _nkv * _nkv);
+ return size;
+}
+
+void Swap1::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 Swap1::write_hdf5(string file_name) {
+ List<string> rec_name_list(1);
+ List<string> rec_type_list(1);
+ List<void *> rec_ptr_list(1);
+ herr_t status;
+ string str_type;
+ int num_elements = 2 * nlmmt * nkv * nkv;
+ rec_name_list.set(0, "NLMMT");
+ rec_type_list.set(0, "INT32_(1)");
+ rec_ptr_list.set(0, &nlmmt);
+ rec_name_list.append("NKV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nkv);
+ rec_name_list.append("WK");
+ str_type = "FLOAT64_(" + to_string(num_elements) + ")";
+ rec_type_list.append(str_type);
+ double *ptr_elements = new double[num_elements]();
+ for (int wi = 0; wi < num_elements / 2; wi++) {
+ ptr_elements[2 * wi] = wk[wi].real();
+ ptr_elements[2 * wi + 1] = wk[wi].imag();
+ }
+ rec_ptr_list.append(ptr_elements);
+
+ 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++)
+ status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ status = hdf_file->close();
+
+ delete[] ptr_elements;
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete hdf_file;
+}
+
+void Swap1::write_legacy(string file_name) {
+ fstream output;
+ double rval, ival;
+ output.open(file_name.c_str(), ios::out | ios::binary);
+ if (output.is_open()) {
+ int num_elements = nlmmt * nkv * nkv;
+ output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
+ output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ for (int j = 0; j < num_elements; j++) {
+ rval = wk[j].real();
+ ival = wk[j].imag();
+ output.write(reinterpret_cast<char *>(&rval), sizeof(double));
+ output.write(reinterpret_cast<char *>(&ival), sizeof(double));
+ }
+ output.close();
+ } else { // Should never occur
+ printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
+ }
+}
+
+bool Swap1::operator ==(Swap1 &other) {
+ if (nlmmt != other.nlmmt) {
+ return false;
+ }
+ if (nkv != other.nkv) {
+ return false;
+ }
+ const int num_elements = nlmmt * nkv * nkv;
+ for (int i = 0; i < num_elements; i++) {
+ if (wk[i] != other.wk[i]) {
+ return false;
+ }
+ }
+ return true;
+}
+// >>> END OF Swap1 CLASS IMPLEMENTATION <<<
+
+// >>> START OF Swap2 CLASS IMPLEMENTATION <<<
+Swap2::Swap2(int _nkv) {
+ nkv = _nkv;
+ vkv = new double[nkv]();
+ vkzm = new double*[nkv];
+ for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv]();
+ last_vector = 0;
+ last_matrix = 0;
+}
+
+Swap2::~Swap2() {
+ delete[] vkv;
+ for (int vi = nkv - 1; vi > -1; vi--) delete[] vkzm[vi];
+ delete[] vkzm;
+}
+
+Swap2* Swap2::from_binary(string file_name, string mode) {
+ Swap2 *instance = NULL;
+ if (mode.compare("LEGACY") == 0) {
+ instance = from_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ instance = from_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+ return instance;
+}
+
+Swap2* Swap2::from_hdf5(string file_name) {
+ Swap2 *instance = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ string str_type;
+ int _nkv, _nlmmt, _nrvc;
+ double value;
+ if (status == 0) {
+ status = hdf_file->read("NKV", "INT32", &_nkv);
+ instance = new Swap2(_nkv);
+ str_type = "FLOAT64_(" + to_string(_nkv) + ")";
+ status = hdf_file->read("VKV", str_type, instance->vkv);
+ str_type = "FLOAT64_(" + to_string(_nkv) + "," + to_string(_nkv) + ")";
+ status = hdf_file->read("VKZM", str_type, instance->vkzm);
+ status = hdf_file->read("APFAFA", "FLOAT64", &value);
+ instance->set_param("apfafa", value);
+ status = hdf_file->read("PMF", "FLOAT64", &value);
+ instance->set_param("pmf", value);
+ status = hdf_file->read("SPD", "FLOAT64", &value);
+ instance->set_param("spd", value);
+ status = hdf_file->read("RIR", "FLOAT64", &value);
+ instance->set_param("rir", value);
+ status = hdf_file->read("FTCN", "FLOAT64", &value);
+ instance->set_param("ftcn", value);
+ status = hdf_file->read("FSHMX", "FLOAT64", &value);
+ instance->set_param("fshmx", value);
+ status = hdf_file->read("VXYZMX", "FLOAT64", &value);
+ instance->set_param("vxyzmx", value);
+ status = hdf_file->read("DELXYZ", "FLOAT64", &value);
+ instance->set_param("delxyz", value);
+ status = hdf_file->read("VKNMX", "FLOAT64", &value);
+ instance->set_param("vknmx", value);
+ status = hdf_file->read("DELK", "FLOAT64", &value);
+ instance->set_param("delk", value);
+ status = hdf_file->read("DELKS", "FLOAT64", &value);
+ instance->set_param("delks", value);
+ status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
+ instance->set_param("nlmmt", 1.0 * _nlmmt);
+ status = hdf_file->read("NRVC", "INT32", &_nrvc);
+ instance->set_param("nrvc", 1.0 * _nrvc);
+ status = hdf_file->close();
+ delete hdf_file;
+ }
+ return instance;
+}
+
+Swap2* Swap2::from_legacy(string file_name) {
+ fstream input;
+ Swap2 *instance = NULL;
+ int _nkv, _nlmmt, _nrvc;
+ double **_vkzm = NULL;
+ double *_vkv = NULL;
+ double value;
+ input.open(file_name.c_str(), ios::in | ios::binary);
+ if (input.is_open()) {
+ input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
+ instance = new Swap2(_nkv);
+ _vkzm = instance->get_matrix();
+ _vkv = instance->get_vector();
+ for (int vj = 0; vj < _nkv; vj++) {
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ _vkv[vj] = value;
+ }
+ for (int mi = 0; mi < _nkv; mi++) {
+ for (int mj = 0; mj < _nkv; mj++) {
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ _vkzm[mi][mj] = value;
+ }
+ }
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("apfafa", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("pmf", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("spd", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("rir", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("ftcn", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("fshmx", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("vxyzmx", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("delxyz", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("vknmx", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("delk", value);
+ input.read(reinterpret_cast<char *>(&value), sizeof(double));
+ instance->set_param("delks", value);
+ input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
+ instance->set_param("nlmmt", 1.0 * _nlmmt);
+ input.read(reinterpret_cast<char *>(&_nrvc), sizeof(int));
+ instance->set_param("nrvc", 1.0 * _nrvc);
+ input.close();
+ } else {
+ printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
+ }
+ return instance;
+}
+
+long Swap2::get_memory_requirement(int _nkv) {
+ long size = (long)(5 * sizeof(int) + 11 * sizeof(double));
+ size += (long)(sizeof(double) * _nkv * (_nkv + 1));
+ return size;
+}
+
+double Swap2::get_param(string param_name) {
+ double value;
+ if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
+ else if (param_name.compare("apfafa") == 0) value = apfafa;
+ else if (param_name.compare("pmf") == 0) value = pmf;
+ else if (param_name.compare("spd") == 0) value = spd;
+ else if (param_name.compare("rir") == 0) value = rir;
+ else if (param_name.compare("ftcn") == 0) value = ftcn;
+ else if (param_name.compare("fshmx") == 0) value = fshmx;
+ else if (param_name.compare("vxyzmx") == 0) value = vxyzmx;
+ else if (param_name.compare("delxyz") == 0) value = delxyz;
+ else if (param_name.compare("vknmx") == 0) value = vknmx;
+ else if (param_name.compare("delk") == 0) value = delk;
+ else if (param_name.compare("delks") == 0) value = delks;
+ else if (param_name.compare("nlmmt") == 0) value = 1.0 * nlmmt;
+ else if (param_name.compare("nrvc") == 0) value = 1.0 * nrvc;
+ else {
+ string message = "Unrecognized parameter name \"" + param_name + "\"";
+ throw UnrecognizedParameterException(message);
+ }
+ return value;
+}
+
+void Swap2::push_matrix(double value) {
+ int col = last_matrix % (nkv - 1);
+ int row = last_matrix - (nkv * row);
+ vkzm[row][col] = value;
+ last_matrix++;
+}
+
+void Swap2::set_param(string param_name, double value) {
+ if (param_name.compare("nkv") == 0) nkv = (int)value;
+ else if (param_name.compare("apfafa") == 0) apfafa = value;
+ else if (param_name.compare("pmf") == 0) pmf = value;
+ else if (param_name.compare("spd") == 0) spd = value;
+ else if (param_name.compare("rir") == 0) rir = value;
+ else if (param_name.compare("ftcn") == 0) ftcn = value;
+ else if (param_name.compare("fshmx") == 0) fshmx = value;
+ else if (param_name.compare("vxyzmx") == 0) vxyzmx = value;
+ else if (param_name.compare("delxyz") == 0) delxyz = value;
+ else if (param_name.compare("vknmx") == 0) vknmx = value;
+ else if (param_name.compare("delk") == 0) delk = value;
+ else if (param_name.compare("delks") == 0) delks = value;
+ else if (param_name.compare("nlmmt") == 0) nlmmt = (int)value;
+ else if (param_name.compare("nrvc") == 0) nrvc = (int)value;
+ else {
+ string message = "Unrecognized parameter name \"" + param_name + "\"";
+ throw UnrecognizedParameterException(message);
+ }
+}
+
+void Swap2::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 Swap2::write_hdf5(string file_name) {
+ List<string> rec_name_list(1);
+ List<string> rec_type_list(1);
+ List<void *> rec_ptr_list(1);
+ herr_t status;
+ string str_type;
+ rec_name_list.set(0, "NKV");
+ rec_type_list.set(0, "INT32_(1)");
+ rec_ptr_list.set(0, &nkv);
+ rec_name_list.append("VKV");
+ str_type = "FLOAT64_(" + to_string(nkv) + ")";
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(vkv);
+ rec_name_list.append("VKZM");
+ str_type = "FLOAT64_(" + to_string(nkv) + "," + to_string(nkv) + ")";
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(vkzm);
+ rec_name_list.append("APFAFA");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&apfafa);
+ rec_name_list.append("PMF");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&pmf);
+ rec_name_list.append("SPD");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&spd);
+ rec_name_list.append("RIR");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&rir);
+ rec_name_list.append("FTCN");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&ftcn);
+ rec_name_list.append("FSHMX");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&fshmx);
+ rec_name_list.append("VXYZMX");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&vxyzmx);
+ rec_name_list.append("delxyz");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&delxyz);
+ rec_name_list.append("VKNMX");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&vknmx);
+ rec_name_list.append("DELK");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&delk);
+ rec_name_list.append("DELKS");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&delks);
+ rec_name_list.append("NLMMT");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nlmmt);
+ rec_name_list.append("NRVC");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nrvc);
+
+ 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++)
+ status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ status = hdf_file->close();
+
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete hdf_file;
+}
+
+void Swap2::write_legacy(string file_name) {
+ fstream output;
+ double value;
+ output.open(file_name.c_str(), ios::out | ios::binary);
+ if (output.is_open()) {
+ output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ for (int j = 0; j < nkv; j++){
+ value = vkv[j];
+ output.write(reinterpret_cast<const char*>(&value), sizeof(double));
+ }
+ for (int mi = 0; mi < nkv; mi++) {
+ for (int mj = 0; mj < nkv; mj++) {
+ value = vkzm[mi][mj];
+ output.write(reinterpret_cast<const char*>(&value), sizeof(double));
+ }
+ }
+ output.write(reinterpret_cast<const char*>(&apfafa), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&pmf), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&spd), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&rir), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&ftcn), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&fshmx), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&vxyzmx), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&delxyz), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&vknmx), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&delk), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&delks), sizeof(double));
+ output.write(reinterpret_cast<const char*>(&nlmmt), sizeof(int));
+ output.write(reinterpret_cast<const char*>(&nrvc), sizeof(int));
+ output.close();
+ } else { // Should never occur
+ printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
+ }
+}
+
+bool Swap2::operator ==(Swap2 &other) {
+ if (nlmmt != other.nlmmt) {
+ return false;
+ }
+ if (nrvc != other.nrvc) {
+ return false;
+ }
+ if (nkv != other.nkv) {
+ return false;
+ }
+ if (apfafa != other.apfafa) {
+ return false;
+ }
+ if (pmf != other.pmf) {
+ return false;
+ }
+ if (spd != other.spd) {
+ return false;
+ }
+ if (rir != other.rir) {
+ return false;
+ }
+ if (ftcn != other.ftcn) {
+ return false;
+ }
+ if (fshmx != other.fshmx) {
+ return false;
+ }
+ if (vxyzmx != other.vxyzmx) {
+ return false;
+ }
+ if (delxyz != other.delxyz) {
+ return false;
+ }
+ if (vknmx != other.vknmx) {
+ return false;
+ }
+ if (delk != other.delk) {
+ return false;
+ }
+ if (delks != other.delks) {
+ return false;
+ }
+ for (int vi = 0; vi < nkv; vi++) {
+ if (vkv[vi] != other.vkv[vi]) {
+ return false;
+ }
+ }
+ for (int mi = 0; mi < nkv; mi++) {
+ for (int mj = 0; mj < nkv; mj++) {
+ if (vkzm[mi][mj] != other.vkzm[mi][mj]) {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+// >>> END OF Swap2 CLASS IMPLEMENTATION <<<
+
+// >>> START OF TFRFME CLASS IMPLEMENTATION <<<
+TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) {
+ lmode = _lmode;
+ lm = _lm;
+ nkv = _nkv;
+ nxv = _nxv;
+ nyv = _nyv;
+ nzv = _nzv;
+ vk = 0.0;
+ exri = 0.0;
+ an = 0.0;
+ ff = 0.0;
+ tra = 0.0;
+ spd = 0.0;
+ frsh = 0.0;
+ exril = 0.0;
+
+ // Array initialization
+ xv = new double[nxv]();
+ yv = new double[nyv]();
+ zv = new double[nzv]();
+ nlmmt = lm * (lm + 2) * 2;
+ nrvc = nxv * nyv * nzv;
+ wsum = new complex<double>*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new complex<double>[nrvc]();
+}
+
+TFRFME::~TFRFME() {
+ delete[] xv;
+ delete[] yv;
+ delete[] zv;
+ for (int wi = nlmmt - 1; wi > -1; wi--) delete[] wsum[wi];
+ delete[] wsum;
+}
+
+TFRFME* TFRFME::from_binary(string file_name, string mode) {
+ TFRFME *instance = NULL;
+ if (mode.compare("LEGACY") == 0) {
+ instance = from_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ instance = from_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+ return instance;
+}
+
+TFRFME* TFRFME::from_hdf5(string file_name) {
+ TFRFME *instance = NULL;
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(file_name, flags);
+ herr_t status = hdf_file->get_status();
+ double *elements;
+ string str_type;
+ int _nlmmt, _nrvc, num_elements;
+ complex<double> value;
+ complex<double> **_wsum = NULL;
+ if (status == 0) {
+ int lmode, lm, nkv, nxv, nyv, nzv;
+ double vk, exri, an, ff, tra, spd, frsh, exril;
+ status = hdf_file->read("LMODE", "INT32", &lmode);
+ status = hdf_file->read("LM", "INT32", &lm);
+ status = hdf_file->read("NKV", "INT32", &nkv);
+ status = hdf_file->read("NXV", "INT32", &nxv);
+ status = hdf_file->read("NYV", "INT32", &nyv);
+ status = hdf_file->read("NZV", "INT32", &nzv);
+ status = hdf_file->read("VK", "FLOAT64", &vk);
+ status = hdf_file->read("EXRI", "FLOAT64", &exri);
+ status = hdf_file->read("AN", "FLOAT64", &an);
+ status = hdf_file->read("FF", "FLOAT64", &ff);
+ status = hdf_file->read("TRA", "FLOAT64", &tra);
+ status = hdf_file->read("SPD", "FLOAT64", &spd);
+ status = hdf_file->read("FRSH", "FLOAT64", &frsh);
+ status = hdf_file->read("EXRIL", "FLOAT64", &exril);
+ instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ _wsum = instance->get_matrix();
+ instance->set_param("vk", vk);
+ instance->set_param("exri", exri);
+ instance->set_param("an", an);
+ instance->set_param("ff", ff);
+ instance->set_param("tra", tra);
+ instance->set_param("spd", spd);
+ instance->set_param("frsh", frsh);
+ instance->set_param("exril", exril);
+ str_type = "FLOAT64_(" + to_string(nxv) + ")";
+ status = hdf_file->read("XVEC", str_type, instance->xv);
+ str_type = "FLOAT64_(" + to_string(nyv) + ")";
+ status = hdf_file->read("YVEC", str_type, instance->yv);
+ str_type = "FLOAT64_(" + to_string(nzv) + ")";
+ status = hdf_file->read("ZVEC", str_type, instance->zv);
+ _nlmmt = 2 * lm * (lm + 2);
+ _nrvc = nxv * nyv * nzv;
+ num_elements = 2 * _nlmmt * _nrvc;
+ elements = new double[num_elements]();
+ str_type = "FLOAT64_(" + to_string(num_elements) + ")";
+ status = hdf_file->read("WSUM", str_type, elements);
+ int index = 0;
+ for (int wj = 0; wj < _nrvc; wj++) {
+ for (int wi = 0; wi < _nlmmt; wi++) {
+ value = complex<double>(elements[index], elements[index + 1]);
+ _wsum[wi][wj] = value;
+ index += 2;
+ } // wi loop
+ } // wj loop
+ delete[] elements;
+ status = hdf_file->close();
+ delete hdf_file;
+ }
+ return instance;
+}
+
+TFRFME* TFRFME::from_legacy(string file_name) {
+ fstream input;
+ TFRFME *instance = NULL;
+ complex<double> **_wsum = NULL;
+ double *coord_vec = NULL;
+ input.open(file_name.c_str(), ios::in | ios::binary);
+ if (input.is_open()) {
+ int lmode, lm, nkv, nxv, nyv, nzv;
+ double vk, exri, an, ff, tra, spd, frsh, exril;
+ double dval, rval, ival;
+ input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+ input.read(reinterpret_cast<char *>(&lm), sizeof(int));
+ input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+ input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+ input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+ input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+ input.read(reinterpret_cast<char *>(&vk), sizeof(double));
+ input.read(reinterpret_cast<char *>(&exri), sizeof(double));
+ input.read(reinterpret_cast<char *>(&an), sizeof(double));
+ input.read(reinterpret_cast<char *>(&ff), sizeof(double));
+ input.read(reinterpret_cast<char *>(&tra), sizeof(double));
+ input.read(reinterpret_cast<char *>(&spd), sizeof(double));
+ input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+ input.read(reinterpret_cast<char *>(&exril), sizeof(double));
+ instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ _wsum = instance->get_matrix();
+ instance->set_param("vk", vk);
+ instance->set_param("exri", exri);
+ instance->set_param("an", an);
+ instance->set_param("ff", ff);
+ instance->set_param("tra", tra);
+ instance->set_param("spd", spd);
+ instance->set_param("frsh", frsh);
+ instance->set_param("exril", exril);
+ coord_vec = instance->get_x();
+ for (int xi = 0; xi < nxv; xi++) {
+ input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+ coord_vec[xi] = dval;
+ }
+ coord_vec = instance->get_y();
+ for (int yi = 0; yi < nyv; yi++) {
+ input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+ coord_vec[yi] = dval;
+ }
+ coord_vec = instance->get_z();
+ for (int zi = 0; zi < nzv; zi++) {
+ input.read(reinterpret_cast<char *>(&dval), sizeof(double));
+ coord_vec[zi] = dval;
+ }
+ int _nlmmt = 2 * lm * (lm + 2);
+ int _nrvc = nxv * nyv * nzv;
+ for (int wj = 0; wj < _nrvc; wj++) {
+ for (int wi = 0; wi < _nlmmt; wi++) {
+ input.read(reinterpret_cast<char *>(&rval), sizeof(double));
+ input.read(reinterpret_cast<char *>(&ival), sizeof(double));
+ complex<double> value(rval, ival);
+ _wsum[wi][wj] = value;
+ } // wi loop
+ } // wj loop
+ input.close();
+ } else {
+ printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
+ }
+ return instance;
+}
+
+long TFRFME::get_memory_requirement(
+ int _lmode, int _lm, int _nkv, int _nxv,
+ int _nyv, int _nzv
+) {
+ int _nlmmt = 2 * _lm * (_lm + 2);
+ int _nrvc = _nxv * _nyv * _nzv;
+ long size = (long)(8 * sizeof(int) + 8 * sizeof(double));
+ size += (long)((_nxv + _nyv + _nzv) * sizeof(double));
+ size += (long)(_nlmmt * _nrvc * sizeof(complex<double>));
+ return size;
+}
+
+double TFRFME::get_param(string param_name) {
+ double value;
+ if (param_name.compare("vk") == 0) value = vk;
+ else if (param_name.compare("exri") == 0) value = exri;
+ else if (param_name.compare("an") == 0) value = an;
+ else if (param_name.compare("ff") == 0) value = ff;
+ else if (param_name.compare("tra") == 0) value = tra;
+ else if (param_name.compare("spd") == 0) value = spd;
+ else if (param_name.compare("frsh") == 0) value = frsh;
+ else if (param_name.compare("exril") == 0) value = exril;
+ else if (param_name.compare("lmode") == 0) value = 1.0 * lmode;
+ else if (param_name.compare("lm") == 0) value = 1.0 * lm;
+ else if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
+ else if (param_name.compare("nxv") == 0) value = 1.0 * nxv;
+ else if (param_name.compare("nyv") == 0) value = 1.0 * nyv;
+ else if (param_name.compare("nzv") == 0) value = 1.0 * nzv;
+ else {
+ string message = "Unrecognized parameter name \"" + param_name + "\"";
+ throw UnrecognizedParameterException(message);
+ }
+ return value;
+}
+
+void TFRFME::set_param(string param_name, double value) {
+ if (param_name.compare("vk") == 0) vk = value;
+ else if (param_name.compare("exri") == 0) exri = value;
+ else if (param_name.compare("an") == 0) an = value;
+ else if (param_name.compare("ff") == 0) ff = value;
+ else if (param_name.compare("tra") == 0) tra = value;
+ else if (param_name.compare("spd") == 0) spd = value;
+ else if (param_name.compare("frsh") == 0) frsh = value;
+ else if (param_name.compare("exril") == 0) exril = value;
+ else {
+ string message = "Unrecognized parameter name \"" + param_name + "\"";
+ throw UnrecognizedParameterException(message);
+ }
+}
+
+void TFRFME::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 TFRFME::write_hdf5(string file_name) {
+ List<string> rec_name_list(1);
+ List<string> rec_type_list(1);
+ List<void *> rec_ptr_list(1);
+ herr_t status;
+ string str_type;
+ rec_name_list.set(0, "LMODE");
+ rec_type_list.set(0, "INT32_(1)");
+ rec_ptr_list.set(0, &lmode);
+ rec_name_list.append("LM");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&lm);
+ rec_name_list.append("NKV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nkv);
+ rec_name_list.append("NXV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nxv);
+ rec_name_list.append("NYV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nyv);
+ rec_name_list.append("NZV");
+ rec_type_list.append("INT32_(1)");
+ rec_ptr_list.append(&nzv);
+ 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("AN");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&an);
+ rec_name_list.append("FF");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&ff);
+ rec_name_list.append("TRA");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&tra);
+ rec_name_list.append("SPD");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&spd);
+ rec_name_list.append("FRSH");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&frsh);
+ rec_name_list.append("EXRIL");
+ rec_type_list.append("FLOAT64_(1)");
+ rec_ptr_list.append(&exril);
+ str_type = "FLOAT64_(" + to_string(nxv) + ")";
+ rec_name_list.append("XVEC");
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(xv);
+ str_type = "FLOAT64_(" + to_string(nyv) + ")";
+ rec_name_list.append("YVEC");
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(yv);
+ str_type = "FLOAT64_(" + to_string(nzv) + ")";
+ rec_name_list.append("ZVEC");
+ rec_type_list.append(str_type);
+ rec_ptr_list.append(zv);
+ rec_name_list.append("WSUM");
+ int num_elements = 2 * nlmmt * nrvc;
+ str_type = "FLOAT64_(" + to_string(num_elements) + ")";
+ rec_type_list.append(str_type);
+ // The (N x M) matrix of complex is converted to a vector of double
+ // with length (2 * N * M)
+ double *ptr_elements = new double[num_elements]();
+ int index = 0;
+ for (int wj = 0; wj < nrvc; wj++) {
+ for (int wi = 0; wi < nlmmt; wi++) {
+ ptr_elements[index++] = wsum[wi][wj].real();
+ ptr_elements[index++] = wsum[wi][wj].imag();
+ } // wi loop
+ } // wj loop
+ rec_ptr_list.append(ptr_elements);
+
+ 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++)
+ status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ status = hdf_file->close();
+
+ delete[] ptr_elements;
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete hdf_file;
+}
+
+void TFRFME::write_legacy(string file_name) {
+ fstream output;
+ output.open(file_name.c_str(), ios::out | ios::binary);
+ if (output.is_open()) {
+ output.write(reinterpret_cast<char *>(&lmode), sizeof(int));
+ output.write(reinterpret_cast<char *>(&lm), sizeof(int));
+ output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ output.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+ output.write(reinterpret_cast<char *>(&nyv), sizeof(int));
+ output.write(reinterpret_cast<char *>(&nzv), sizeof(int));
+ output.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ output.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ output.write(reinterpret_cast<char *>(&an), sizeof(double));
+ output.write(reinterpret_cast<char *>(&ff), sizeof(double));
+ output.write(reinterpret_cast<char *>(&tra), sizeof(double));
+ output.write(reinterpret_cast<char *>(&spd), sizeof(double));
+ output.write(reinterpret_cast<char *>(&frsh), sizeof(double));
+ output.write(reinterpret_cast<char *>(&exril), sizeof(double));
+ for (int xi = 0; xi < nxv; xi++)
+ output.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
+ for (int yi = 0; yi < nyv; yi++)
+ output.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
+ for (int zi = 0; zi < nzv; zi++)
+ output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
+ for (int wj = 0; wj < nrvc; wj++) {
+ for (int wi = 0; wi < nlmmt; wi++) {
+ double rval = wsum[wi][wj].real();
+ double ival = wsum[wi][wj].imag();
+ output.write(reinterpret_cast<char *>(&rval), sizeof(double));
+ output.write(reinterpret_cast<char *>(&ival), sizeof(double));
+ } // wi loop
+ } // wj loop
+ output.close();
+ } else { // Should never occur
+ printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
+ }
+}
+
+bool TFRFME::operator ==(TFRFME &other) {
+ if (lmode != other.lmode) {
+ return false;
+ }
+ if (lm != other.lm) {
+ return false;
+ }
+ if (nkv != other.nkv) {
+ return false;
+ }
+ if (nxv != other.nxv) {
+ return false;
+ }
+ if (nyv != other.nyv) {
+ return false;
+ }
+ if (nzv != other.nzv) {
+ return false;
+ }
+ if (vk != other.vk) {
+ return false;
+ }
+ if (exri != other.exri) {
+ return false;
+ }
+ if (an != other.an) {
+ return false;
+ }
+ if (ff != other.ff) {
+ return false;
+ }
+ if (tra != other.tra) {
+ return false;
+ }
+ if (spd != other.spd) {
+ return false;
+ }
+ if (frsh != other.frsh) {
+ return false;
+ }
+ if (exril != other.exril) {
+ return false;
+ }
+ for (int xi = 0; xi < nxv; xi++) {
+ if (xv[xi] != other.xv[xi]) {
+ return false;
+ }
+ }
+ for (int yi = 0; yi < nyv; yi++) {
+ if (yv[yi] != other.yv[yi]) {
+ return false;
+ }
+ }
+ for (int zi = 0; zi < nzv; zi++) {
+ if (zv[zi] != other.zv[zi]) {
+ return false;
+ }
+ }
+ for (int wi = 0; wi < nlmmt; wi++) {
+ for (int wj = 0; wj < nrvc; wj++) {
+ if (wsum[wi][wj] != other.wsum[wi][wj]) {
+ return false;
+ }
+ } // wj loop
+ } // wi loop
+ return true;
+}
+// >>> END OF TFRFME CLASS IMPLEMENTATION <<<
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index aef2de470a5c3fe0b68aa95f05c0dbdac9d63615..47c2b28edfc1fb784f06712f8f70baff9830b5ad 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -14,6 +14,10 @@
#include "../include/sph_subs.h"
#endif
+#ifndef INCLUDE_TFRFME_H_
+#include "../include/tfrfme.h"
+#endif
+
#ifndef INCLUDE_TRA_SUBS_H_
#include "../include/tra_subs.h"
#endif
@@ -77,7 +81,8 @@ void ffrf(
int ltpo = lpo + l40;
int imm = l40 * lpo;
for (int ilmp40 = 1; ilmp40 <= 3; ilmp40++) {
- if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) continue; // ilmp40 loop
+ // NOTE: if trapping execution diverges, replace "break" with "continue"
+ if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) break; // ilmp40 loop
int lmpml = ilmp40 - 2;
int lmp = l40 + lmpml;
uimmp = uim * (-1.0 * lmpml);
@@ -130,7 +135,8 @@ void ffrf(
int ltpo = lpo + l80;
int imm = l80 * lpo;
for (int ilmp80 = 1; ilmp80 <= 3; ilmp80++) {
- if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) continue; // ilmp80 loop
+ // NOTE: if trapping execution diverges, replace "break" with "continue"
+ if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) break; // ilmp80 loop
int lmpml = ilmp80 - 2;
int lmp = l80 + lmpml;
uimmp = uim * (-1.0 * lmpml);
@@ -240,46 +246,40 @@ void ffrt(
delete[] ctqcs;
}
-void frfmer(
- int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
+complex<double> *frfmer(
+ int nkv, double vkm, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
- std::fstream &tt1, std::fstream &tt2
+ Swap1 *tt1, Swap2 *tt2
) {
const int nlemt = le * (le + 2) * 2;
const complex<double> cc0(0.0, 0.0);
- complex<double> *wk = new complex<double>[nlemt]();
double sq = vkm * vkm;
+ double *_vkv = tt2->get_vector();
+ double **_vkzm = tt2->get_matrix();
+ complex<double> *wk = new complex<double>[nlemt]();
for (int jy90 = 0; jy90 < nkv; jy90++) {
- double vky = vkv[jy90];
+ double vky = _vkv[jy90];
double sqy = vky * vky;
for (int jx80 = 0; jx80 < nkv; jx80++) {
- double vkx = vkv[jx80];
+ double vkx = _vkv[jx80];
double sqx = vkx * vkx;
double sqn = sqx + sqy;
double vkn = sqrt(sqn);
if (vkn <= vknmx) {
double vkz = sqrt(sq - sqn);
wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
- for (int j = 0; j < nlemt; j++) {
- double vreal = wk[j].real();
- double vimag = wk[j].imag();
- tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ for (int j = 0; j < nlemt; j++) tt1->append(wk[j]);
+ //tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ _vkzm[jx80][jy90] = vkz;
} else { // label 50
- for (int j = 0; j < nlemt; j++) {
- double vreal = 0.0;
- double vimag = 0.0;
- tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- double vkz = 0.0;
- tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ for (int j = 0; j < nlemt; j++) tt1->append(cc0);
+ //double vkz = 0.0;
+ //tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ _vkzm[jx80][jy90] = 0.0;
}
} // jx80 loop
} // jy90 loop
- delete[] wk;
+ return wk;
}
void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
@@ -384,27 +384,29 @@ void wamff(
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
- if (!(onz && onx && ony)) {
- rho = sqrt(x * x + y * y);
- cph = x / rho;
- sph = y / rho;
- } else {
- if (onz) { // label 10
- cph = 1.0;
- sph = 0.0;
- } else {
- if (onx) { // label 12
- rho = sqrt(x * x);
- cph = (x < 0.0)? -1.0 : 1.0;
- sph = 0.0;
- } else {
- if (ony) { // label 13
- rho = sqrt(y * y);
- cph = 0.0;
- sph = (y < 0.0)? -1.0: 1.0;
- }
+ if (!onz) {
+ if (!onx) {
+ if (!ony) {
+ rho = sqrt(x * x + y * y);
+ cph = x / rho;
+ sph = y / rho;
+ // goes to 15
+ } else { // label 13
+ rho = (y < 0.0) ? -y : y;
+ cph = 0.0;
+ sph = (y < 0.0) ? -1.0 : 1.0;
+ // goes to 15
}
+ } else { // label 12
+ rho = (x < 0.0) ? -x : x;
+ cph = (x < 0.0) ? -1.0 : 1.0;
+ sph = 0.0;
+ // goes to 15
}
+ } else { // label 10
+ cph = 1.0;
+ sph = 0.0;
+ // goes to 15
}
// label 15
if (z == 0.0) {
diff --git a/src/make.inc b/src/make.inc
index c56fd1334e83be3d50ed2363a3ccfc8d2df375b2..dcdd896f0619de451fcf37eaeb62ec664f21c8e6 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -40,8 +40,8 @@ endif
# CXXFLAGS defines the default compilation options for the C++ compiler
ifndef CXXFLAGS
-#override CXXFLAGS=-O0 -ggdb -pg -coverage -I$(HDF5_INCLUDE)
-override CXXFLAGS=-O3 -I$(HDF5_INCLUDE)
+override CXXFLAGS=-O3 -ggdb -pg -coverage -I$(HDF5_INCLUDE)
+#override CXXFLAGS=-O3 -I$(HDF5_INCLUDE)
endif
# HDF5_LIB defines the default path to the HDF5 libraries to use
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index d5f9e499c7d92ae1c3bf732152ee14c8c397442d..be6a0d95582085bebf7815c25aeb23c56633ced1 100644
--- a/src/sphere/Makefile
+++ b/src/sphere/Makefile
@@ -24,6 +24,8 @@ 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 $(OBJDIR)/TransitionMatrix.o
CXX_SPH_OBJS=$(OBJDIR)/np_sphere.o $(OBJDIR)/sphere.o
+CXX_SPH_DEBUG=$(OBJDIR)/np_sphere.g* $(OBJDIR)/sphere.g*
+
all: $(LIBNPTM) $(BUILDDIR_SPH)/edfb_sph $(BUILDDIR_SPH)/sph $(BUILDDIR_SPH)/np_sphere
$(OBJDIR):
@@ -46,8 +48,8 @@ $(BUILDDIR_SPH)/np_sphere: $(OBJDIR) $(CXX_SPH_OBJS) $(BUILDDIR_SPH) $(LIBNPTM)
$(CXX) $(CXXFLAGS) -o $(BUILDDIR_SPH)/np_sphere $(CXX_SPH_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
clean:
- rm -f $(F_SPH_OBJS) $(CXX_SPH_OBJS)
+ rm -f $(F_SPH_OBJS) $(CXX_SPH_OBJS) $(CXX_SPH_DEBUG)
wipe:
- rm -f $(BUILDDIR_SPH)/edfb_sph $(BUILDDIR_SPH)/sph $(BUILDDIR_SPH)/np_sphere $(F_SPH_OBJS) $(CXX_SPH_OBJS)
+ rm -f $(BUILDDIR_SPH)/edfb_sph $(BUILDDIR_SPH)/sph $(BUILDDIR_SPH)/np_sphere $(F_SPH_OBJS) $(CXX_SPH_OBJS) $(CXX_SPH_DEBUG)
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index efad7ad1819f879ab0128b2ed08e10ca7fd4b33f..7016d8361eb8d61fef8437f504516352f0444948 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -8,6 +8,10 @@
#include <fstream>
#include <string>
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
diff --git a/src/testing/Makefile b/src/testing/Makefile
index 88d265d815863bf2eba17ca56eeafc2969506be0..8598b2ef955c0d025b31a7e3a32ab99406ec04ab 100644
--- a/src/testing/Makefile
+++ b/src/testing/Makefile
@@ -21,7 +21,9 @@ include ../make.inc
CXX_TEDF_OBJS=$(OBJDIR)/test_TEDF.o
-CXX_TTMS_OBJS=$(OBJDIR)/test_TTMS.o
+CXX_TTMS_OBJS=$(OBJDIR)/test_TTMS.o
+
+CXX_TESTING_DEBUG=$(OBJDIR)/test_TEDF.g* $(OBJDIR)/test_TTMS.g*
all: $(LIBNPTM) $(BUILDDIR_TEST)/test_TEDF $(BUILDDIR_TEST)/test_TTMS
@@ -43,7 +45,7 @@ $(BUILDDIR_TEST)/test_TTMS: $(OBJDIR) $(CXX_TTMS_OBJS) $(BUILDDIR_TEST) $(LIBNPT
$(CXX) $(CXXFLAGS) -o $(BUILDDIR_TEST)/test_TTMS $(CXX_TTMS_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
clean:
- rm -f $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS)
+ rm -f $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS) $(CXX_TESTING_DEBUG)
wipe:
- rm -f $(BUILDDIR_TEST)/test_TEDF $(BUILDDIR_TEST)/test_TTMS $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS)
+ rm -f $(BUILDDIR_TEST)/test_TEDF $(BUILDDIR_TEST)/test_TTMS $(CXX_TEDF_OBJS) $(CXX_TTMS_OBJS) $(CXX_TESTING_DEBUG)
diff --git a/src/testing/test_TEDF.cpp b/src/testing/test_TEDF.cpp
index b384cc1df0d0faa9982ec12a45a3b8322c615d24..98a681d3895a63659978cd130352083a7028c6e7 100644
--- a/src/testing/test_TEDF.cpp
+++ b/src/testing/test_TEDF.cpp
@@ -2,11 +2,25 @@
#include <complex>
#include <cstdio>
+#include <exception>
#include <hdf5.h>
#include <string>
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
#include "../include/file_io.h"
+#endif
+
+#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
+#endif
using namespace std;
diff --git a/src/testing/test_TTMS.cpp b/src/testing/test_TTMS.cpp
index 3b054359c0ef366f74b9f10b5b023f0bca45f7f4..7c55d16a7939a67ba211f90644b23ba2d329852e 100644
--- a/src/testing/test_TTMS.cpp
+++ b/src/testing/test_TTMS.cpp
@@ -2,11 +2,25 @@
#include <complex>
#include <cstdio>
+#include <exception>
#include <hdf5.h>
#include <string>
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
#include "../include/file_io.h"
+#endif
+
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
#include "../include/TransitionMatrix.h"
+#endif
using namespace std;
diff --git a/src/trapping/Makefile b/src/trapping/Makefile
index ec6059e8eb2c1e71aead588a9d0ab090a563cab5..e7df7dd73631be2bee66add65385f55335103354 100644
--- a/src/trapping/Makefile
+++ b/src/trapping/Makefile
@@ -23,6 +23,8 @@ F_TRAP_OBJS=$(OBJDIR)/frfme.o $(OBJDIR)/lffft.o
CXX_TRAP_OBJS=$(OBJDIR)/np_trapping.o $(OBJDIR)/cfrfme.o $(OBJDIR)/clffft.o
+CXX_TRAP_DEBUG=$(OBJDIR)/np_trapping.g* $(OBJDIR)/cfrfme.g* $(OBJDIR)/clffft.g*
+
all: $(LIBNPTM) $(BUILDDIR_TRA)/frfme $(BUILDDIR_TRA)/lffft $(BUILDDIR_TRA)/np_trapping
$(OBJDIR):
@@ -45,8 +47,8 @@ $(BUILDDIR_TRA)/np_trapping: $(OBJDIR) $(CXX_TRAP_OBJS) $(BUILDDIR_TRA) $(LIBNPT
$(CXX) $(CXXFLAGS) -o $(BUILDDIR_TRA)/np_trapping $(CXX_TRAP_OBJS) $(LIBNPTM) $(CXXLDFLAGS)
clean:
- rm -f $(F_TRAP_OBJS) $(CXX_TRAP_OBJS)
+ rm -f $(F_TRAP_OBJS) $(CXX_TRAP_OBJS) $(CXX_TRAP_DEBUG)
wipe:
- rm -f $(BUILDDIR_TRA)/frfme $(BUILDDIR_TRA)/lffft $(BUILDDIR_TRA)/np_trapping $(F_TRAP_OBJS) $(CXX_TRAP_OBJS)
+ rm -f $(BUILDDIR_TRA)/frfme $(BUILDDIR_TRA)/lffft $(BUILDDIR_TRA)/np_trapping $(F_TRAP_OBJS) $(CXX_TRAP_OBJS) $(CXX_TRAP_DEBUG)
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index f5d7af1cfefb7a5a1d41591962acbc9fbf3a0ae4..c4866e1f65073bcc2f61ec54d7a23a125ac315ca 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -4,6 +4,7 @@
*/
#include <complex>
#include <cstdio>
+#include <exception>
#include <fstream>
#include <regex>
#include <string>
@@ -16,10 +17,18 @@
#include "../include/Commons.h"
#endif
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+
#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif
+#ifndef INCLUDE_TFRFME_H_
+#include "../include/tfrfme.h"
+#endif
+
#ifndef INCLUDE_TRA_SUBS_H_
#include "../include/tra_subs.h"
#endif
@@ -32,13 +41,14 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void frfme(string data_file, string output_path) {
- string tfrfme_name = output_path + "/c_TFRFME";
- fstream tfrfme;
+ string tfrfme_name = output_path + "/c_TFRFME.hd5";
+ TFRFME *tfrfme = NULL;
+ Swap1 *tt1 = NULL;
+ Swap2 *tt2 = NULL;
char namef[7];
char more;
- double *xv = NULL, *yv = NULL, *zv = NULL;
- double *vkv = NULL, **vkzm = NULL;
- complex<double> *wk = NULL, **w = NULL, **wsum = NULL;
+ complex<double> **w = NULL;
+ complex<double> *wk = NULL;
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
int line_count = 0, last_read_line = 0;
@@ -58,73 +68,51 @@ void frfme(string data_file, string output_path) {
double apfafa = 0.0, pmf = 0.0, spd = 0.0, rir = 0.0, ftcn = 0.0, fshmx = 0.0;
double vxyzmx = 0.0, delxyz = 0.0, vknmx = 0.0, delk = 0.0, delks = 0.0;
double frsh = 0.0, exril = 0.0;
+ double **vkzm = NULL;
+ double *vkv = NULL;
int nlmmt = 0, nrvc = 0;
// Vector size variables
int wsum_size;
// End of vector size variables
if (jlmf != 1) {
int nxv, nyv, nzv;
- tfrfme.open(tfrfme_name, ios::in | ios::binary);
- if (tfrfme.is_open()) {
- tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nkv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nxv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nyv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nzv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&vk), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&exri), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&an), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&ff), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&tra), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
- vkv = new double[nkv]();
- xv = new double[nxv]();
- yv = new double[nyv]();
- zv = new double[nzv]();
- for (int xi = 0; xi < nxv; xi++) tfrfme.read(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
- for (int yi = 0; yi < nxv; yi++) tfrfme.read(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
- for (int zi = 0; zi < nxv; zi++) tfrfme.read(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
- fstream temptape2;
- string tempname2 = output_path + "/c_TEMPTAPE2";
- temptape2.open(tempname2.c_str(), ios::in | ios::binary);
- if (temptape2.is_open()) {
- for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
- vkzm = new double*[nkv];
- for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();
- for (int jy10 = 0; jy10 < nkv; jy10++) {
- for (int jx10 = 0; jx10 < nkv; jx10++) {
- temptape2.read(reinterpret_cast<char *>(&(vkzm[jx10][jy10])), sizeof(double));
- } //jx10 loop
- } // jy10 loop
- temptape2.read(reinterpret_cast<char *>(&apfafa), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&pmf), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&spd), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&rir), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&ftcn), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&fshmx), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&vxyzmx), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&delxyz), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&vknmx), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&delk), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&delks), sizeof(double));
- temptape2.read(reinterpret_cast<char *>(&nlmmt), sizeof(int));
- temptape2.read(reinterpret_cast<char *>(&nrvc), sizeof(int));
- temptape2.close();
+ if (tfrfme == NULL) tfrfme = TFRFME::from_binary(tfrfme_name, "HDF5");
+ if (tfrfme != NULL) {
+ lmode = (int)tfrfme->get_param("lmode");
+ lm = (int)tfrfme->get_param("lm");
+ nkv = (int)tfrfme->get_param("nkv");
+ nxv = (int)tfrfme->get_param("nxv");
+ nyv = (int)tfrfme->get_param("nyv");
+ nzv = (int)tfrfme->get_param("nzv");
+ vk = tfrfme->get_param("vk");
+ exri = tfrfme->get_param("exri");
+ an = tfrfme->get_param("an");
+ ff = tfrfme->get_param("ff");
+ tra = tfrfme->get_param("tra");
+ spd = tfrfme->get_param("spd");
+ frsh = tfrfme->get_param("frsh");
+ exril = tfrfme->get_param("exril");
+ string tempname2 = output_path + "/c_TEMPTAPE2.hd5";
+ if (tt2 == NULL) tt2 = Swap2::from_binary(tempname2, "HDF5");
+ if (tt2 != NULL) {
+ vkv = tt2->get_vector();
+ vkzm = tt2->get_matrix();
+ apfafa = tt2->get_param("apfafa");
+ pmf = tt2->get_param("pmf");
+ spd = tt2->get_param("spd");
+ rir = tt2->get_param("rir");
+ ftcn = tt2->get_param("ftcn");
+ fshmx = tt2->get_param("fshmx");
+ vxyzmx = tt2->get_param("vxyzmx");
+ delxyz = tt2->get_param("delxyz");
+ vknmx = tt2->get_param("vknmx");
+ delk = tt2->get_param("delk");
+ delks = tt2->get_param("delks");
+ nlmmt = (int)tt2->get_param("nlmmt");
+ nrvc = (int)tt2->get_param("nrvc");
} else {
printf("ERROR: could not open TEMPTAPE2 file.\n");
}
- for (int ixyz12 = 0; ixyz12 < nrvc; ixyz12++) {
- for (int j12 = 0; j12 < jlmf - 1; j12++) {
- double vreal, vimag;
- tfrfme.read(reinterpret_cast<char *>(&vreal), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- wsum[j12][ixyz12] = complex<double>(vreal, vimag);
- } // j12 loop
- } // ixyz12 loop
- tfrfme.close();
} else {
printf("ERROR: could not open TFRFME file.\n");
}
@@ -179,23 +167,22 @@ void frfme(string data_file, string output_path) {
re = regex("[0-9]+");
regex_search(str_target, m, re);
int ixi = stoi(m.str());
- fstream tedf;
- string tedf_name = output_path + "/" + namef;
- tedf.open(tedf_name.c_str(), ios::in | ios::binary);
- if (tedf.is_open()) {
+ string tedf_name = output_path + "/" + namef + ".hd5";
+ ScattererConfiguration *tedf = ScattererConfiguration::from_binary(tedf_name, "HDF5");
+ if (tedf != NULL) {
int iduml, idum;
- tedf.read(reinterpret_cast<char *>(&iduml), sizeof(int));
- for (int i = 0; i < iduml; i++) tedf.read(reinterpret_cast<char *>(&idum), sizeof(int));
- tedf.read(reinterpret_cast<char *>(&exdc), sizeof(double));
- tedf.read(reinterpret_cast<char *>(&wp), sizeof(double));
- tedf.read(reinterpret_cast<char *>(&xip), sizeof(double));
- tedf.read(reinterpret_cast<char *>(&idfc), sizeof(int));
- tedf.read(reinterpret_cast<char *>(&nxi), sizeof(int));
+ iduml = (int)tedf->get_param("nsph");
+ idum = tedf->get_iog(iduml - 1);
+ exdc = tedf->get_param("exdc");
+ wp = tedf->get_param("wp");
+ xip = tedf->get_param("xip");
+ idfc = (int)tedf->get_param("idfc");
+ nxi = (int)tedf->get_param("nxi");
if (idfc >= 0) {
if (ixi <= nxi) {
- for (int i = 0; i < ixi; i++) tedf.read(reinterpret_cast<char *>(&xi), sizeof(double));
+ xi = tedf->get_scale(ixi - 1);
} else { // label 96
- tedf.close();
+ delete tedf;
// label 98
string output_name = output_path + "/c_OFRFME";
FILE *output = fopen(output_name.c_str(), "w");
@@ -206,7 +193,7 @@ void frfme(string data_file, string output_path) {
xi = xip;
}
// label 20
- tedf.close();
+ delete tedf;
double wn = wp / 3.0e8;
vk = xi * wn;
exri = sqrt(exdc);
@@ -230,9 +217,18 @@ void frfme(string data_file, string output_path) {
nks = nksh * 2;
nkv = nks + 1;
// Array initialization
- vkv = new double[nkv]();
- vkzm = new double*[nkv];
- for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv];
+ long swap1_size, swap2_size, tfrfme_size;
+ double size_mb;
+ printf("INFO: calculation memory requirements\n");
+ swap1_size = Swap1::get_memory_requirement(lm, nkv);
+ size_mb = 1.0 * swap1_size / 1024.0 / 1024.0;
+ printf("Swap 1: %.2lg MB\n", size_mb);
+ swap2_size = Swap2::get_memory_requirement(nkv);
+ size_mb = 1.0 * swap2_size / 1024.0 / 1024.0;
+ printf("Swap 2: %.2lg MB\n", size_mb);
+ tt2 = new Swap2(nkv);
+ vkv = tt2->get_vector();
+ vkzm = tt2->get_matrix();
// End of array initialization
double vkm = vk * exri;
vknmx = vk * an;
@@ -244,26 +240,33 @@ void frfme(string data_file, string output_path) {
int nxs = nxsh * 2;
int nxv = nxs + 1;
int nxshpo = nxsh + 1;
- xv = new double[nxv]();
- for (int i24 = nxshpo; i24 <= nxs; i24++) {
- xv[i24] = xv[i24 - 1] + delxyz;
- xv[nxv - i24 - 1] = -xv[i24];
- } // i24 loop
int nys = nysh * 2;
int nyv = nys + 1;
int nyshpo = nysh + 1;
- yv = new double[nyv]();
- for (int i25 = nyshpo; i25 <= nys; i25++) {
- yv[i25] = yv[i25 - 1] + delxyz;
- yv[nyv - i25 - 1] = -yv[i25];
- } // i25 loop
int nzs = nzsh * 2;
int nzv = nzs + 1;
int nzshpo = nzsh + 1;
- zv = new double[nzv]();
+ tfrfme_size = TFRFME::get_memory_requirement(lmode, lm, nkv, nxv, nyv, nzv);
+ size_mb = 1.0 * tfrfme_size / 1024.0 / 1024.0;
+ printf("TFRFME: %.2lg MB\n", size_mb);
+ size_mb = 1.0 * (swap1_size + swap2_size + tfrfme_size) / 1024.0 / 1024.0;
+ printf("TOTAL: %.2lg MB\n", size_mb);
+ tfrfme = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+ double *_xv = tfrfme->get_x();
+ double *_yv = tfrfme->get_y();
+ double *_zv = tfrfme->get_z();
+ complex<double> **_wsum = tfrfme->get_matrix();
+ for (int i24 = nxshpo; i24 <= nxs; i24++) {
+ _xv[i24] = _xv[i24 - 1] + delxyz;
+ _xv[nxv - i24 - 1] = -_xv[i24];
+ } // i24 loop
+ for (int i25 = nyshpo; i25 <= nys; i25++) {
+ _yv[i25] = _yv[i25 - 1] + delxyz;
+ _yv[nyv - i25 - 1] = -_yv[i25];
+ } // i25 loop
for (int i27 = nzshpo; i27 <= nzs; i27++) {
- zv[i27] = zv[i27 - 1] + delxyz;
- zv[nzv - i27 - 1] = -zv[i27];
+ _zv[i27] = _zv[i27 - 1] + delxyz;
+ _zv[nzv - i27 - 1] = -_zv[i27];
} // i27 loop
int nrvc = nxv * nyv * nzv;
int nkshpo = nksh + 1;
@@ -271,71 +274,41 @@ void frfme(string data_file, string output_path) {
vkv[i28] = vkv[i28 - 1] + delk;
vkv[nkv - i28 - 1] = -vkv[i28];
} // i28 loop
- tfrfme.open(tfrfme_name.c_str(), ios::out | ios::binary);
- if (tfrfme.is_open()) {
- tfrfme.write(reinterpret_cast<char *>(&lmode), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&lm), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&nkv), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&nxv), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&nyv), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&nxv), sizeof(int));
- tfrfme.write(reinterpret_cast<char *>(&vk), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&exri), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&an), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&ff), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&tra), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&spd), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&frsh), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&exril), sizeof(double));
- for (int xi = 0; xi < nxv; xi++)
- tfrfme.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
- for (int yi = 0; yi < nyv; yi++)
- tfrfme.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
- for (int zi = 0; zi < nzv; zi++)
- tfrfme.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
- fstream temptape1, temptape2;
- string temp_name1 = output_path + "/c_TEMPTAPE1";
- string temp_name2 = output_path + "/c_TEMPTAPE2";
- temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
- temptape2.open(temp_name2.c_str(), ios::out | ios::binary);
- for (int jx = 0; jx < nkv; jx++)
- temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
- frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, temptape1, temptape2);
- temptape1.close();
- temptape2.write(reinterpret_cast<char *>(&apfafa), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&pmf), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&spd), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&rir), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&ftcn), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&fshmx), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&vxyzmx), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&delxyz), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&vknmx), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&delk), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&delks), sizeof(double));
- temptape2.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
- temptape2.write(reinterpret_cast<char *>(&nrvc), sizeof(int));
- temptape2.close();
- temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);
- for (int jx = 0; jx < nkv; jx++) {
- double value = 0.0;
- temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
- vkv[jx] = value;
- }
- for (int jy40 = 0; jy40 < nkv; jy40++) {
- for (int jx40 = 0; jx40 < nkv; jx40++) {
- double value = 0.0;
- temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
- vkzm[jx40][jy40] = value;
- }
- } // jy40 loop
- temptape2.close();
- if (wsum != NULL) {
- for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
- delete[] wsum;
- }
- wsum = new complex<double>*[nlmmt];
- for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+ if (tfrfme != NULL) {
+ tfrfme->set_param("vk", vk);
+ tfrfme->set_param("exri", exri);
+ tfrfme->set_param("an", an);
+ tfrfme->set_param("ff", ff);
+ tfrfme->set_param("tra", tra);
+ tfrfme->set_param("spd", spd);
+ tfrfme->set_param("frsh", frsh);
+ tfrfme->set_param("exril", exril);
+ //fstream temptape1, temptape2;
+ string temp_name1 = output_path + "/c_TEMPTAPE1.hd5";
+ string temp_name2 = output_path + "/c_TEMPTAPE2.hd5";
+ //temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
+ tt1 = new Swap1(lm, nkv);
+ if (wk != NULL) delete[] wk;
+ wk = frfmer(nkv, vkm, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, tt1, tt2);
+ tt1->write_binary(temp_name1, "HDF5");
+ //delete tt1;
+ tt2->set_param("apfafa", apfafa);
+ tt2->set_param("pmf", pmf);
+ tt2->set_param("spd", spd);
+ tt2->set_param("rir", rir);
+ tt2->set_param("ftcn", ftcn);
+ tt2->set_param("fshmx", fshmx);
+ tt2->set_param("vxyzmx", vxyzmx);
+ tt2->set_param("delxyz", delxyz);
+ tt2->set_param("vknmx", vknmx);
+ tt2->set_param("delk", delk);
+ tt2->set_param("delks", delks);
+ tt2->set_param("nlmmt", 1.0 * nlmmt);
+ tt2->set_param("nrvc", 1.0 * nrvc);
+ tt2->write_binary(temp_name2, "HDF5");
+ for (int j80 = jlmf; j80 <= jlml; j80++) {
+ complex<double> *tt1_wk = tt1->get_vector();
+ int wk_index = 0;
// w matrix
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
@@ -343,36 +316,27 @@ void frfme(string data_file, string output_path) {
}
w = new complex<double>*[nkv];
for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
- if (wk != NULL) delete[] wk;
- wk = new complex<double>[nlmmt]();
- temptape1.open(temp_name1.c_str(), ios::in | ios::binary);
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
- for (int i = 0; i < nlmmt; i++) {
- double vreal, vimag;
- temptape1.read(reinterpret_cast<char *>(&vreal), sizeof(double));
- temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double));
- wk[i] = complex<double>(vreal, vimag);
- }
- w[jx50][jy50] = wk[j80];
+ for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
+ w[jx50][jy50] = wk[j80 - 1];
} // jx50
} // jy50 loop
- temptape1.close();
int ixyz = 0;
- wsum[j80] = new complex<double>[nrvc]();
+ for (int wj = 0; wj < nrvc; wj++) _wsum[j80 - 1][wj] = cc0;
for (int iz75 = 0; iz75 < nzv; iz75++) {
- double z = zv[iz75] + frsh;
+ double z = _zv[iz75] + frsh;
for (int iy70 = 0; iy70 < nyv; iy70++) {
- double y = yv[iy70];
+ double y = _yv[iy70];
for (int ix65 = 0; ix65 < nxv; ix65++) {
- double x = xv[ix65];
+ double x = _xv[ix65];
ixyz++;
complex<double> sumy = cc0;
for (int jy60 = 0; jy60 < nkv; jy60++) {
double vky = vkv[jy60];
double vkx = vkv[nkv - 1];
double vkzf = vkzm[0][jy60];
- complex<double> phasf = exp(uim * (-vkx * x + vky * y +vkzf * z));
+ complex<double> phasf = exp(uim * (-vkx * x + vky * y + vkzf * z));
double vkzl = vkzm[nkv - 1][jy60];
complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));
complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
@@ -385,41 +349,13 @@ void frfme(string data_file, string output_path) {
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
sumy += sumx;
} // jy60 loop
- wsum[j80][ixyz - 1] = sumy * delks;
+ _wsum[j80 - 1][ixyz - 1] = sumy * delks;
} // ix65 loop
} // iy70 loop
} // iz75 loop
} // j80 loop
- if (jlmf != 1) {
- tfrfme.open(tfrfme_name, ios::in | ios::out | ios::binary);
- tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nkv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nxv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nyv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&nzv), sizeof(int));
- tfrfme.read(reinterpret_cast<char *>(&vk), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&exri), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&an), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&ff), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&tra), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
- tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
- for (int i = 0; i < nxv; i++) tfrfme.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
- for (int i = 0; i < nyv; i++) tfrfme.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
- for (int i = 0; i < nzv; i++) tfrfme.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
- }
// label 88
- for (int ixyz = 0; ixyz < nrvc; ixyz++) {
- for (int j = 0; j< jlml; j++) {
- double vreal = wsum[j][ixyz].real();
- double vimag = wsum[j][ixyz].imag();
- tfrfme.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tfrfme.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- } // j loop
- } // ixyz loop
- tfrfme.close();
+ tfrfme->write_binary(tfrfme_name, "HDF5");
string output_name = output_path + "/c_OFRFME";
FILE *output = fopen(output_name.c_str(), "w");
fprintf(output, " IF JLML < NLMMT, PRESERVE TEMPTAPE1, TEMPTAPE2, AND TFRFRME,\n");
@@ -441,24 +377,14 @@ void frfme(string data_file, string output_path) {
}
}
// label 45
- if (tfrfme.is_open()) tfrfme.close();
+ if (tfrfme != NULL) delete tfrfme;
delete[] file_lines;
- if (xv != NULL) delete[] xv;
- if (yv != NULL) delete[] yv;
- if (zv != NULL) delete[] zv;
- if (vkv != NULL) delete[] vkv;
- if (vkzm != NULL) {
- for (int vki = nkv - 1; vki > -1; vki--) delete[] vkzm[vki];
- delete[] vkzm;
- }
+ if (tt2 != NULL) delete tt2;
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
- if (wsum != NULL) {
- for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
- delete[] wsum;
- }
if (wk != NULL) delete[] wk;
- printf("Done.\n");
+ if (tt1 != NULL) delete tt1;
+ printf("FRFME: Done.\n");
}
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index 9cd9b50b409b48383fb02d74983ef761a5c71506..7e3619a252ada2d7102a4d88b3b2a2ba97e60f19 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -4,6 +4,7 @@
*/
#include <complex>
#include <cstdio>
+#include <exception>
#include <fstream>
#include <regex>
#include <string>
@@ -20,6 +21,10 @@
#include "../include/sph_subs.h"
#endif
+#ifndef INCLUDE_TFRFME_H_
+#include "../include/tfrfme.h"
+#endif
+
#ifndef INCLUDE_TRA_SUBS_H_
#include "../include/tra_subs.h"
#endif
@@ -38,12 +43,12 @@ void lffft(string data_file, string output_path) {
fstream tlfff, tlfft;
double ****zpv = NULL;
- double *xv = NULL, *yv = NULL, *zv = NULL;
complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;
complex<double> **am0m = NULL;
complex<double> **amd = NULL;
int **indam = NULL;
complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
+ complex<double> **_wsum = NULL;
int jft, jss, jtw;
int is, le, nvam = 0;
double vks, exris;
@@ -140,37 +145,35 @@ void lffft(string data_file, string output_path) {
}
// label 150
ttms.close();
- fstream binary_input;
+ TFRFME *tfrfme = NULL;
string binary_name;
- if (jss != 1) binary_name = output_path + "/c_TFRFME";
+ if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";
else binary_name = output_path + "/c_TWS";
- binary_input.open(binary_name, ios::in | ios::binary);
- if (binary_input.is_open()) {
+ tfrfme = TFRFME::from_binary(binary_name, "HDF5");
+ if (tfrfme != NULL) {
int lmode, lm, nkv, nxv, nyv, nzv;
double vk, exri, an, ff, tra;
double spd, frsh, exril;
- binary_input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
- binary_input.read(reinterpret_cast<char *>(&lm), sizeof(int));
- binary_input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
- binary_input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
- binary_input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
- binary_input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+ lmode = (int)tfrfme->get_param("lmode");
+ lm = (int)tfrfme->get_param("lm");
+ nkv = (int)tfrfme->get_param("nkv");
+ nxv = (int)tfrfme->get_param("nxv");
+ nyv = (int)tfrfme->get_param("nyv");
+ nzv = (int)tfrfme->get_param("nzv");
+ _wsum = tfrfme->get_matrix();
+ double *_xv = tfrfme->get_x();
+ double *_yv = tfrfme->get_y();
+ double *_zv = tfrfme->get_z();
if (lm >= le) {
- binary_input.read(reinterpret_cast<char *>(&vk), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&exri), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&an), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&ff), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&tra), sizeof(double));
+ vk = tfrfme->get_param("vk");
+ exri = tfrfme->get_param("exri");
+ an = tfrfme->get_param("an");
+ ff = tfrfme->get_param("ff");
+ tra = tfrfme->get_param("tra");
if (vk == vks && exri == exris) {
- binary_input.read(reinterpret_cast<char *>(&spd), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&exril), sizeof(double));
- xv = new double[nxv];
- for (int i = 0; i < nxv; i++) binary_input.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
- yv = new double[nyv];
- for (int i = 0; i < nyv; i++) binary_input.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
- zv = new double[nzv];
- for (int i = 0; i < nzv; i++) binary_input.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ spd = tfrfme->get_param("spd");
+ frsh = tfrfme->get_param("frsh");
+ exril = tfrfme->get_param("exril");
bool goto160 = false;
if (jft <= 0) {
zpv = new double***[le];
@@ -204,12 +207,18 @@ void lffft(string data_file, string output_path) {
tlfff.write(reinterpret_cast<char *>(&spd), sizeof(double));
tlfff.write(reinterpret_cast<char *>(&frsh), sizeof(double));
tlfff.write(reinterpret_cast<char *>(&exril), sizeof(double));
- for (int i = 0; i < nxv; i++)
- tlfff.write(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
- for (int i = 0; i < nyv; i++)
- tlfff.write(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
- for (int i = 0; i < nzv; i++)
- tlfff.write(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ for (int i = 0; i < nxv; i++) {
+ double x = _xv[i];
+ tlfff.write(reinterpret_cast<char *>(&x), sizeof(double));
+ }
+ for (int i = 0; i < nyv; i++) {
+ double y = _yv[i];
+ tlfff.write(reinterpret_cast<char *>(&y), sizeof(double));
+ }
+ for (int i = 0; i < nzv; i++) {
+ double z = _zv[i];
+ tlfff.write(reinterpret_cast<char *>(&z), sizeof(double));
+ }
if (jft < 0) goto160 = true;
} else { // Should never happen.
printf("ERROR: could not open TLFFF file.\n");
@@ -236,12 +245,18 @@ void lffft(string data_file, string output_path) {
tlfft.write(reinterpret_cast<char *>(&spd), sizeof(double));
tlfft.write(reinterpret_cast<char *>(&frsh), sizeof(double));
tlfft.write(reinterpret_cast<char *>(&exril), sizeof(double));
- for (int i = 0; i < nxv; i++)
- tlfft.write(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
- for (int i = 0; i < nyv; i++)
- tlfft.write(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
- for (int i = 0; i < nzv; i++)
- tlfft.write(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ for (int i = 0; i < nxv; i++) {
+ double x = _xv[i];
+ tlfft.write(reinterpret_cast<char *>(&x), sizeof(double));
+ }
+ for (int i = 0; i < nyv; i++) {
+ double y = _yv[i];
+ tlfft.write(reinterpret_cast<char *>(&y), sizeof(double));
+ }
+ for (int i = 0; i < nzv; i++) {
+ double z = _zv[i];
+ tlfft.write(reinterpret_cast<char *>(&z), sizeof(double));
+ }
} else { // Should never happen.
printf("ERROR: could not open TLFFT file.\n");
}
@@ -262,13 +277,18 @@ void lffft(string data_file, string output_path) {
for (int iy475 = 0; iy475 < nyv; iy475++) {
for (int ix475 = 0; ix475 < nxv; ix475++) {
for (int i = 0; i < nlmmt; i++) {
- double vreal, vimag;
- binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));
- binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ //double vreal, vimag;
+ //binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ //binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ int row = i;
+ int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
+ complex<double> value = _wsum[row][col];
if (lm <= le) {
- ws[i] = complex<double>(vreal, vimag);
+ //ws[i] = complex<double>(vreal, vimag);
+ ws[i] = value;
} else { // label 170
- wsl[i] = complex<double>(vreal, vimag);
+ //wsl[i] = complex<double>(vreal, vimag);
+ wsl[i] = value;
for (int i175 = 0; i175 < cil->nlem; i175++) {
int ie = i175 + cil->nlem;
int iel = i175 + nlmm;
@@ -384,7 +404,7 @@ void lffft(string data_file, string output_path) {
fclose(output67);
}
}
- binary_input.close();
+ delete tfrfme;
} else {
printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());
}
@@ -395,9 +415,6 @@ void lffft(string data_file, string output_path) {
// Clean up memory
if (ac != NULL) delete[] ac;
if (ws != NULL) delete[] ws;
- if (xv != NULL) delete[] xv;
- if (yv != NULL) delete[] yv;
- if (zv != NULL) delete[] zv;
if (wsl != NULL) delete[] wsl;
if (tmsm != NULL) delete[] tmsm;
if (tmse != NULL) delete[] tmse;
@@ -430,5 +447,5 @@ void lffft(string data_file, string output_path) {
delete cil;
delete ccr;
delete[] file_lines;
- printf("Done.\n");
+ printf("LFFT: Done.\n");
}