diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 49e0c3b3e5585b36e762a8eb83a3b6a26e7f102f..c8aa1da3c685a489a4fcb2857c592723a2773ff0 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -2,8 +2,8 @@
*
* \brief Implementation of the calculation for a cluster of spheres.
*/
-#include <cstdio>
#include <complex>
+#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
diff --git a/src/include/TransitionMatrix.h b/src/include/TransitionMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..3e08852e2a302778a78c1f16f225100a21b12729
--- /dev/null
+++ b/src/include/TransitionMatrix.h
@@ -0,0 +1,99 @@
+/*! \file TransitionMatrix.h
+ *
+ * \brief Representation of the Transition Matrix.
+ */
+
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#define INCLUDE_TRANSITIONMATRIX_H_
+
+/**
+ * \brief Exception for unrecognized file formats.
+ */
+class UnrecognizedFormatException: public std::exception {
+protected:
+ //! Description of the problem.
+ std::string message;
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param problem: `string` Description of the problem that occurred.
+ */
+ UnrecognizedFormatException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/*! \brief Class to represent the Transition Matrix.
+ */
+class TransitionMatrix {
+ protected:
+ //! Matrix type identifier.
+ int is;
+ //! Maximum field expansion order.
+ int l_max;
+ //! Wave number in scale units.
+ double vk;
+ //! External medium refractive index.
+ double exri;
+ //! Vectorized matrix elements.
+ std::complex<double> *elements;
+ //! Sphere radius.
+ double sphere_radius;
+ //! Matrix shape
+ int *shape;
+
+ /*! \brief Write the Transition Matrix to HDF5 binary output.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_hdf5(std::string file_name);
+
+ /*! \brief Write the Transition Matrix to legacy binary output.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ */
+ void write_legacy(std::string file_name);
+ public:
+ /*! \brief Transition Matrix instance constructor for single sphere.
+ *
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _vk: `double` Wave number in scale units.
+ * \param _exri: `double` External medium refractive index.
+ * \param _rmi: Matrix of complex.
+ * \param _rei: Matrix of complex.
+ * \param _sphere_radius: `double` Radius of the sphere.
+ */
+ TransitionMatrix(
+ int _lm, double _vk, double _exri, std::complex<double> **_rmi,
+ std::complex<double> **_rei, double _sphere_radius
+ );
+
+ /*! \brief Transition Matrix instance constructor for a cluster of spheres.
+ *
+ * \param _nlemt: `int` Size of the matrix (2 * LE * (LE + 2)).
+ * \param _lm: `int` Maximum field expansion order.
+ * \param _vk: `double`
+ * \param _exri: `double`
+ * \param _am0m: Matrix of complex.
+ */
+ TransitionMatrix(int _nlemt, int _lm, double _vk, double _exri, std::complex<double> **am0m);
+
+ /*! \brief Transition Matrix instance destroyer.
+ */
+ ~TransitionMatrix();
+
+ /*! \brief Write the Transition Matrix to a binary file.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
+ * (default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode="LEGACY");
+};
+
+#endif
diff --git a/src/libnptm/TransitionMatrix.cpp b/src/libnptm/TransitionMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9a08cf77e890f57607ee800a7cdfe82846a3d10
--- /dev/null
+++ b/src/libnptm/TransitionMatrix.cpp
@@ -0,0 +1,104 @@
+/*! \file TransitionMatrix.cpp
+ *
+ * \brief Implementation of the Transition Matrix structure.
+ */
+#include <complex>
+#include <exception>
+#include <fstream>
+
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#include "../include/TransitionMatrix.h"
+#endif
+
+using namespace std;
+
+TransitionMatrix::~TransitionMatrix() {
+ if (elements != NULL) delete[] elements;
+ if (shape != NULL) delete[] shape;
+}
+
+TransitionMatrix::TransitionMatrix(
+ int _lm, double _vk, double _exri, complex<double> **_rmi,
+ complex<double> **_rei, double _sphere_radius
+) {
+ is = 1111;
+ shape = new int[2];
+ shape[0] = _lm;
+ shape[1] = 2;
+ l_max = _lm;
+ vk = _vk;
+ exri = _exri;
+ sphere_radius = _sphere_radius;
+ elements = new complex<double>[2 * l_max]();
+ for (int ei = 0; ei < l_max; ei++) {
+ elements[2 * ei] = -1.0 / _rmi[ei][0];
+ elements[2 * ei + 1] = -1.0 / _rei[ei][0];
+ }
+}
+
+TransitionMatrix::TransitionMatrix(
+ int _nlemt, int _lm, double _vk, double _exri,
+ std::complex<double> **am0m
+) {
+ is = 1;
+ shape = new int[2];
+ shape[0] = _nlemt;
+ shape[1] = _nlemt;
+ l_max = _lm;
+ vk = _vk;
+ exri = _exri;
+ sphere_radius = 0.0;
+ elements = new complex<double>[_nlemt * _nlemt]();
+ for (int ei = 0; ei < _nlemt; ei++) {
+ for (int ej = 0; ej < _nlemt; ej++) elements[_nlemt * ei + ej] = am0m[ei][ej];
+ }
+}
+
+void TransitionMatrix::write_binary(string file_name, string mode) {
+ if (mode.compare("LEGACY") == 0) {
+ write_legacy(file_name);
+ } else if (mode.compare("HDF5") == 0) {
+ write_hdf5(file_name);
+ } else {
+ string message = "Unknown format mode: \"" + mode + "\"";
+ throw UnrecognizedFormatException(message);
+ }
+}
+
+void TransitionMatrix::write_hdf5(string file_name) {
+ // TODO: needs implementation.
+ return;
+}
+
+void TransitionMatrix::write_legacy(string file_name) {
+ fstream ttms;
+ if (is == 1111 || is == 1) {
+ ttms.open(file_name, ios::binary | ios::out);
+ if (ttms.is_open()) {
+ ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&l_max), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ }
+ } else {
+ string message = "Unrecognized matrix data.";
+ throw UnrecognizedFormatException(message);
+ }
+ if (ttms.is_open()) {
+ int num_elements = shape[0] * shape[1];
+ for (int ei = 0; ei < num_elements; ei++) {
+ complex<double> element1 = elements[ei];
+ double vreal, vimag;
+ vreal = element1.real();
+ vimag = element1.imag();
+ ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+ }
+ if (is == 1111) {
+ ttms.write(reinterpret_cast<char *>(&sphere_radius), sizeof(double));
+ }
+ ttms.close();
+ } else { // Failed to open output file. Should never happen.
+ printf("ERROR: could not open Transition Matrix file for writing.\n");
+ }
+}
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index 384a85b9ede135f920d78585c0f924495e9d13f9..a6410192e143414085bc1fd1a35073dde2dc78e7 100644
--- a/src/sphere/Makefile
+++ b/src/sphere/Makefile
@@ -17,8 +17,8 @@ edfb: edfb.o
sph: sph.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/sph $(BUILDDIR)/sph.o $(LDFLAGS)
-np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o
- $(CXX) $(CXXFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o $(CXXLDFLAGS)
+np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o $(BUILDDIR)/TransitionMatrix.o
+ $(CXX) $(CXXFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o $(BUILDDIR)/TransitionMatrix.o $(CXXLDFLAGS)
#$(BUILDDIR)/np_sphere.o:
# $(CXX) $(CXXFLAGS) -c np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index f3264ce867404ce120770ca43759ebf5592c172f..129b2f37a0378450faae5113694f84f83083e6be 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -2,10 +2,11 @@
*
* \brief Implementation of the single sphere calculation.
*/
+#include <complex>
#include <cstdio>
+#include <exception>
#include <fstream>
#include <string>
-#include <complex>
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
@@ -19,6 +20,10 @@
#include "../include/sph_subs.h"
#endif
+#ifndef INCLUDE_TRANSITIONMATRIX_H_
+#include "../include/TransitionMatrix.h"
+#endif
+
using namespace std;
/*! \brief C++ implementation of SPH
@@ -286,34 +291,9 @@ void sphere(string config_file, string data_file, string output_path) {
} // i132
if (sconf->idfc >= 0 and nsph == 1 and jxi == gconf->jwtm) {
// This is the condition that writes the transition matrix to output.
- int is = 1111;
- fstream ttms;
+ TransitionMatrix ttms(gconf->l_max, vk, exri, c1->rmi, c1->rei, sconf->radii_of_spheres[0]);
string ttms_name = output_path + "/c_TTMS";
- ttms.open(ttms_name.c_str(), ios::binary | ios::out);
- if (ttms.is_open()) {
- ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
- ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
- ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
- for (int lmi = 0; lmi < gconf->l_max; lmi++) {
- complex<double> element1 = -1.0 / c1->rmi[lmi][0];
- complex<double> element2 = -1.0 / c1->rei[lmi][0];
- double vreal, vimag;
- vreal = element1.real();
- vimag = element1.imag();
- ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- vreal = element2.real();
- vimag = element2.imag();
- ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
- ttms.close();
- } else { // Failed to open output file. Should never happen.
- printf("ERROR: could not open TTMS file.\n");
- tppoan.close();
- }
+ ttms.write_binary(ttms_name, "LEGACY");
}
double cs0 = 0.25 * vk * vk * vk / half_pi;
//printf("DEBUG: cs0 = %lE\n", cs0);