Define a TransitionMatrix class

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>

src/include/TransitionMatrix.h

+/*! \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

src/libnptm/TransitionMatrix.cpp

+/*! \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");
+  }
+}

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

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);