From 90ee07a04953a90c19a62de6cf4da496eccb39f8 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 13 Feb 2025 12:39:08 +0100
Subject: [PATCH] Implement InclusionOutputInfo class

---
 src/include/outputs.h   |  442 ++++++++-
 src/libnptm/outputs.cpp | 2037 ++++++++++++++++++++++++++++++++++++++-
 2 files changed, 2467 insertions(+), 12 deletions(-)

diff --git a/src/include/outputs.h b/src/include/outputs.h
index 1370745b..b508e79d 100644
--- a/src/include/outputs.h
+++ b/src/include/outputs.h
@@ -21,6 +21,7 @@
 #ifndef INCLUDE_OUTPUTS_H_
 #define INCLUDE_OUTPUTS_H_
 
+// >>> OUTPUT FOR CLUSTER <<<
 /*! \brief Class to collect output information for scattering from clusters.
  *
  * The results of the calculation can be saved in different formats.
@@ -54,12 +55,7 @@ protected:
    *
    * This function takes care of writing the output using the legacy
    * formatted ASCII structure. If the output file does not exist, it
-   * is created. If it exists, the new content is appended to the old
-   * one, without overwriting. This is useful for a simpler handling
-   * of multi-process output combination, provided that only one process
-   * writes to the output file at any given time. Simultaneous write
-   * operations by different processes is not supported by the legacy
-   * ASCII format.
+   * is created. If it exists, the new content overwritten.
    *
    * \param output: `const string &` Path to the output to be written.
    * \return result: `int` Exit code (0 if successful).
@@ -500,7 +496,7 @@ public:
 
   /*! \brief Insert in the current output data the data of another block.
    *
-   * \param rhs: `const ClusterOutputData &` Reference to the source data block.
+   * \param rhs: `const ClusterOutputInfo &` Reference to the source data block.
    * \return result: `int` Exit code (0 if successful).
    */
   int insert(const ClusterOutputInfo &rhs);
@@ -541,4 +537,436 @@ public:
   int mpisend(const mixMPI *mpidata);
 #endif // MPI_VERSION
 };
+// >>> END OF OUTPUT FOR CLUSTER <<<
+
+// >>> OUTPUT FOR INCLUSION <<<
+/*! \brief Class to collect output information for scattering from particle with inclusions.
+ *
+ * The results of the calculation can be saved in different formats.
+ * It is therefore convenient to have a proper memory structure that
+ * allows for storing the results and flushing them in any of the
+ * permitted formats with just one operation. The purpose of the
+ * `InclusionOutputInfo` class is to provide a wrapper for the output
+ * of the particle with inclusions scattering solver.
+ */
+class InclusionOutputInfo {
+protected:
+  //! \brief Number of incident azimuth calculations.
+  int _num_theta;
+  //! \brief Number of scattered azimuth calculations.
+  int _num_thetas;
+  //! \brief Number of incident elevation calculations.
+  int _num_phi;
+  //! \brief Number of scattered elevation calculations.
+  int _num_phis;
+  //! \brief ID of the first computed wavelength
+  int _first_xi;
+  
+  /*! \brief Write the output to a HDF5 file.
+   *
+   * \param file_name: `const string &` Path to the output to be written.
+   * \return result: `int` Exit code (0 if successful).
+   */
+  int write_hdf5(const std::string &file_name);
+  
+  /*! \brief Write the output to a legacy text file.
+   *
+   * This function takes care of writing the output using the legacy
+   * formatted ASCII structure. If the output file does not exist, it
+   * is created. If it exists, the new content is overwritten.
+   *
+   * \param output: `const string &` Path to the output to be written.
+   * \return result: `int` Exit code (0 if successful).
+   */
+  int write_legacy(const std::string &output);
+  
+public:
+  //! \brief Read-only view on the ID of the first scale
+  const int &first_xi = _first_xi;
+  //! \brief Number of spheres in the aggregate.
+  int nsph;
+  //! \brief Maximum internal field expansion order.
+  int li;
+  //! \brief Maximum external field expansion order.
+  int le;
+  //! \brief Maximum field expansion order.
+  int lm;
+  //! \brief Maximum coefficient matrix dimension.
+  np_int mxndm;
+  //! \brief Incident polarization flag.
+  int inpol;
+  //! \brief Number of points for transition layer integration.
+  int npnt;
+  //! \brief Number of points for non-transition layer integration.
+  int npntts;
+  //! \brief Flag for intensity.
+  int iavm;
+  //! \brief Flag for reference to meridional plane.
+  int isam;
+  //! \brief Flag for dielectric function definition.
+  int idfc;
+  //! \brief Vector of spherical components X Cartesian coordinates.
+  double *vec_x_coords;
+  //! \brief Vector of spherical components Y Cartesian coordinates.
+  double *vec_y_coords;
+  //! \brief Vector of spherical components Z Cartesian coordinates.
+  double *vec_z_coords;
+  //! \brief First incident radiation azimuth angle.
+  double th;
+  //! \brief Incident radiation azimuth angle step.
+  double thstp;
+  //! \brief Last incident radiation azimuth angle.
+  double thlst;
+  //! \brief First scattered radiation azimuth angle.
+  double ths;
+  //! \brief Scattered radiation azimuth angle step.
+  double thsstp;
+  //! \brief Last scattered radiation azimuth angle.
+  double thslst;
+  //! \brief First incident radiation elevation angle.
+  double ph;
+  //! \brief Incident radiation elevation angle step.
+  double phstp;
+  //! \brief Last incident radiation elevation angle.
+  double phlst;
+  //! \brief First scattered radiation elevation angle.
+  double phs;
+  //! \brief Scattered radiation elevation angle step.
+  double phsstp;
+  //! \brief Last scattered radiation elevation angle.
+  double phslst;
+  //! \brief Number of directions to be explicitly solved.
+  int ndirs;
+  //! \brief Refractive index of external medium.
+  double exri;
+  //! \brief Number of scales (wavelengths)
+  int nxi;
+  //! \brief Number of scales handled by the current process.
+  int xi_block_size;
+  //! \brief Index of the wavelength for T-matrix output.
+  int jwtm;
+  //! \brief Vector of scale (wavelength) indices.
+  int *vec_jxi;
+  //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).
+  short *vec_ier;
+  //! \brief Vector of vacuum wave numbers.
+  double *vec_vk;
+  //! \brief Vector of computed scales.
+  double *vec_xi;
+  //! \brief Number of sphere configurations.
+  int configurations;
+  //! \brief Vector of sphere sizes (all configurations for every scale).
+  double *vec_sphere_sizes;
+  //! \brief Vector of sphere refractive indices  (all configurations for every scale).
+  dcomplex *vec_sphere_ref_indices;
+  //! \brief Vector of particle scattering cross-sections (parallel polarization).
+  double *vec_scs1;
+  //! \brief Vector of particle scattering cross-sections (perpendicular polarization).
+  double *vec_scs2;
+  //! \brief Vector of particle absorption cross-sections (parallel polarization).
+  double *vec_abs1;
+  //! \brief Vector of particle absorption cross-sections (perpendicular polarization).
+  double *vec_abs2;
+  //! \brief Vector of particle extinction cross-sections (parallel polarization).
+  double *vec_exs1;
+  //! \brief Vector of particle extinction cross-sections (perpendicular polarization).
+  double *vec_exs2;
+  //! \brief Vector of particle albedos (parallel polarization).
+  double *vec_albeds1;
+  //! \brief Vector of particle albedos (perpendicular polarization).
+  double *vec_albeds2;
+  //! \brief Vector of particle scattering-to-geometric cross-sections (parallel polarization).
+  double *vec_scsrt1;
+  //! \brief Vector of particle scattering-to-geometric cross-sections (perpendicular polarization).
+  double *vec_scsrt2;
+  //! \brief Vector of particle absorption-to-geometric cross-sections (parallel polarization).
+  double *vec_absrt1;
+  //! \brief Vector of particle absorption-to-geometric cross-sections (perpendicular polarization).
+  double *vec_absrt2;
+  //! \brief Vector of particle extinction-to-geometric cross-sections (parallel polarization).
+  double *vec_exsrt1;
+  //! \brief Vector of particle extinction-to-geometric cross-sections (perpendicular polarization).
+  double *vec_exsrt2;
+  //! \brief Vector of particle QSCHU (parallel polarization).
+  double *vec_qschu1;
+  //! \brief Vector of particle QSCHU (perpendicular polarization).
+  double *vec_qschu2;
+  //! \brief Vector of particle PSCHU (parallel polarization).
+  double *vec_pschu1;
+  //! \brief Vector of particle PSCHU (perpendicular polarization).
+  double *vec_pschu2;
+  //! \brief Vector of particle S0MAG (parallel polarization).
+  double *vec_s0mag1;
+  //! \brief Vector of particle S0MAG (perpendicular polarization).
+  double *vec_s0mag2;
+  //! \brief Vector of particle average asymmetry parameter (parallel polarization).
+  double *vec_cosav1;
+  //! \brief Vector of particle average asymmetry parameter (perpendicular polarization).
+  double *vec_cosav2;
+  //! \brief Vector of particle average radiation pressure force (N - parallel polarization).
+  double *vec_raprs1;
+  //! \brief Vector of particle average radiation pressure force (N - perpendicular polarization).
+  double *vec_raprs2;
+  //! \brief Vector of particle average radiation force along incidence direction (N - parallel polarization).
+  double *vec_fk1;
+  //! \brief Vector of particle average radiation force along incidence direction (N - perpendicular polarization).
+  double *vec_fk2;
+  //! \brief Vector of forward scattering amplitudes for polarization parallel to incidence (one per scale).
+  dcomplex *vec_fsas11;
+  //! \brief Vector of forward scattering amplitudes for polarization perpendicular to incidence (one per scale).
+  dcomplex *vec_fsas21;
+  //! \brief Vector of forward scattering amplitudes for polarization parallel to incidence (one per scale).
+  dcomplex *vec_fsas22;
+  //! \brief Vector of forward scattering amplitudes for polarization perpendicular to incidence (one per scale).
+  dcomplex *vec_fsas12;
+  //! \brief Vector of incidence azimuth directions (one per incidence azimuth).
+  double *vec_dir_tidg;
+  //! \brief Vector of incidence elevation directions (one per incidence elevation).
+  double *vec_dir_pidg;
+  //! \brief Vector of scattering azimuth directions (one per scattering azimuth).
+  double *vec_dir_tsdg;
+  //! \brief Vector of scattering elevation directions (one per scattering elevation).
+  double *vec_dir_psdg;
+  //! \brief Vector of scattering angles (one per direction).
+  double *vec_dir_scand;
+  //! \brief Control parameter for incidence plane referred to meridional plane (one per direction).
+  double *vec_dir_cfmp;
+  //! \brief Control parameter for scattering plane referred to meridional plane (one per direction).
+  double *vec_dir_sfmp;
+  //! \brief Control parameter for incidence plane referred to scattering plane (one per direction).
+  double *vec_dir_cfsp;
+  //! \brief Control parameter for scattering plane referred to scattering plane (one per direction).
+  double *vec_dir_sfsp;
+  //! \brief Components of the unitary vector perpendicular to incidence plane (three per direction).
+  double *vec_dir_un;
+  //! \brief Components of the unitary vector perpendicular to scattering plane (three per direction).
+  double *vec_dir_uns;
+  //! \brief Vector of particle differential scattering cross-sections (parallel polarization).
+  double *vec_dir_scs1;
+  //! \brief Vector of particle differential scattering cross-sections (perpendicular polarization).
+  double *vec_dir_scs2;
+  //! \brief Vector of particle differential absorption cross-sections (parallel polarization).
+  double *vec_dir_abs1;
+  //! \brief Vector of particle differential absorption cross-sections (perpendicular polarization).
+  double *vec_dir_abs2;
+  //! \brief Vector of particle differential extinction cross-sections (parallel polarization).
+  double *vec_dir_exs1;
+  //! \brief Vector of particle differential extinction cross-sections (perpendicular polarization).
+  double *vec_dir_exs2;
+  //! \brief Vector of particle differential albedos (parallel polarization).
+  double *vec_dir_albeds1;
+  //! \brief Vector of particle differential albedos (perpendicular polarization).
+  double *vec_dir_albeds2;
+  //! \brief Vector of particle differential scattering-to-geometric cross-sections (parallel polarization).
+  double *vec_dir_scsrt1;
+  //! \brief Vector of particle differential scattering-to-geometric cross-sections (perpendicular polarization).
+  double *vec_dir_scsrt2;
+  //! \brief Vector of particle differential absorption-to-geometric cross-sections (parallel polarization).
+  double *vec_dir_absrt1;
+  //! \brief Vector of particle differential absorption-to-geometric cross-sections (perpendicular polarization).
+  double *vec_dir_absrt2;
+  //! \brief Vector of particle differential extinction-to-geometric cross-sections (parallel polarization).
+  double *vec_dir_exsrt1;
+  //! \brief Vector of particle differential extinction-to-geometric cross-sections (perpendicular polarization).
+  double *vec_dir_exsrt2;
+  //! \brief Vector of particle differential forward scattering amplitude with polarization parallel to parallel incidence field.
+  dcomplex *vec_dir_fsas11;
+  //! \brief Vector of particle differential forward scattering amplitude with polarization perpendicular to the parallel incidence field.
+  dcomplex *vec_dir_fsas21;
+  //! \brief Vector of particle differential forward scattering amplitude with polarization perpendicular to perpendicular incidence field.
+  dcomplex *vec_dir_fsas12;
+  //! \brief Vector of particle differential forward scattering amplitude with polarization parallel the perpendicular incidence field.
+  dcomplex *vec_dir_fsas22;
+  //! \brief Vector of particle differential scattering amplitude with polarization parallel to parallel incidence field.
+  dcomplex *vec_dir_sas11;
+  //! \brief Vector of particle differential scattering amplitude with polarization perpendicular to the parallel incidence field.
+  dcomplex *vec_dir_sas21;
+  //! \brief Vector of particle differential scattering amplitude with polarization perpendicular to perpendicular incidence field.
+  dcomplex *vec_dir_sas12;
+  //! \brief Vector of particle differential scattering amplitude with polarization parallel the perpendicular incidence field.
+  dcomplex *vec_dir_sas22;
+  //! \brief Vector of differential particle QSCHU (parallel polarization).
+  double *vec_dir_qschu1;
+  //! \brief Vector of differential particle QSCHU (perpendicular polarization).
+  double *vec_dir_qschu2;
+  //! \brief Vector of differential particle PSCHU (parallel polarization).
+  double *vec_dir_pschu1;
+  //! \brief Vector of differential particle PSCHU (perpendicular polarization).
+  double *vec_dir_pschu2;
+  //! \brief Vector of particle differential S0MAG (parallel polarization).
+  double *vec_dir_s0mag1;
+  //! \brief Vector of particle differential S0MAG (perpendicular polarization).
+  double *vec_dir_s0mag2;
+  //! \brief Vector of differential particle asymmetry parameters (parallel polarization).
+  double *vec_dir_cosav1;
+  //! \brief Vector of differential particle asymmetry parameters (perpendicular polarization).
+  double *vec_dir_cosav2;
+  //! \brief Vector of differential particle radiation pressure forces (1).
+  double *vec_dir_rapr1;
+  //! \brief Vector of differential particle radiation pressure forces (1).
+  double *vec_dir_rapr2;
+  //! \brief Vector of differential radiation pressure force components along the polarization direction (parallel polarization).
+  double *vec_dir_fl1;
+  //! \brief Vector of differential radiation pressure force components along the polarization direction (perpendicular polarization).
+  double *vec_dir_fl2;
+  //! \brief Vector of differential radiation pressure force components perpendicular to the polarization direction (parallel polarization).
+  double *vec_dir_fr1;
+  //! \brief Vector of differential radiation pressure force components perpendicular to the polarization direction (perpendicular polarization).
+  double *vec_dir_fr2;
+  //! \brief Vector of differential radiation pressure force components along the incidence direction (parallel polarization).
+  double *vec_dir_fk1;
+  //! \brief Vector of differential radiation pressure force components along the incidence direction (perpendicular polarization).
+  double *vec_dir_fk2;
+  //! \brief Vector of differential radiation pressure force components along the X axis (parallel polarization).
+  double *vec_dir_fx1;
+  //! \brief Vector of differential radiation pressure force components along the X axis (perpendicular polarization).
+  double *vec_dir_fx2;
+  //! \brief Vector of differential radiation pressure force components along the Y axis (parallel polarization).
+  double *vec_dir_fy1;
+  //! \brief Vector of differential radiation pressure force components along the Y axis (perpendicular polarization).
+  double *vec_dir_fy2;
+  //! \brief Vector of differential radiation pressure force components along the Z axis (parallel polarization).
+  double *vec_dir_fz1;
+  //! \brief Vector of differential radiation pressure force components along the Z axis (perpendicular polarization).
+  double *vec_dir_fz2;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the polarization direction (parallel polarization).
+  double *vec_dir_tqel1;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the polarization direction (perpendicular polarization).
+  double *vec_dir_tqel2;
+  //! \brief Vector of differential extinction contribution to radiation torque components perpendicular to the polarization direction (parallel polarization).
+  double *vec_dir_tqer1;
+  //! \brief Vector of differential extinction contribution to radiation torque components perpendicular to the polarization direction (perpendicular polarization).
+  double *vec_dir_tqer2;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the incidence direction (parallel polarization).
+  double *vec_dir_tqek1;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the incidence direction (perpendicular polarization).
+  double *vec_dir_tqek2;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the X axis (parallel polarization).
+  double *vec_dir_tqex1;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the X axis (perpendicular polarization).
+  double *vec_dir_tqex2;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the Y axis (parallel polarization).
+  double *vec_dir_tqey1;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the Y axis (perpendicular polarization).
+  double *vec_dir_tqey2;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the Z axis (parallel polarization).
+  double *vec_dir_tqez1;
+  //! \brief Vector of differential extinction contribution to radiation torque components along the Z axis (perpendicular polarization).
+  double *vec_dir_tqez2;
+  //! \brief Vector of differential scattering contribution to radiation torque components along the polarization direction (parallel polarization).
+  double *vec_dir_tqsl1;
+  //! \brief Vector of differential scattering contribution to radiation torque components along the polarization direction (perpendicular polarization).
+  double *vec_dir_tqsl2;
+  //! \brief Vector of differential scattering contribution to radiation torque components perpendicular to the polarization direction (parallel polarization).
+  double *vec_dir_tqsr1;
+  //! \brief Vector of differential scattering contribution to radiation torque components perpendicular to the polarization direction (perpendicular polarization).
+  double *vec_dir_tqsr2;
+  //! \brief Vector of differential scattering contribution to radiation torque components along the incidence direction (parallel polarization).
+  double *vec_dir_tqsk1;
+  //! \brief Vector of differential scattering contribution to radiation torque components along the incidence direction (perpendicular polarization).
+  double *vec_dir_tqsk2;
+  //! \brief Vector of differential scattering contribution to radiation torque components along X axis (parallel polarization).
+  double *vec_dir_tqsx1;
+  //! \brief Vector of differential scattering contribution to radiation torque components along X axis (perpendicular polarization).
+  double *vec_dir_tqsx2;
+  //! \brief Vector of differential scattering contribution to radiation torque components along Y axis (parallel polarization).
+  double *vec_dir_tqsy1;
+  //! \brief Vector of differential scattering contribution to radiation torque components along Y axis (perpendicular polarization).
+  double *vec_dir_tqsy2;
+  //! \brief Vector of differential scattering contribution to radiation torque components along Z axis (parallel polarization).
+  double *vec_dir_tqsz1;
+  //! \brief Vector of differential scattering contribution to radiation torque components along Z axis (perpendicular polarization).
+  double *vec_dir_tqsz2;
+  //! \brief Vector of cluster Mueller transormation matrices referred to meridional plane (16 per direction per scale).
+  double *vec_dir_mull;
+  //! \brief Vector of cluster Mueller transormation matrices referred to scattering plane (16 per direction per scale).
+  double *vec_dir_mulllr;
+  
+  /*! \brief `InclusionOutputInfo` default instance constructor.
+   *
+   * \param sc: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` instance.
+   * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.
+   * \param mpidata: `const mixMPI*` Pointer to a mixMPI instance.
+   * \param first_xi: `int` Index of the first scale in output (optional, default is 1).
+   * \param xi_length: `int` Number of scales tobe included in output (optional, default is 0, meaning all).
+   */   
+  InclusionOutputInfo(
+    ScattererConfiguration *sc, GeometryConfiguration *gc,
+    const mixMPI *mpidata, int first_xi = 1, int xi_length = 0
+  );
+
+  /*! \brief `InclusionOutputInfo` constructor from HDF5 input.
+   *
+   * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.
+   */   
+  InclusionOutputInfo(const std::string &hdf5_file_name);
+
+  /*! \brief `InclusionOutputInfo` instance destroyer.
+   */
+  ~InclusionOutputInfo();
+
+  /*! \brief Estimate the size of the structure that would be built for given input.
+   *
+   * \param sc: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` instance.
+   * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.
+   * \param first_xi: `int` Index of the first scale in output (optional, default is 0).
+   * \param xi_length: `int` Number of scales tobe included in output (optional, default is all).
+   * \return size: `long` Estimated instance size in bytes.
+   */
+  static long compute_size(
+    ScattererConfiguration *sc, GeometryConfiguration *gc,
+    int first_xi = 0, int xi_length = 0
+  );
+  
+  /*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
+   *
+   * \return size: `long` Estimated instance size in bytes.
+   */
+  long compute_size();
+
+  /*! \brief Insert in the current output data the data of another block.
+   *
+   * \param rhs: `const InclusionOutputInfo &` Reference to the source data block.
+   * \return result: `int` Exit code (0 if successful).
+   */
+  int insert(const InclusionOutputInfo &rhs);
+
+  /*! \brief Write the output to a file.
+   *
+   * \param output: `const string &` Path to the output to be written.
+   * \param format: `const string &` Output format (one of LEGACY or HDF5).
+   * \return result: `int` Exit code (0 if successful).
+   */
+  int write(const std::string &output, const std::string &format);
+
+#ifdef MPI_VERSION
+  /*! \brief Receive output data from worker processes.
+   *
+   * This function is invoked by the MPI rank-0 process to fetch the
+   * output data produced by higher rank processes. When calling this
+   * function, process 0 halts until a valid data chunk is transmitted
+   * by the queried process.
+   *
+   * \param mpidata: `const mixMPI*` Pointer to a `mixMPI` instance.
+   * \param pid: `int` Rank of the process that is transmitting data.
+   * \return result: `int` An exit code (0 for success).
+   */
+  int mpireceive(const mixMPI* mpidata, int pid);
+
+  /*! \brief Send output data to process 0.
+   *
+   * This function is invoked by non-zero ranked MPI processes when
+   * they are ready to send back the output data. When a process meets
+   * this function call, it halts until MPI process 0 asks for the
+   * data transmission.
+   *
+   * \param mpidata: `const mixMPI*` Pointer to a `mixMPI` instance.
+   * \param pid: `int` Rank of the process that is transmitting data.
+   * \return result: `int` An exit code (0 for success).
+   */
+  int mpisend(const mixMPI *mpidata);
+#endif // MPI_VERSION
+};
+// >>> END OF OUTPUT FOR INCLUSION <<<
+
 #endif // INCLUDE_OUTPUTS_H_
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 1f50f1a0..e121e95f 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -56,7 +56,7 @@
 
 using namespace std;
 
-// >> ClusterOutputInfo CLASS IMPLEMENTATION <<
+// >>> ClusterOutputInfo CLASS IMPLEMENTATION <<<
 ClusterOutputInfo::ClusterOutputInfo(
   ScattererConfiguration *sc, GeometryConfiguration *gc,
   const mixMPI *mpidata, int first_xi, int xi_length
@@ -873,7 +873,7 @@ long ClusterOutputInfo::compute_size(
 ) {
   long result = sizeof(np_int);
   result += 21 * sizeof(int);
-  result += 56 * sizeof(double);
+  result += 14 * sizeof(double);
   result += 121 * sizeof(long);
   int _nsph = gc->number_of_spheres;
   double _th = gc->in_theta_start;
@@ -928,7 +928,7 @@ long ClusterOutputInfo::compute_size(
 long ClusterOutputInfo::compute_size() {
   long result = sizeof(np_int);
   result += 21 * sizeof(int);
-  result += 56 * sizeof(double);
+  result += 14 * sizeof(double);
   result += 121 * sizeof(long);
   result += 3 * nsph * sizeof(double); // sphere coordinate vectors
   result += xi_block_size * sizeof(int); // scale index vector
@@ -2739,8 +2739,2035 @@ int ClusterOutputInfo::mpisend(const mixMPI *mpidata) {
     MPI_Send(vec_dir_mulc, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
     MPI_Send(vec_dir_mulclr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
   }
-  result = 0;
   return result;
 }
 #endif //MPI_VERSION
-// >> END OF ClusterOutputInfo CLASS IMPLEMENTATION <<
+// >>> END OF ClusterOutputInfo CLASS IMPLEMENTATION <<<
+
+// >>> InclusionOutputInfo CLASS IMPLEMENTATION <<<
+InclusionOutputInfo::InclusionOutputInfo(
+  ScattererConfiguration *sc, GeometryConfiguration *gc,
+  const mixMPI *mpidata, int first_xi, int xi_length
+) {
+  nsph = gc->number_of_spheres;
+  li = gc->li;
+  le = gc->le;
+  lm = gc->l_max;
+  mxndm = gc->mxndm;
+  inpol = gc->in_pol;
+  npnt = gc->npnt;
+  npntts = gc->npntts;
+  iavm = gc->iavm;
+  isam = gc->isam;
+  jwtm = gc->jwtm;
+  // Get spherical constituent coordinates
+  vec_x_coords = new double[nsph];
+  vec_y_coords = new double[nsph];
+  vec_z_coords = new double[nsph];
+  for (int nsi = 0; nsi < nsph; nsi++) {
+    vec_x_coords[nsi] = gc->get_sph_x(nsi);
+    vec_y_coords[nsi] = gc->get_sph_y(nsi);
+    vec_z_coords[nsi] = gc->get_sph_z(nsi);
+  }
+  // Get directional information
+  th = gc->in_theta_start;
+  thstp = gc->in_theta_step;
+  thlst = gc->in_theta_end;
+  ths = gc->sc_theta_start;
+  thsstp = gc->sc_theta_step;
+  thslst = gc->sc_theta_end;
+  _num_theta = (thstp == 0.0) ? 1 : 1 + int((thlst - th) / thstp);
+  _num_thetas = (thsstp == 0.0) ? 1 : 1 + int((thslst - ths) / thsstp);
+  ph = gc->in_phi_start;
+  phstp = gc->in_phi_step;
+  phlst = gc->in_phi_end;
+  phs = gc->sc_phi_start;
+  phsstp = gc->sc_phi_step;
+  phslst = gc->sc_phi_end;
+  _num_phi = (phstp == 0.0) ? 1 : 1 + int((phlst - ph) / phstp);
+  _num_phis = (phsstp == 0.0) ? 1 : 1 + int((phslst - phs) / phsstp);
+  ndirs = _num_theta * _num_thetas * _num_phi * _num_phis;
+  // Get scattering problem configuration
+  double exdc = sc->exdc;
+  exri = sqrt(exdc);
+  idfc = sc->idfc;
+  nxi = sc->number_of_scales;
+  _first_xi = first_xi;
+  xi_block_size = (xi_length == 0) ? nxi - (_first_xi - 1) : xi_length;
+  vec_jxi = new int[xi_block_size]();
+  vec_ier = new short[xi_block_size]();
+  vec_vk = new double[xi_block_size]();
+  vec_xi = new double[xi_block_size]();
+  configurations = sc->configurations;
+  vec_sphere_sizes = new double[xi_block_size * (configurations + 1)]();
+  vec_sphere_ref_indices = new dcomplex[xi_block_size * (configurations + 1)]();
+  vec_scs1 = new double[xi_block_size]();
+  vec_scs2 = new double[xi_block_size]();
+  vec_abs1 = new double[xi_block_size]();
+  vec_abs2 = new double[xi_block_size]();
+  vec_exs1 = new double[xi_block_size]();
+  vec_exs2 = new double[xi_block_size]();
+  vec_albeds1 = new double[xi_block_size]();
+  vec_albeds2 = new double[xi_block_size]();
+  vec_scsrt1 = new double[xi_block_size]();
+  vec_scsrt2 = new double[xi_block_size]();
+  vec_absrt1 = new double[xi_block_size]();
+  vec_absrt2 = new double[xi_block_size]();
+  vec_exsrt1 = new double[xi_block_size]();
+  vec_exsrt2 = new double[xi_block_size]();
+  vec_qschu1 = new double[xi_block_size]();
+  vec_qschu2 = new double[xi_block_size]();
+  vec_pschu1 = new double[xi_block_size]();
+  vec_pschu2 = new double[xi_block_size]();
+  vec_s0mag1 = new double[xi_block_size]();
+  vec_s0mag2 = new double[xi_block_size]();
+  vec_cosav1 = new double[xi_block_size]();
+  vec_cosav2 = new double[xi_block_size]();
+  vec_raprs1 = new double[xi_block_size]();
+  vec_raprs2 = new double[xi_block_size]();
+  vec_fk1 = new double[xi_block_size]();
+  vec_fk2 = new double[xi_block_size]();
+  vec_fsas11 = new dcomplex[xi_block_size]();
+  vec_fsas21 = new dcomplex[xi_block_size]();
+  vec_fsas22 = new dcomplex[xi_block_size]();
+  vec_fsas12 = new dcomplex[xi_block_size]();
+  vec_dir_tidg = new double[_num_theta];
+  vec_dir_pidg = new double[_num_phi];
+  vec_dir_tsdg = new double[_num_thetas];
+  vec_dir_psdg = new double[_num_phis];
+  // Initialize directions (they are scale-independent)
+  double cti = th, cpi = ph, cts = ths, cps = phs;
+  for (int di = 0; di < _num_theta; di++) {
+    vec_dir_tidg[di] = cti;
+    cti += thstp;
+  }
+  for (int di = 0; di < _num_thetas; di++) {
+    vec_dir_tsdg[di] = cts;
+    cts += thsstp;
+  }
+  for (int di = 0; di < _num_phi; di++) {
+    vec_dir_pidg[di] = cpi;
+    cpi += phstp;
+  }
+  for (int di = 0; di < _num_phis; di++) {
+    vec_dir_psdg[di] = cps;
+    cps += phsstp;
+  }
+  vec_dir_scand = new double[ndirs]();
+  vec_dir_cfmp = new double[ndirs]();
+  vec_dir_sfmp = new double[ndirs]();
+  vec_dir_cfsp = new double[ndirs]();
+  vec_dir_sfsp = new double[ndirs]();
+  vec_dir_un = new double[3 * ndirs]();
+  vec_dir_uns = new double[3 * ndirs]();
+  vec_dir_scs1 = new double[ndirs * xi_block_size]();
+  vec_dir_scs2 = new double[ndirs * xi_block_size]();
+  vec_dir_abs1 = new double[ndirs * xi_block_size]();
+  vec_dir_abs2 = new double[ndirs * xi_block_size]();
+  vec_dir_exs1 = new double[ndirs * xi_block_size]();
+  vec_dir_exs2 = new double[ndirs * xi_block_size]();
+  vec_dir_albeds1 = new double[ndirs * xi_block_size]();
+  vec_dir_albeds2 = new double[ndirs * xi_block_size]();
+  vec_dir_scsrt1 = new double[ndirs * xi_block_size]();
+  vec_dir_scsrt2 = new double[ndirs * xi_block_size]();
+  vec_dir_absrt1 = new double[ndirs * xi_block_size]();
+  vec_dir_absrt2 = new double[ndirs * xi_block_size]();
+  vec_dir_exsrt1 = new double[ndirs * xi_block_size]();
+  vec_dir_exsrt2 = new double[ndirs * xi_block_size]();
+  vec_dir_fsas11 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_fsas21 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_fsas12 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_fsas22 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_sas11 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_sas21 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_sas12 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_sas22 = new dcomplex[ndirs * xi_block_size]();
+  vec_dir_qschu1 = new double[ndirs * xi_block_size]();
+  vec_dir_qschu2 = new double[ndirs * xi_block_size]();
+  vec_dir_pschu1 = new double[ndirs * xi_block_size]();
+  vec_dir_pschu2 = new double[ndirs * xi_block_size]();
+  vec_dir_s0mag1 = new double[ndirs * xi_block_size]();
+  vec_dir_s0mag2 = new double[ndirs * xi_block_size]();
+  vec_dir_cosav1 = new double[ndirs * xi_block_size]();
+  vec_dir_cosav2 = new double[ndirs * xi_block_size]();
+  vec_dir_rapr1 = new double[ndirs * xi_block_size]();
+  vec_dir_rapr2 = new double[ndirs * xi_block_size]();
+  vec_dir_fl1 = new double[ndirs * xi_block_size]();
+  vec_dir_fl2 = new double[ndirs * xi_block_size]();
+  vec_dir_fr1 = new double[ndirs * xi_block_size]();
+  vec_dir_fr2 = new double[ndirs * xi_block_size]();
+  vec_dir_fk1 = new double[ndirs * xi_block_size]();
+  vec_dir_fk2 = new double[ndirs * xi_block_size]();
+  vec_dir_fx1 = new double[ndirs * xi_block_size]();
+  vec_dir_fx2 = new double[ndirs * xi_block_size]();
+  vec_dir_fy1 = new double[ndirs * xi_block_size]();
+  vec_dir_fy2 = new double[ndirs * xi_block_size]();
+  vec_dir_fz1 = new double[ndirs * xi_block_size]();
+  vec_dir_fz2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqel1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqel2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqer1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqer2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqek1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqek2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqex1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqex2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqey1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqey2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqez1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqez2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsl1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsl2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsr1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsr2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsk1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsk2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsx1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsx2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsy1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsy2 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsz1 = new double[ndirs * xi_block_size]();
+  vec_dir_tqsz2 = new double[ndirs * xi_block_size]();
+  vec_dir_mull = new double[16 * ndirs * xi_block_size]();
+  vec_dir_mulllr = new double[16 * ndirs * xi_block_size]();
+}
+
+InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_file_name) {
+  unsigned int flags = H5F_ACC_RDONLY;
+  HDFFile *hdf_file = new HDFFile(hdf5_file_name, flags);
+  herr_t status = hdf_file->get_status();
+  string str_name, str_type;
+  if (status == 0) {
+    status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
+    status = hdf_file->read("LI", "INT32_(1)", &li);
+    status = hdf_file->read("LE", "INT32_(1)", &le);
+    status = hdf_file->read("LM", "INT32_(1)", &lm);
+    long tmp;
+    status = hdf_file->read("MXNDM", "INT64_(1)", &tmp);
+    mxndm = (np_int)tmp;
+    status = hdf_file->read("INPOL", "INT32_(1)", &inpol);
+    status = hdf_file->read("NPNT", "INT32_(1)", &npnt);
+    status = hdf_file->read("NPNTTS", "INT32_(1)", &npntts);
+    status = hdf_file->read("IAVM", "INT32_(1)", &iavm);
+    status = hdf_file->read("ISAM", "INT32_(1)", &isam);
+    status = hdf_file->read("JWTM", "INT32_(1)", &jwtm);
+    str_type = "FLOAT64_(" + to_string(nsph) + ")";
+    vec_x_coords = new double[nsph];
+    vec_y_coords = new double[nsph];
+    vec_z_coords = new double[nsph];
+    status = hdf_file->read("VEC_SPH_X", str_type, vec_x_coords);
+    status = hdf_file->read("VEC_SPH_Y", str_type, vec_y_coords);
+    status = hdf_file->read("VEC_SPH_Z", str_type, vec_z_coords);
+    status = hdf_file->read("TH_START", "FLOAT64_(1)", &th);
+    status = hdf_file->read("TH_STEP", "FLOAT64_(1)", &thstp);
+    status = hdf_file->read("TH_END", "FLOAT64_(1)", &thlst);
+    _num_theta = (thstp == 0.0) ? 1 : 1 + int((thlst - th) / thstp);
+    status = hdf_file->read("THS_START", "FLOAT64_(1)", &ths);
+    status = hdf_file->read("THS_STEP", "FLOAT64_(1)", &thsstp);
+    status = hdf_file->read("THS_END", "FLOAT64_(1)", &thslst);
+    _num_thetas = (thsstp == 0.0) ? 1 : 1 + int((thslst - ths) / thsstp);
+    status = hdf_file->read("PH_START", "FLOAT64_(1)", &ph);
+    status = hdf_file->read("PH_STEP", "FLOAT64_(1)", &phstp);
+    status = hdf_file->read("PH_END", "FLOAT64_(1)", &phlst);
+    _num_phi = (phstp == 0.0) ? 1 : 1 + int((phlst - ph) / phstp);
+    status = hdf_file->read("PHS_START", "FLOAT64_(1)", &phs);
+    status = hdf_file->read("PHS_STEP", "FLOAT64_(1)", &phsstp);
+    status = hdf_file->read("PHS_END", "FLOAT64_(1)", &phslst);
+    _num_phis = (phsstp == 0.0) ? 1 : 1 + int((phslst - phs) / phsstp);
+    ndirs = _num_theta * _num_thetas * _num_phi * _num_phis;
+    status = hdf_file->read("EXRI", "FLOAT64_(1)", &exri);
+    status = hdf_file->read("IDFC", "INT32_(1)", &idfc);
+    status = hdf_file->read("XI1", "INT32_(1)", &_first_xi);
+    status = hdf_file->read("NXI", "INT32_(1)", &xi_block_size);
+    nxi = (_first_xi == 1) ? xi_block_size : xi_block_size + _first_xi;
+    str_type = "INT32_(" + to_string(xi_block_size) + ")";
+    vec_jxi = new int[xi_block_size];
+    status = hdf_file->read("VEC_JXI", str_type, vec_jxi);
+    str_type = "INT16_(" + to_string(xi_block_size) + ")";
+    vec_ier = new short[xi_block_size];
+    status = hdf_file->read("VEC_IER", str_type, vec_ier);
+    str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+    vec_vk = new double[xi_block_size];
+    status = hdf_file->read("VEC_VK", str_type, vec_vk);
+    vec_xi = new double[xi_block_size];
+    status = hdf_file->read("VEC_VK", str_type, vec_xi);
+    status = hdf_file->read("NCONF", "INT32_(1)", &configurations);
+    str_type = "FLOAT64_(" + to_string(xi_block_size * (configurations + 1)) + ")";
+    vec_sphere_sizes = new double[xi_block_size * (configurations + 1)];
+    status = hdf_file->read("VEC_SPH_SIZES", str_type, vec_sphere_sizes);
+    str_type = "FLOAT64_(" + to_string(2 * xi_block_size * (configurations + 1)) + ")";
+    vec_sphere_ref_indices = new dcomplex[xi_block_size * (configurations + 1)];
+    status = hdf_file->read("VEC_SPH_REFRI", str_type, vec_sphere_ref_indices);
+    str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+    vec_scs1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_SCS1", str_type, vec_scs1);
+    vec_scs2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_SCS2", str_type, vec_scs2);
+    vec_abs1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ABS1", str_type, vec_abs1);
+    vec_abs2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ABS2", str_type, vec_abs2);
+    vec_exs1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_EXS1", str_type, vec_exs1);
+    vec_exs2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_EXS2", str_type, vec_exs2);
+    vec_albeds1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ALBEDS1", str_type, vec_albeds1);
+    vec_albeds2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ALBEDS2", str_type, vec_albeds2);
+    vec_scsrt1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_SCSRT1", str_type, vec_scsrt1);
+    vec_scsrt2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_SCSRT2", str_type, vec_scsrt2);
+    vec_absrt1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ABSRT1", str_type, vec_absrt1);
+    vec_absrt2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_ABSRT2", str_type, vec_absrt2);
+    vec_exsrt1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_EXSRT1", str_type, vec_exsrt1);
+    vec_exsrt2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_EXSRT2", str_type, vec_exsrt2);
+    vec_qschu1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_QSCHU1", str_type, vec_qschu1);
+    vec_qschu2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_QSCHU2", str_type, vec_qschu2);
+    vec_pschu1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_PSCHU1", str_type, vec_pschu1);
+    vec_pschu2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_PSCHU2", str_type, vec_pschu2);
+    vec_s0mag1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_S0MAG1", str_type, vec_s0mag1);
+    vec_s0mag2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_S0MAG2", str_type, vec_s0mag2);
+    vec_cosav1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_COSAV1", str_type, vec_cosav1);
+    vec_cosav2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_COSAV2", str_type, vec_cosav2);
+    vec_raprs1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_RAPRS1", str_type, vec_raprs1);
+    vec_raprs2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_RAPRS2", str_type, vec_raprs2);
+    vec_fk1 = new double[xi_block_size];
+    status = hdf_file->read("VEC_FK1", str_type, vec_fk1);
+    vec_fk2 = new double[xi_block_size];
+    status = hdf_file->read("VEC_FK2", str_type, vec_fk2);
+    str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+    vec_fsas11 = new dcomplex[xi_block_size];
+    status = hdf_file->read("VEC_FSAS11", str_type, vec_fsas11);
+    vec_fsas21 = new dcomplex[xi_block_size];
+    status = hdf_file->read("VEC_FSAS21", str_type, vec_fsas21);
+    vec_fsas22 = new dcomplex[xi_block_size];
+    status = hdf_file->read("VEC_FSAS22", str_type, vec_fsas22);
+    vec_fsas12 = new dcomplex[xi_block_size];
+    status = hdf_file->read("VEC_FSAS12", str_type, vec_fsas12);
+    str_type = "FLOAT64_(" + to_string(_num_theta) + ")";
+    vec_dir_tidg = new double[_num_theta];
+    status = hdf_file->read("VEC_DIR_TIDG", str_type, vec_dir_tidg);
+    str_type = "FLOAT64_(" + to_string(_num_phi) + ")";
+    vec_dir_pidg = new double[_num_phi];
+    status = hdf_file->read("VEC_DIR_PIDG", str_type, vec_dir_pidg);
+    str_type = "FLOAT64_(" + to_string(_num_thetas) + ")";
+    vec_dir_tsdg = new double[_num_thetas];
+    status = hdf_file->read("VEC_DIR_TSDG", str_type, vec_dir_tsdg);
+    str_type = "FLOAT64_(" + to_string(_num_phis) + ")";
+    vec_dir_psdg = new double[_num_phis];
+    status = hdf_file->read("VEC_DIR_PSDG", str_type, vec_dir_psdg);
+    str_type = "FLOAT64_(" + to_string(ndirs) + ")";
+    vec_dir_scand = new double[ndirs];
+    status = hdf_file->read("VEC_DIR_SCAND", str_type, vec_dir_scand);
+    vec_dir_cfmp = new double[ndirs];
+    status = hdf_file->read("VEC_DIR_CFMP", str_type, vec_dir_cfmp);
+    vec_dir_sfmp = new double[ndirs];
+    status = hdf_file->read("VEC_DIR_SFMP", str_type, vec_dir_sfmp);
+    vec_dir_cfsp = new double[ndirs];
+    status = hdf_file->read("VEC_DIR_CFSP", str_type, vec_dir_cfsp);
+    vec_dir_sfsp = new double[ndirs];
+    status = hdf_file->read("VEC_DIR_SFSP", str_type, vec_dir_sfsp);
+    str_type = "FLOAT64_(" + to_string(3 * ndirs) + ")";
+    vec_dir_un = new double[3 * ndirs];
+    status = hdf_file->read("VEC_DIR_UN", str_type, vec_dir_un);
+    vec_dir_uns = new double[3 * ndirs];
+    status = hdf_file->read("VEC_DIR_UNS", str_type, vec_dir_uns);
+    str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+    vec_dir_scs1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SCS1", str_type, vec_dir_scs1);
+    vec_dir_scs2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SCS2", str_type, vec_dir_scs2);
+    vec_dir_abs1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ABS1", str_type, vec_dir_abs1);
+    vec_dir_abs2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ABS2", str_type, vec_dir_abs2);
+    vec_dir_exs1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_EXS1", str_type, vec_dir_exs1);
+    vec_dir_exs2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_EXS2", str_type, vec_dir_exs2);
+    vec_dir_albeds1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ALBEDS1", str_type, vec_dir_albeds1);
+    vec_dir_albeds2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ALBEDS2", str_type, vec_dir_albeds2);
+    vec_dir_scsrt1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SCSRT1", str_type, vec_dir_scsrt1);
+    vec_dir_scsrt2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SCSRT2", str_type, vec_dir_scsrt2);
+    vec_dir_absrt1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ABSRT1", str_type, vec_dir_absrt1);
+    vec_dir_absrt2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_ABSRT2", str_type, vec_dir_absrt2);
+    vec_dir_exsrt1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_EXSRT1", str_type, vec_dir_exsrt1);
+    vec_dir_exsrt2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_EXSRT2", str_type, vec_dir_exsrt2);
+    str_type = "FLOAT64_(" + to_string(2 * ndirs * xi_block_size) + ")";
+    vec_dir_fsas11 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FSAS11", str_type, vec_dir_fsas11);
+    vec_dir_fsas21 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FSAS21", str_type, vec_dir_fsas21);
+    vec_dir_fsas12 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FSAS12", str_type, vec_dir_fsas12);
+    vec_dir_fsas22 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FSAS22", str_type, vec_dir_fsas22);
+    vec_dir_sas11 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SAS11", str_type, vec_dir_sas11);
+    vec_dir_sas21 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SAS21", str_type, vec_dir_sas21);
+    vec_dir_sas12 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SAS12", str_type, vec_dir_sas12);
+    vec_dir_sas22 = new dcomplex[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_SAS22", str_type, vec_dir_sas22);
+    str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+    vec_dir_qschu1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_QSCHU1", str_type, vec_dir_qschu1);
+    vec_dir_qschu2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_QSCHU2", str_type, vec_dir_qschu2);
+    vec_dir_pschu1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_PSCHU1", str_type, vec_dir_pschu1);
+    vec_dir_pschu2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_PSCHU2", str_type, vec_dir_pschu2);
+    vec_dir_s0mag1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_S0MAG1", str_type, vec_dir_s0mag1);
+    vec_dir_s0mag2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_S0MAG2", str_type, vec_dir_s0mag2);
+    vec_dir_cosav1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_COSAV1", str_type, vec_dir_cosav1);
+    vec_dir_cosav2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_COSAV2", str_type, vec_dir_cosav2);
+    vec_dir_rapr1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_RAPR1", str_type, vec_dir_rapr1);
+    vec_dir_rapr2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_RAPR2", str_type, vec_dir_rapr2);
+    vec_dir_fl1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FL1", str_type, vec_dir_fl1);
+    vec_dir_fl2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FL2", str_type, vec_dir_fl2);
+    vec_dir_fr1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FR1", str_type, vec_dir_fr1);
+    vec_dir_fr2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FR2", str_type, vec_dir_fr2);
+    vec_dir_fk1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FK1", str_type, vec_dir_fk1);
+    vec_dir_fk2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FK2", str_type, vec_dir_fk2);
+    vec_dir_fx1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FX1", str_type, vec_dir_fx1);
+    vec_dir_fx2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FX2", str_type, vec_dir_fx2);
+    vec_dir_fy1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FY1", str_type, vec_dir_fy1);
+    vec_dir_fy2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FY2", str_type, vec_dir_fy2);
+    vec_dir_fz1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FZ1", str_type, vec_dir_fz1);
+    vec_dir_fz2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_FZ2", str_type, vec_dir_fz2);
+    vec_dir_tqel1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEL1", str_type, vec_dir_tqel1);
+    vec_dir_tqel2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEL2", str_type, vec_dir_tqel2);
+    vec_dir_tqer1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQER1", str_type, vec_dir_tqer1);
+    vec_dir_tqer2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQER2", str_type, vec_dir_tqer2);
+    vec_dir_tqek1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEK1", str_type, vec_dir_tqek1);
+    vec_dir_tqek2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEK2", str_type, vec_dir_tqek2);
+    vec_dir_tqex1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEX1", str_type, vec_dir_tqex1);
+    vec_dir_tqex2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEX2", str_type, vec_dir_tqex2);
+    vec_dir_tqey1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEY1", str_type, vec_dir_tqey1);
+    vec_dir_tqey2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEY2", str_type, vec_dir_tqey2);
+    vec_dir_tqez1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEZ1", str_type, vec_dir_tqez1);
+    vec_dir_tqez2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQEZ2", str_type, vec_dir_tqez2);
+    vec_dir_tqsl1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSL1", str_type, vec_dir_tqsl1);
+    vec_dir_tqsl2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSL2", str_type, vec_dir_tqsl2);
+    vec_dir_tqsr1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSR1", str_type, vec_dir_tqsr1);
+    vec_dir_tqsr2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSR2", str_type, vec_dir_tqsr2);
+    vec_dir_tqsk1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSK1", str_type, vec_dir_tqsk1);
+    vec_dir_tqsk2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSK2", str_type, vec_dir_tqsk2);
+    vec_dir_tqsx1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSX1", str_type, vec_dir_tqsx1);
+    vec_dir_tqsx2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSX2", str_type, vec_dir_tqsx2);
+    vec_dir_tqsy1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSY1", str_type, vec_dir_tqsy1);
+    vec_dir_tqsy2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSY2", str_type, vec_dir_tqsy2);
+    vec_dir_tqsz1 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSZ1", str_type, vec_dir_tqsz1);
+    vec_dir_tqsz2 = new double[ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_TQSZ2", str_type, vec_dir_tqsz2);
+    str_type = "FLOAT64_(" + to_string(16 * ndirs * xi_block_size) + ")";
+    vec_dir_mull = new double[16 * ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_MULL", str_type, vec_dir_mull);
+    vec_dir_mulllr = new double[16 * ndirs * xi_block_size];
+    status = hdf_file->read("VEC_DIR_MULLLR", str_type, vec_dir_mulllr);
+    status = hdf_file->close();
+    delete hdf_file;
+  } else {
+    if (hdf_file != NULL) delete hdf_file;
+    UnrecognizedFormatException ex("Error: " + hdf5_file_name + " not recognized as a valid HDF5 file!");
+    throw ex;
+  }
+}
+
+InclusionOutputInfo::~InclusionOutputInfo() {
+  delete[] vec_x_coords;
+  delete[] vec_y_coords;
+  delete[] vec_z_coords;
+  delete[] vec_jxi;
+  delete[] vec_ier;
+  delete[] vec_vk;
+  delete[] vec_xi;
+  delete[] vec_sphere_sizes;
+  delete[] vec_sphere_ref_indices;
+  delete[] vec_scs1;
+  delete[] vec_scs2;
+  delete[] vec_abs1;
+  delete[] vec_abs2;
+  delete[] vec_exs1;
+  delete[] vec_exs2;
+  delete[] vec_albeds1;
+  delete[] vec_albeds2;
+  delete[] vec_scsrt1;
+  delete[] vec_scsrt2;
+  delete[] vec_absrt1;
+  delete[] vec_absrt2;
+  delete[] vec_exsrt1;
+  delete[] vec_exsrt2;
+  delete[] vec_qschu1;
+  delete[] vec_qschu2;
+  delete[] vec_pschu1;
+  delete[] vec_pschu2;
+  delete[] vec_s0mag1;
+  delete[] vec_s0mag2;
+  delete[] vec_cosav1;
+  delete[] vec_cosav2;
+  delete[] vec_raprs1;
+  delete[] vec_raprs2;
+  delete[] vec_fk1;
+  delete[] vec_fk2;
+  delete[] vec_fsas11;
+  delete[] vec_fsas21;
+  delete[] vec_fsas22;
+  delete[] vec_fsas12;
+  delete[] vec_dir_tidg;
+  delete[] vec_dir_pidg;
+  delete[] vec_dir_tsdg;
+  delete[] vec_dir_psdg;
+  delete[] vec_dir_scand;
+  delete[] vec_dir_cfmp;
+  delete[] vec_dir_sfmp;
+  delete[] vec_dir_cfsp;
+  delete[] vec_dir_sfsp;
+  delete[] vec_dir_un;
+  delete[] vec_dir_uns;
+  delete[] vec_dir_scs1;
+  delete[] vec_dir_scs2;
+  delete[] vec_dir_abs1;
+  delete[] vec_dir_abs2;
+  delete[] vec_dir_exs1;
+  delete[] vec_dir_exs2;
+  delete[] vec_dir_albeds1;
+  delete[] vec_dir_albeds2;
+  delete[] vec_dir_scsrt1;
+  delete[] vec_dir_scsrt2;
+  delete[] vec_dir_absrt1;
+  delete[] vec_dir_absrt2;
+  delete[] vec_dir_exsrt1;
+  delete[] vec_dir_exsrt2;
+  delete[] vec_dir_fsas11;
+  delete[] vec_dir_fsas21;
+  delete[] vec_dir_fsas12;
+  delete[] vec_dir_fsas22;
+  delete[] vec_dir_sas11;
+  delete[] vec_dir_sas21;
+  delete[] vec_dir_sas12;
+  delete[] vec_dir_sas22;
+  delete[] vec_dir_qschu1;
+  delete[] vec_dir_qschu2;
+  delete[] vec_dir_pschu1;
+  delete[] vec_dir_pschu2;
+  delete[] vec_dir_s0mag1;
+  delete[] vec_dir_s0mag2;
+  delete[] vec_dir_cosav1;
+  delete[] vec_dir_cosav2;
+  delete[] vec_dir_rapr1;
+  delete[] vec_dir_rapr2;
+  delete[] vec_dir_fl1;
+  delete[] vec_dir_fl2;
+  delete[] vec_dir_fr1;
+  delete[] vec_dir_fr2;
+  delete[] vec_dir_fk1;
+  delete[] vec_dir_fk2;
+  delete[] vec_dir_fx1;
+  delete[] vec_dir_fx2;
+  delete[] vec_dir_fy1;
+  delete[] vec_dir_fy2;
+  delete[] vec_dir_fz1;
+  delete[] vec_dir_fz2;
+  delete[] vec_dir_tqel1;
+  delete[] vec_dir_tqel2;
+  delete[] vec_dir_tqer1;
+  delete[] vec_dir_tqer2;
+  delete[] vec_dir_tqek1;
+  delete[] vec_dir_tqek2;
+  delete[] vec_dir_tqex1;
+  delete[] vec_dir_tqex2;
+  delete[] vec_dir_tqey1;
+  delete[] vec_dir_tqey2;
+  delete[] vec_dir_tqez1;
+  delete[] vec_dir_tqez2;
+  delete[] vec_dir_tqsl1;
+  delete[] vec_dir_tqsl2;
+  delete[] vec_dir_tqsr1;
+  delete[] vec_dir_tqsr2;
+  delete[] vec_dir_tqsk1;
+  delete[] vec_dir_tqsk2;
+  delete[] vec_dir_tqsx1;
+  delete[] vec_dir_tqsx2;
+  delete[] vec_dir_tqsy1;
+  delete[] vec_dir_tqsy2;
+  delete[] vec_dir_tqsz1;
+  delete[] vec_dir_tqsz2;
+  delete[] vec_dir_mull;
+  delete[] vec_dir_mulllr;
+}
+
+long InclusionOutputInfo::compute_size() {
+  long result = sizeof(np_int);
+  result += 21 * sizeof(int);
+  result += 13 * sizeof(double);
+  result += 120 * sizeof(long); // vector pointers
+  result += xi_block_size * (sizeof(int) + sizeof(short)); // vec_jxi, vec_ier
+  result += 3 * nsph * sizeof(double); // coordinate vectors
+  result += 29 * xi_block_size * sizeof(double); // wavelength dependent real values
+  result += 5 * xi_block_size * sizeof(dcomplex); // wavelength dependent complex values
+  result += 15 * ndirs * sizeof(double); // values depending only on directions
+  result += 92 * ndirs * xi_block_size * sizeof(double); // values depending on directions and wavelengths
+  result += 8 * ndirs * xi_block_size * sizeof(dcomplex); // values depending on directions and wavelengths
+  return result;
+}
+
+long InclusionOutputInfo::compute_size(
+  ScattererConfiguration *sc, GeometryConfiguration *gc,
+  int first_xi, int xi_length
+) {
+  long result = sizeof(np_int);
+  result += 21 * sizeof(int);
+  result += 13 * sizeof(double);
+  result += 120 * sizeof(long); // vector pointers
+  int _nsph = gc->number_of_spheres;
+  double _th = gc->in_theta_start;
+  double _thstp = gc->in_theta_step;
+  double _thlst = gc->in_theta_end;
+  double _ths = gc->sc_theta_start;
+  double _thsstp = gc->sc_theta_step;
+  double _thslst = gc->sc_theta_end;
+  int num_theta = 1 + int((_thlst - _th) / _thstp);
+  int num_thetas = 1 + int((_thslst - _ths) / _thsstp);
+  double _ph = gc->in_phi_start;
+  double _phstp = gc->in_phi_step;
+  double _phlst = gc->in_phi_end;
+  double _phs = gc->sc_phi_start;
+  double _phsstp = gc->sc_phi_step;
+  double _phslst = gc->sc_phi_end;
+  int num_phi = 1 + int((_phlst - _ph) / _phstp);
+  int num_phis = 1 + int((_phslst - _phs) / _phsstp);
+  int _ndirs = num_theta * num_thetas * num_phi * num_phis;
+  int _nxi = sc->number_of_scales;
+  int _xi_block_size = (xi_length == 0) ? _nxi : xi_length;
+  int _configurations = sc->configurations;
+  result += _xi_block_size * (sizeof(int) + sizeof(short)); // vec_jxi, vec_ier
+  result += 3 * _nsph * sizeof(double); // coordinate vectors
+  result += 29 * _xi_block_size * sizeof(double); // wavelength dependent real values
+  result += 5 * _xi_block_size * sizeof(dcomplex); // wavelength dependent complex values
+  result += 15 * _ndirs * sizeof(double); // values depending only on directions
+  result += 92 * _ndirs * _xi_block_size * sizeof(double); // values depending on directions and wavelengths
+  result += 8 * _ndirs * _xi_block_size * sizeof(dcomplex); // values depending on directions and wavelengths
+  return result;
+}
+
+int InclusionOutputInfo::insert(const InclusionOutputInfo &rhs) {
+  int result = 0;
+  result += (rhs.inpol == inpol) ? 0 : 1;
+  result += (rhs.iavm == iavm) ? 0 : 1;
+  result += (rhs.isam == isam) ? 0 : 1;
+  result += (rhs._num_theta == _num_theta) ? 0 : 1;
+  result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
+  result += (rhs._num_phi == _num_phi) ? 0 : 1;
+  result += (rhs._num_phis == _num_phis) ? 0 : 1;
+  result += (rhs.ndirs == ndirs) ? 0 : 1;
+  result += (rhs.exri == exri) ? 0 : 1;
+  result += (rhs.idfc == idfc) ? 0 : 1;
+  result += (rhs.configurations == configurations) ? 0 : 1;
+  if (result == 0) {
+    int offset, chunk_size, xi1;
+    xi1 = rhs._first_xi;
+    // Insert vectors whose size depends on wavelengths
+    offset = xi1 - _first_xi;
+    chunk_size = rhs.xi_block_size;
+    memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
+    memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
+    memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
+    memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
+    memcpy(vec_scs1 + offset, rhs.vec_scs1, chunk_size * sizeof(double));
+    memcpy(vec_scs2 + offset, rhs.vec_scs2, chunk_size * sizeof(double));
+    memcpy(vec_abs1 + offset, rhs.vec_abs1, chunk_size * sizeof(double));
+    memcpy(vec_abs2 + offset, rhs.vec_abs2, chunk_size * sizeof(double));
+    memcpy(vec_exs1 + offset, rhs.vec_exs1, chunk_size * sizeof(double));
+    memcpy(vec_exs2 + offset, rhs.vec_exs2, chunk_size * sizeof(double));
+    memcpy(vec_albeds1 + offset, rhs.vec_albeds1, chunk_size * sizeof(double));
+    memcpy(vec_albeds2 + offset, rhs.vec_albeds2, chunk_size * sizeof(double));
+    memcpy(vec_scsrt1 + offset, rhs.vec_scsrt1, chunk_size * sizeof(double));
+    memcpy(vec_scsrt2 + offset, rhs.vec_scsrt2, chunk_size * sizeof(double));
+    memcpy(vec_absrt1 + offset, rhs.vec_absrt1, chunk_size * sizeof(double));
+    memcpy(vec_absrt2 + offset, rhs.vec_absrt2, chunk_size * sizeof(double));
+    memcpy(vec_exsrt1 + offset, rhs.vec_exsrt1, chunk_size * sizeof(double));
+    memcpy(vec_exsrt2 + offset, rhs.vec_exsrt2, chunk_size * sizeof(double));
+    memcpy(vec_qschu1 + offset, rhs.vec_qschu1, chunk_size * sizeof(double));
+    memcpy(vec_qschu2 + offset, rhs.vec_qschu2, chunk_size * sizeof(double));
+    memcpy(vec_pschu1 + offset, rhs.vec_pschu1, chunk_size * sizeof(double));
+    memcpy(vec_pschu2 + offset, rhs.vec_pschu2, chunk_size * sizeof(double));
+    memcpy(vec_s0mag1 + offset, rhs.vec_s0mag1, chunk_size * sizeof(double));
+    memcpy(vec_s0mag2 + offset, rhs.vec_s0mag2, chunk_size * sizeof(double));
+    memcpy(vec_cosav1 + offset, rhs.vec_cosav1, chunk_size * sizeof(double));
+    memcpy(vec_cosav2 + offset, rhs.vec_cosav2, chunk_size * sizeof(double));
+    memcpy(vec_raprs1 + offset, rhs.vec_raprs1, chunk_size * sizeof(double));
+    memcpy(vec_raprs2 + offset, rhs.vec_raprs2, chunk_size * sizeof(double));
+    memcpy(vec_fk1 + offset, rhs.vec_fk1, chunk_size * sizeof(double));
+    memcpy(vec_fk2 + offset, rhs.vec_fk2, chunk_size * sizeof(double));
+    memcpy(vec_fsas11 + offset, rhs.vec_fsas11, chunk_size * sizeof(dcomplex));
+    memcpy(vec_fsas21 + offset, rhs.vec_fsas21, chunk_size * sizeof(dcomplex));
+    memcpy(vec_fsas22 + offset, rhs.vec_fsas22, chunk_size * sizeof(dcomplex));
+    memcpy(vec_fsas12 + offset, rhs.vec_fsas12, chunk_size * sizeof(dcomplex));
+    
+    // Insert vectors whose sizes depend on configurations
+    offset = (xi1 - _first_xi) * (configurations + 1);
+    chunk_size = rhs.xi_block_size * (configurations + 1);
+    memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
+    memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
+
+    // Insert vectors whose sizes depend on directions and wavelengths
+    offset = (xi1 - _first_xi) * ndirs;
+    chunk_size = rhs.xi_block_size * ndirs;
+    memcpy(vec_dir_scs1 + offset, rhs.vec_dir_scs1, chunk_size * sizeof(double));
+    memcpy(vec_dir_scs2 + offset, rhs.vec_dir_scs2, chunk_size * sizeof(double));
+    memcpy(vec_dir_abs1 + offset, rhs.vec_dir_abs1, chunk_size * sizeof(double));
+    memcpy(vec_dir_abs2 + offset, rhs.vec_dir_abs2, chunk_size * sizeof(double));
+    memcpy(vec_dir_exs1 + offset, rhs.vec_dir_exs1, chunk_size * sizeof(double));
+    memcpy(vec_dir_exs2 + offset, rhs.vec_dir_exs2, chunk_size * sizeof(double));
+    memcpy(vec_dir_albeds1 + offset, rhs.vec_dir_albeds1, chunk_size * sizeof(double));
+    memcpy(vec_dir_albeds2 + offset, rhs.vec_dir_albeds2, chunk_size * sizeof(double));
+    memcpy(vec_dir_scsrt1 + offset, rhs.vec_dir_scsrt1, chunk_size * sizeof(double));
+    memcpy(vec_dir_scsrt2 + offset, rhs.vec_dir_scsrt2, chunk_size * sizeof(double));
+    memcpy(vec_dir_absrt1 + offset, rhs.vec_dir_absrt1, chunk_size * sizeof(double));
+    memcpy(vec_dir_absrt2 + offset, rhs.vec_dir_absrt2, chunk_size * sizeof(double));
+    memcpy(vec_dir_exsrt1 + offset, rhs.vec_dir_exsrt1, chunk_size * sizeof(double));
+    memcpy(vec_dir_exsrt2 + offset, rhs.vec_dir_exsrt2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fsas11 + offset, rhs.vec_dir_fsas11, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_fsas21 + offset, rhs.vec_dir_fsas21, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_fsas12 + offset, rhs.vec_dir_fsas12, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_fsas22 + offset, rhs.vec_dir_fsas22, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
+    memcpy(vec_dir_qschu1 + offset, rhs.vec_dir_qschu1, chunk_size * sizeof(double));
+    memcpy(vec_dir_qschu2 + offset, rhs.vec_dir_qschu2, chunk_size * sizeof(double));
+    memcpy(vec_dir_pschu1 + offset, rhs.vec_dir_pschu1, chunk_size * sizeof(double));
+    memcpy(vec_dir_pschu2 + offset, rhs.vec_dir_pschu2, chunk_size * sizeof(double));
+    memcpy(vec_dir_s0mag1 + offset, rhs.vec_dir_s0mag1, chunk_size * sizeof(double));
+    memcpy(vec_dir_s0mag2 + offset, rhs.vec_dir_s0mag2, chunk_size * sizeof(double));
+    memcpy(vec_dir_cosav1 + offset, rhs.vec_dir_cosav1, chunk_size * sizeof(double));
+    memcpy(vec_dir_cosav2 + offset, rhs.vec_dir_cosav2, chunk_size * sizeof(double));
+    memcpy(vec_dir_rapr1 + offset, rhs.vec_dir_rapr1, chunk_size * sizeof(double));
+    memcpy(vec_dir_rapr2 + offset, rhs.vec_dir_rapr2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fl1 + offset, rhs.vec_dir_fl1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fl2 + offset, rhs.vec_dir_fl2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fr1 + offset, rhs.vec_dir_fr1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fr2 + offset, rhs.vec_dir_fr2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fk1 + offset, rhs.vec_dir_fk1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fk2 + offset, rhs.vec_dir_fk2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fx1 + offset, rhs.vec_dir_fx1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fx2 + offset, rhs.vec_dir_fx2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fy1 + offset, rhs.vec_dir_fy1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fy2 + offset, rhs.vec_dir_fy2, chunk_size * sizeof(double));
+    memcpy(vec_dir_fz1 + offset, rhs.vec_dir_fz1, chunk_size * sizeof(double));
+    memcpy(vec_dir_fz2 + offset, rhs.vec_dir_fz2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqel1 + offset, rhs.vec_dir_tqel1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqel2 + offset, rhs.vec_dir_tqel2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqer1 + offset, rhs.vec_dir_tqer1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqer2 + offset, rhs.vec_dir_tqer2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqek1 + offset, rhs.vec_dir_tqek1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqek2 + offset, rhs.vec_dir_tqek2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqex1 + offset, rhs.vec_dir_tqex1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqex2 + offset, rhs.vec_dir_tqex2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqey1 + offset, rhs.vec_dir_tqey1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqey2 + offset, rhs.vec_dir_tqey2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqez1 + offset, rhs.vec_dir_tqez1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqez2 + offset, rhs.vec_dir_tqez2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsl1 + offset, rhs.vec_dir_tqsl1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsl2 + offset, rhs.vec_dir_tqsl2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsr1 + offset, rhs.vec_dir_tqsr1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsr2 + offset, rhs.vec_dir_tqsr2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsk1 + offset, rhs.vec_dir_tqsk1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsk2 + offset, rhs.vec_dir_tqsk2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsx1 + offset, rhs.vec_dir_tqsx1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsx2 + offset, rhs.vec_dir_tqsx2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsy1 + offset, rhs.vec_dir_tqsy1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsy2 + offset, rhs.vec_dir_tqsy2, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsz1 + offset, rhs.vec_dir_tqsz1, chunk_size * sizeof(double));
+    memcpy(vec_dir_tqsz2 + offset, rhs.vec_dir_tqsz2, chunk_size * sizeof(double));
+    memcpy(vec_dir_mull + 16 * offset, rhs.vec_dir_mull, 16 * chunk_size * sizeof(double));
+    memcpy(vec_dir_mulllr + 16 * offset, rhs.vec_dir_mulllr, 16 * chunk_size * sizeof(double));
+  }
+  return result;
+}
+
+int InclusionOutputInfo::write(const std::string &output, const std::string &format) {
+  int result = 0;
+  if (format.compare("LEGACY") == 0) {
+    result = write_legacy(output);
+  } else if (format.compare("HDF5") == 0) {
+    result = write_hdf5(output);
+  } else {
+    string message = "Unknown format mode: \"" + format + "\"";
+    throw UnrecognizedConfigurationException(message);
+  }
+  return result;
+}
+
+int InclusionOutputInfo::write_hdf5(const std::string &file_name) {
+  List<string> *rec_name_list = new List<string>(1);
+  List<string> *rec_type_list = new List<string>(1);
+  List<void *> *rec_ptr_list = new List<void *>(1);
+  string str_type, str_name;
+  rec_name_list->set(0, "NSPH");
+  rec_type_list->set(0, "INT32_(1)");
+  rec_ptr_list->set(0, &nsph);
+  rec_name_list->append("LI");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&li);
+  rec_name_list->append("LE");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&le);
+  rec_name_list->append("LM");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&lm);
+  rec_name_list->append("MXNDM");
+  rec_type_list->append("INT64_(1)");
+  rec_ptr_list->append(&mxndm);
+  rec_name_list->append("INPOL");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&inpol);
+  rec_name_list->append("NPNT");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&npnt);
+  rec_name_list->append("NPNTTS");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&npntts);
+  rec_name_list->append("IAVM");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&iavm);
+  rec_name_list->append("ISAM");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&isam);
+  rec_name_list->append("JWTM");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&jwtm);
+  rec_name_list->append("VEC_SPH_X");
+  rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+  rec_ptr_list->append(vec_x_coords);
+  rec_name_list->append("VEC_SPH_Y");
+  rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+  rec_ptr_list->append(vec_y_coords);
+  rec_name_list->append("VEC_SPH_Z");
+  rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+  rec_ptr_list->append(vec_z_coords);
+  rec_name_list->append("TH_START");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&th);
+  rec_name_list->append("TH_STEP");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&thstp);
+  rec_name_list->append("TH_END");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&thlst);
+  rec_name_list->append("THS_START");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&ths);
+  rec_name_list->append("THS_STEP");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&thsstp);
+  rec_name_list->append("THS_END");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&thslst);
+  rec_name_list->append("PH_START");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&ph);
+  rec_name_list->append("PH_STEP");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&phstp);
+  rec_name_list->append("PH_END");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&phlst);
+  rec_name_list->append("PHS_START");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&phs);
+  rec_name_list->append("PHS_STEP");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&phsstp);
+  rec_name_list->append("PHS_END");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&phslst);
+  rec_name_list->append("EXRI");
+  rec_type_list->append("FLOAT64_(1)");
+  rec_ptr_list->append(&exri);
+  rec_name_list->append("IDFC");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&idfc);
+  rec_name_list->append("XI1");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&_first_xi);
+  rec_name_list->append("NXI");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&xi_block_size);
+  rec_name_list->append("VEC_JXI");
+  rec_type_list->append("INT32_(" + to_string(xi_block_size) + ")");
+  rec_ptr_list->append(vec_jxi);
+  rec_name_list->append("VEC_IER");
+  rec_type_list->append("INT16_(" + to_string(xi_block_size) + ")");
+  rec_ptr_list->append(vec_ier);
+  rec_name_list->append("VEC_VK");
+  rec_type_list->append("FLOAT64_(" + to_string(xi_block_size) + ")");
+  rec_ptr_list->append(vec_vk);
+  rec_name_list->append("VEC_XI");
+  rec_type_list->append("FLOAT64_(" + to_string(xi_block_size) + ")");
+  rec_ptr_list->append(vec_xi);
+  rec_name_list->append("NCONF");
+  rec_type_list->append("INT32_(1)");
+  rec_ptr_list->append(&configurations);
+  rec_name_list->append("VEC_SPH_SIZES");
+  rec_type_list->append("FLOAT64_(" + to_string(xi_block_size * (configurations + 1)) + ")");
+  rec_ptr_list->append(vec_sphere_sizes);
+  rec_name_list->append("VEC_SPH_REFRI");
+  rec_type_list->append("FLOAT64_(" + to_string(2 * xi_block_size * (configurations+ 1)) + ")");
+  rec_ptr_list->append(vec_sphere_ref_indices);
+  str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+  rec_name_list->append("VEC_SCS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_scs1);
+  rec_name_list->append("VEC_SCS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_scs2);
+  rec_name_list->append("VEC_ABS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_abs1);
+  rec_name_list->append("VEC_ABS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_abs2);
+  rec_name_list->append("VEC_EXS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_exs1);
+  rec_name_list->append("VEC_EXS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_exs2);
+  rec_name_list->append("VEC_ALBEDS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_albeds1);
+  rec_name_list->append("VEC_ALBEDS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_albeds2);
+  rec_name_list->append("VEC_SCSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_scsrt1);
+  rec_name_list->append("VEC_SCSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_scsrt2);
+  rec_name_list->append("VEC_ABSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_absrt1);
+  rec_name_list->append("VEC_ABSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_absrt2);
+  rec_name_list->append("VEC_EXSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_exsrt1);
+  rec_name_list->append("VEC_EXSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_exsrt2);
+  rec_name_list->append("VEC_QSCHU1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_qschu1);
+  rec_name_list->append("VEC_QSCHU2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_qschu2);
+  rec_name_list->append("VEC_PSCHU1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_pschu1);
+  rec_name_list->append("VEC_PSCHU2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_pschu2);
+  rec_name_list->append("VEC_S0MAG1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_s0mag1);
+  rec_name_list->append("VEC_S0MAG2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_s0mag2);
+  rec_name_list->append("VEC_COSAV1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_cosav1);
+  rec_name_list->append("VEC_COSAV2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_cosav2);
+  rec_name_list->append("VEC_RAPRS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_raprs1);
+  rec_name_list->append("VEC_RAPRS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_raprs2);
+  rec_name_list->append("VEC_FK1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fk1);
+  rec_name_list->append("VEC_FK2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fk2);
+  str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+  rec_name_list->append("VEC_FSAS11");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fsas11);
+  rec_name_list->append("VEC_FSAS21");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fsas21);
+  rec_name_list->append("VEC_FSAS12");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fsas12);
+  rec_name_list->append("VEC_FSAS22");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_fsas22);
+  str_type = "FLOAT64_(" + to_string(_num_theta) + ")";
+  rec_name_list->append("VEC_DIR_TIDG");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tidg);
+  str_type = "FLOAT64_(" + to_string(_num_phi) + ")";
+  rec_name_list->append("VEC_DIR_PIDG");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_pidg);
+  str_type = "FLOAT64_(" + to_string(_num_thetas) + ")";
+  rec_name_list->append("VEC_DIR_TSDG");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tsdg);
+  str_type = "FLOAT64_(" + to_string(_num_phis) + ")";
+  rec_name_list->append("VEC_DIR_PSDG");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_psdg);
+  str_type = "FLOAT64_(" + to_string(ndirs) + ")";
+  rec_name_list->append("VEC_DIR_SCAND");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_scand);
+  rec_name_list->append("VEC_DIR_CFMP");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_cfmp);
+  rec_name_list->append("VEC_DIR_CFSP");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_cfsp);
+  rec_name_list->append("VEC_DIR_SFMP");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sfmp);
+  rec_name_list->append("VEC_DIR_SFSP");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sfsp);
+  str_type = "FLOAT64_(" + to_string(3 * ndirs) + ")";
+  rec_name_list->append("VEC_DIR_UN");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_un);
+  rec_name_list->append("VEC_DIR_UNS");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_uns);
+  str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+  rec_name_list->append("VEC_DIR_SCS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_scs1);
+  rec_name_list->append("VEC_DIR_SCS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_scs2);
+  rec_name_list->append("VEC_DIR_ABS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_abs1);
+  rec_name_list->append("VEC_DIR_ABS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_abs2);
+  rec_name_list->append("VEC_DIR_EXS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_exs1);
+  rec_name_list->append("VEC_DIR_EXS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_exs2);
+  rec_name_list->append("VEC_DIR_ALBEDS1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_albeds1);
+  rec_name_list->append("VEC_DIR_ALBEDS2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_albeds2);
+  rec_name_list->append("VEC_DIR_SCSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_scsrt1);
+  rec_name_list->append("VEC_DIR_SCSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_scsrt2);
+  rec_name_list->append("VEC_DIR_ABSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_absrt1);
+  rec_name_list->append("VEC_DIR_ABSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_absrt2);
+  rec_name_list->append("VEC_DIR_EXSRT1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_exsrt1);
+  rec_name_list->append("VEC_DIR_EXSRT2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_exsrt2);
+  str_type = "FLOAT64_(" + to_string(2 * ndirs * xi_block_size) + ")"; 
+  rec_name_list->append("VEC_DIR_FSAS11");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fsas11);
+  rec_name_list->append("VEC_DIR_FSAS21");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fsas21);
+  rec_name_list->append("VEC_DIR_FSAS12");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fsas12);
+  rec_name_list->append("VEC_DIR_FSAS22");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fsas22);
+  rec_name_list->append("VEC_DIR_SAS11");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sas11);
+  rec_name_list->append("VEC_DIR_SAS21");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sas21);
+  rec_name_list->append("VEC_DIR_SAS12");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sas12);
+  rec_name_list->append("VEC_DIR_SAS22");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_sas22);
+  str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+  rec_name_list->append("VEC_DIR_QSCHU1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_qschu1);
+  rec_name_list->append("VEC_DIR_QSCHU2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_qschu2);
+  rec_name_list->append("VEC_DIR_PSCHU1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_pschu1);
+  rec_name_list->append("VEC_DIR_PSCHU2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_pschu2);
+  rec_name_list->append("VEC_DIR_S0MAG1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_s0mag1);
+  rec_name_list->append("VEC_DIR_S0MAG2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_s0mag2);
+  rec_name_list->append("VEC_DIR_COSAV1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_cosav1);
+  rec_name_list->append("VEC_DIR_COSAV2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_cosav2);
+  rec_name_list->append("VEC_DIR_RAPR1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_rapr1);
+  rec_name_list->append("VEC_DIR_RAPR2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_rapr2);
+  rec_name_list->append("VEC_DIR_FL1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fl1);
+  rec_name_list->append("VEC_DIR_FL2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fl2);
+  rec_name_list->append("VEC_DIR_FR1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fr1);
+  rec_name_list->append("VEC_DIR_FR2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fr2);
+  rec_name_list->append("VEC_DIR_FK1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fk1);
+  rec_name_list->append("VEC_DIR_FK2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fk2);
+  rec_name_list->append("VEC_DIR_FX1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fx1);
+  rec_name_list->append("VEC_DIR_FX2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fx2);
+  rec_name_list->append("VEC_DIR_FY1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fy1);
+  rec_name_list->append("VEC_DIR_FY2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fy2);
+  rec_name_list->append("VEC_DIR_FZ1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fz1);
+  rec_name_list->append("VEC_DIR_FZ2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_fz2);
+  rec_name_list->append("VEC_DIR_TQEL1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqel1);
+  rec_name_list->append("VEC_DIR_TQEL2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqel2);
+  rec_name_list->append("VEC_DIR_TQER1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqer1);
+  rec_name_list->append("VEC_DIR_TQER2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqer2);
+  rec_name_list->append("VEC_DIR_TQEK1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqek1);
+  rec_name_list->append("VEC_DIR_TQEK2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqek2);
+  rec_name_list->append("VEC_DIR_TQEX1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqex1);
+  rec_name_list->append("VEC_DIR_TQEX2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqex2);
+  rec_name_list->append("VEC_DIR_TQEY1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqey1);
+  rec_name_list->append("VEC_DIR_TQEY2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqey2);
+  rec_name_list->append("VEC_DIR_TQEZ1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqez1);
+  rec_name_list->append("VEC_DIR_TQEZ2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqez2);
+  rec_name_list->append("VEC_DIR_TQSL1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsl1);
+  rec_name_list->append("VEC_DIR_TQSL2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsl2);
+  rec_name_list->append("VEC_DIR_TQSR1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsr1);
+  rec_name_list->append("VEC_DIR_TQSR2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsr2);
+  rec_name_list->append("VEC_DIR_TQSK1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsk1);
+  rec_name_list->append("VEC_DIR_TQSK2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsk2);
+  rec_name_list->append("VEC_DIR_TQSX1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsx1);
+  rec_name_list->append("VEC_DIR_TQSX2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsx2);
+  rec_name_list->append("VEC_DIR_TQSY1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsy1);
+  rec_name_list->append("VEC_DIR_TQSY2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsy2);
+  rec_name_list->append("VEC_DIR_TQSZ1");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsz1);
+  rec_name_list->append("VEC_DIR_TQSZ2");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_tqsz2);
+  str_type = "FLOAT64_(" + to_string(16 * ndirs * xi_block_size) + ")";
+  rec_name_list->append("VEC_DIR_MULL");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_mull);
+  rec_name_list->append("VEC_DIR_MULLLR");
+  rec_type_list->append(str_type);
+  rec_ptr_list->append(vec_dir_mulllr);
+  
+  // Convert the lists to arrays and write them to HDF5
+  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 = new FileSchema(rec_num, rec_types, rec_names);
+  HDFFile *hdf_file = HDFFile::from_schema(*schema, file_name, H5F_ACC_TRUNC);
+  for (int ri = 0; ri < rec_num; ri++)
+    hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+  hdf_file->close();
+  
+  // Clean memory
+  delete rec_name_list;
+  delete rec_type_list;
+  delete rec_ptr_list;
+  delete[] rec_names;
+  delete[] rec_types;
+  delete[] rec_pointers;
+  delete schema;
+  delete hdf_file;
+  return 0;
+}
+
+int InclusionOutputInfo::write_legacy(const std::string &output) {
+  const dcomplex cc0 = 0.0 + I * 0.0;
+  int result = 0;
+  FILE *p_outfile = fopen(output.c_str(), "w");
+  if (p_outfile != NULL) {
+    if (vec_xi[0] == 1) {
+      fprintf(p_outfile, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
+#ifdef USE_ILP64
+      fprintf(
+	p_outfile, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
+	nsph, li, le, mxndm, inpol, npnt, npntts,
+	iavm, iavm
+      );
+#else
+      fprintf(
+	p_outfile, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
+	nsph, li, le, mxndm, inpol, npnt, npntts,
+	iavm, iavm
+      );
+#endif // USE_ILP64
+      fprintf(p_outfile, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
+      for (int ri = 0; ri < nsph; ri++) {
+	fprintf(
+	  p_outfile, "%17.8lE%17.8lE%17.8lE\n",
+	  vec_x_coords[ri], vec_y_coords[ri], vec_z_coords[ri]
+	);
+      }
+      fprintf(p_outfile, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+      fprintf(
+        p_outfile, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+	th, thstp, thlst, ths, thsstp, thslst
+      );
+      fprintf(p_outfile, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+      fprintf(
+        p_outfile, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+	ph, phstp, phlst, phs, phsstp, phslst
+      );
+      fprintf(p_outfile, " READ(IR,*)JWTM\n");
+      fprintf(p_outfile, " %5d\n", jwtm);
+      fprintf(p_outfile, "  READ(ITIN)NSPHT\n");
+      fprintf(p_outfile, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");
+      fprintf(p_outfile, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
+      fprintf(p_outfile, "  READ(ITIN)(XIV(I),I=1,NXI)\n");
+      fprintf(p_outfile, "  READ(ITIN)NSHL(I),ROS(I)\n");
+      fprintf(p_outfile, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n \n");
+      fprintf(p_outfile, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+      if (idfc < 0) {
+	fprintf( p_outfile, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n \n", vec_vk[0]);
+      }
+      for (int jxi = 0; jxi < xi_block_size; jxi++) {
+	int done_dirs = 0;
+	double alamb = 2.0 * 3.141592653589793238 / vec_vk[jxi];
+	fprintf(p_outfile, "========== JXI =%3d ====================\n", vec_jxi[jxi]);
+	if (idfc >= 0) {
+	  fprintf(p_outfile, "  VK=%15.7lE, XI=%15.7lE\n", vec_vk[jxi], vec_xi[jxi]);
+	} else {
+	  fprintf(p_outfile, "  XI=%15.7lE\n", vec_xi[jxi]);
+	}
+	if (vec_ier[jxi] == 1) {
+	  fprintf(p_outfile, "  STOP IN INDME\n");
+	  break;
+	} else if (vec_ier[jxi] == 2) {
+	  fprintf(p_outfile, "  STOP IN OSPV\n");
+	  break;
+	}
+	for (int i168 = 1; i168 <= configurations; i168++) {
+	  if (vec_sphere_ref_indices[i168 - 1] == cc0) {
+	    fprintf(p_outfile, "  SPHERE N.%2d: SIZE=%15.7lE\n", i168, vec_sphere_sizes[i168 - 1]);
+	  } else {
+	    fprintf(
+	      p_outfile, "  SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
+	      i168, vec_sphere_sizes[i168 - 1], real(vec_sphere_ref_indices[i168 - 1]),
+	      imag(vec_sphere_ref_indices[i168 - 1])
+	    );
+	  }
+	} // i168 configuration loop
+	fprintf(
+	  p_outfile, "  EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
+	  vec_sphere_sizes[configurations], real(vec_sphere_ref_indices[configurations]),
+	  imag(vec_sphere_ref_indices[configurations])
+	);
+	fprintf(p_outfile, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);
+	if (inpol == 0) fprintf(p_outfile, "   LIN -1\n");
+	else fprintf(p_outfile, "  CIRC -1\n");
+	fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+	fprintf(
+	  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+	  vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi], vec_albeds1[jxi]
+	);
+	fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
+		-1, alamb, vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi]
+	);
+	fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
+	fprintf(
+	  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+	  vec_scsrt1[jxi], vec_absrt1[jxi], vec_exsrt1[jxi]
+	);
+	fprintf(
+	  p_outfile, "  FSAS(1,1)=%15.7lE%15.7lE   FSAS(2,1)=%15.7lE%15.7lE\n",
+	  real(vec_fsas11[jxi]), imag(vec_fsas11[jxi]),
+	  real(vec_fsas21[jxi]), imag(vec_fsas21[jxi])
+	);
+	fprintf(
+	  p_outfile, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+	  vec_qschu1[jxi], vec_pschu1[jxi], vec_s0mag1[jxi]
+	);
+	fprintf(
+	  p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
+	  vec_cosav1[jxi], vec_raprs1[jxi]
+	);
+	fprintf(p_outfile, "  Fk=%15.7lE\n", vec_fk1[jxi]);
+	if (inpol == 0) fprintf(p_outfile, "   LIN  1\n");
+	else fprintf(p_outfile, "  CIRC  1\n");
+	fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+	fprintf(
+	  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+	  vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi], vec_albeds2[jxi]
+	);
+	fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
+		1, alamb, vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi]
+	);
+	fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
+	fprintf(
+	  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+	  vec_scsrt2[jxi], vec_absrt2[jxi], vec_exsrt2[jxi]
+	);
+	fprintf(
+	  p_outfile, "  FSAS(2,2)=%15.7lE%15.7lE   FSAS(1,2)=%15.7lE%15.7lE\n",
+	  real(vec_fsas22[jxi]), imag(vec_fsas22[jxi]),
+	  real(vec_fsas12[jxi]), imag(vec_fsas12[jxi])
+	);
+	fprintf(
+	  p_outfile, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+	  vec_qschu2[jxi], vec_pschu2[jxi], vec_s0mag2[jxi]
+	);
+	fprintf(
+	  p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
+	  vec_cosav2[jxi], vec_raprs2[jxi]
+	);
+	fprintf(p_outfile, "  Fk=%15.7lE\n", vec_fk2[jxi]);
+	fprintf(
+	  p_outfile, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n",
+	  (real(vec_fsas11[jxi]) - real(vec_fsas22[jxi])) / real(vec_fsas11[jxi])
+	);
+	fprintf(
+	  p_outfile, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n",
+	  (imag(vec_fsas11[jxi]) - imag(vec_fsas22[jxi])) / imag(vec_fsas11[jxi])
+	);
+	for (int jth = 0; jth < _num_theta; jth++) {
+	  for (int jph = 0; jph < _num_phi; jph++) {
+	    for (int jths = 0; jths < _num_thetas; jths++) {
+	      for (int jphs = 0; jphs < _num_phis; jphs++) {
+		bool goto290 = isam >= 0 && (jths > 0 || jphs > 0);
+		int dir_index = ndirs * jxi + done_dirs;
+		fprintf(
+			p_outfile, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
+			jth + 1, jph + 1, jths + 1, jphs +1
+		);
+		fprintf(
+		  p_outfile, "  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",
+		  th + jth * thstp, ph + jph * phstp, ths + jths * thsstp, phs + jphs * phsstp
+		);
+		fprintf(p_outfile, "  SCAND=%10.3lE\n", vec_dir_scand[done_dirs]);
+		fprintf(
+		  p_outfile, "  CFMP=%15.7lE, SFMP=%15.7lE\n",
+		  vec_dir_cfmp[done_dirs], vec_dir_sfmp[done_dirs]
+		);
+		fprintf(
+		  p_outfile, "  CFSP=%15.7lE, SFSP=%15.7lE\n",
+		  vec_dir_cfsp[done_dirs], vec_dir_sfsp[done_dirs]
+		);
+		if (isam >= 0) {
+		  fprintf(
+		    p_outfile, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n",
+		    vec_dir_un[done_dirs],
+		    vec_dir_un[done_dirs + 1],
+		    vec_dir_un[done_dirs + 2]
+		  );
+		  fprintf(
+		    p_outfile, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n",
+		    vec_dir_uns[done_dirs],
+		    vec_dir_uns[done_dirs + 1],
+		    vec_dir_uns[done_dirs + 2]
+		  );
+		} else { // label 214
+		  fprintf(
+		    p_outfile, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n",
+		    vec_dir_un[done_dirs],
+		    vec_dir_un[done_dirs + 1],
+		    vec_dir_un[done_dirs + 2]
+		  );
+		}
+		fprintf(p_outfile, "     SINGLE SCATTERER\n");
+		if (inpol == 0) {
+		  fprintf(p_outfile, "   LIN -1\n");
+		} else {
+		  fprintf(p_outfile, "  CIRC -1\n");
+		}
+		// label 275
+		fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+		fprintf(
+		  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+		  vec_dir_scs1[dir_index],
+		  vec_dir_abs1[dir_index],
+		  vec_dir_exs1[dir_index],
+		  vec_dir_albeds1[dir_index]
+		);
+		fprintf(
+		  p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
+		  -1, alamb, vec_dir_scs1[dir_index], vec_dir_abs1[dir_index],
+		  vec_dir_exs1[dir_index]
+		);
+		fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+		fprintf(
+		  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+		  vec_dir_scsrt1[dir_index],
+		  vec_dir_absrt1[dir_index],
+		  vec_dir_exsrt1[dir_index]
+		);
+		fprintf(
+		  p_outfile, "  FSAS(1,1)=%15.7lE%15.7lE   FSAS(2,1)=%15.7lE%15.7lE\n",
+		  real(vec_dir_fsas11[dir_index]), imag(vec_dir_fsas11[dir_index]),
+		  real(vec_dir_fsas21[dir_index]), imag(vec_dir_fsas21[dir_index])
+		);
+		fprintf(
+		  p_outfile, "   SAS(1,1)=%15.7lE%15.7lE    SAS(2,1)=%15.7lE%15.7lE\n",
+		  real(vec_dir_sas11[dir_index]), imag(vec_dir_sas11[dir_index]),
+		  real(vec_dir_sas21[dir_index]), imag(vec_dir_sas21[dir_index])
+		);
+		fprintf(
+		  p_outfile, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+		  vec_dir_qschu1[dir_index],
+		  vec_dir_pschu1[dir_index],
+		  vec_dir_s0mag1[dir_index]
+		);
+		if (!goto290) {
+		  fprintf(
+		    p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
+		    vec_dir_cosav1[dir_index], vec_dir_rapr1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n",
+		    vec_dir_fl1[dir_index], vec_dir_fr1[dir_index],
+		    vec_dir_fk1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+		    vec_dir_fx1[dir_index], vec_dir_fy1[dir_index],
+		    vec_dir_fz1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQEl=%15.7lE,  TQEr=%15.7lE,  TQEk=%15.7lE\n",
+		    vec_dir_tqel1[dir_index], vec_dir_tqer1[dir_index],
+		    vec_dir_tqek1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n",
+		    vec_dir_tqsl1[dir_index], vec_dir_tqsr1[dir_index],
+		    vec_dir_tqsk1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n",
+		    vec_dir_tqex1[dir_index], vec_dir_tqey1[dir_index],
+		    vec_dir_tqez1[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n",
+		    vec_dir_tqsx1[dir_index], vec_dir_tqsy1[dir_index],
+		    vec_dir_tqsz1[dir_index]
+		  );
+		} // goto290 switch
+		if (inpol == 0) {
+		  fprintf(p_outfile, "   LIN  1\n");
+		} else {
+		  fprintf(p_outfile, "  CIRC  1\n");
+		}
+		// label 275
+		fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+		fprintf(
+		  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+		  vec_dir_scs2[dir_index],
+		  vec_dir_abs2[dir_index],
+		  vec_dir_exs2[dir_index],
+		  vec_dir_albeds2[dir_index]
+		);
+		fprintf(
+		  p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
+		  1, alamb, vec_dir_scs2[dir_index], vec_dir_abs2[dir_index],
+		  vec_dir_exs2[dir_index]
+		);
+		fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+		fprintf(
+		  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+		  vec_dir_scsrt2[dir_index],
+		  vec_dir_absrt2[dir_index],
+		  vec_dir_exsrt2[dir_index]
+		);
+		fprintf(
+		  p_outfile, "  FSAS(2,2)=%15.7lE%15.7lE   FSAS(1,2)=%15.7lE%15.7lE\n",
+		  real(vec_dir_fsas22[dir_index]), imag(vec_dir_fsas22[dir_index]),
+		  real(vec_dir_fsas12[dir_index]), imag(vec_dir_fsas12[dir_index])
+		);
+		fprintf(
+		  p_outfile, "   SAS(2,2)=%15.7lE%15.7lE    SAS(1,2)=%15.7lE%15.7lE\n",
+		  real(vec_dir_sas22[dir_index]), imag(vec_dir_sas22[dir_index]),
+		  real(vec_dir_sas12[dir_index]), imag(vec_dir_sas12[dir_index])
+		);
+		fprintf(
+		  p_outfile, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+		  vec_dir_qschu2[dir_index],
+		  vec_dir_pschu2[dir_index],
+		  vec_dir_s0mag2[dir_index]
+		);
+		if (!goto290) {
+		  fprintf(
+		    p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
+		    vec_dir_cosav2[dir_index], vec_dir_rapr2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n",
+		    vec_dir_fl2[dir_index], vec_dir_fr2[dir_index],
+		    vec_dir_fk2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+		    vec_dir_fx2[dir_index], vec_dir_fy2[dir_index],
+		    vec_dir_fz2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQEl=%15.7lE,  TQEr=%15.7lE,  TQEk=%15.7lE\n",
+		    vec_dir_tqel2[dir_index], vec_dir_tqer2[dir_index],
+		    vec_dir_tqek2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n",
+		    vec_dir_tqsl2[dir_index], vec_dir_tqsr2[dir_index],
+		    vec_dir_tqsk2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n",
+		    vec_dir_tqex2[dir_index], vec_dir_tqey2[dir_index],
+		    vec_dir_tqez2[dir_index]
+		  );
+		  fprintf(
+		    p_outfile, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n",
+		    vec_dir_tqsx2[dir_index], vec_dir_tqsy2[dir_index],
+		    vec_dir_tqsz2[dir_index]
+		  );
+		} // goto290 switch
+		fprintf(
+		  p_outfile, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n",
+		  (real(vec_dir_fsas11[dir_index]) - real(vec_dir_fsas22[dir_index])) / real(vec_dir_fsas11[dir_index])
+		);
+		fprintf(
+		  p_outfile, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n",
+		  (imag(vec_dir_fsas11[dir_index]) - imag(vec_dir_fsas22[dir_index])) / imag(vec_dir_fsas11[dir_index])
+		);
+		fprintf(p_outfile, "  MULL\n");
+		for (int i = 0; i < 4; i++) {
+		  int mull_index = 16 * dir_index + 4 * i;
+		  fprintf(
+		    p_outfile, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		    vec_dir_mull[mull_index],
+		    vec_dir_mull[mull_index + 1],
+		    vec_dir_mull[mull_index + 2],
+		    vec_dir_mull[mull_index + 3]
+		  );
+		} // i MULL loop
+		fprintf(p_outfile, "  MULLLR\n");
+		for (int i = 0; i < 4; i++) {
+		  int mull_index = 16 * dir_index + 4 * i;
+		  fprintf(
+		    p_outfile, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		    vec_dir_mulllr[mull_index],
+		    vec_dir_mulllr[mull_index + 1],
+		    vec_dir_mulllr[mull_index + 2],
+		    vec_dir_mulllr[mull_index + 3]
+		  );
+		} // i MULLR loop
+		if (iavm != 0) {
+		  fprintf(p_outfile, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);
+		  if (inpol == 0) fprintf(p_outfile, "   LIN\n");
+		  else fprintf(p_outfile, "  CIRC\n");
+		  // label 318
+		  fprintf(p_outfile, "  MULL\n");
+		  for (int i = 0; i < 4; i++) {
+		    int mul_dir_index = 16 * dir_index + 4 * i;
+		    fprintf(
+		      p_outfile, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		      vec_dir_mull[mul_dir_index],
+		      vec_dir_mull[mul_dir_index + 1],
+		      vec_dir_mull[mul_dir_index + 2],
+		      vec_dir_mull[mul_dir_index + 3]
+		    );
+		  } // i MULL loop
+		  fprintf(p_outfile, "  MULLLR\n");
+		  for (int i = 0; i < 4; i++) {
+		    int mul_dir_index = 16 * dir_index + 4 * i;
+		    fprintf(
+		      p_outfile, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		      vec_dir_mulllr[mul_dir_index],
+		      vec_dir_mulllr[mul_dir_index + 1],
+		      vec_dir_mulllr[mul_dir_index + 2],
+		      vec_dir_mulllr[mul_dir_index + 3]
+		    );
+		  } // i MULLR loop
+		} // iavm != 0 block
+		done_dirs++;
+	      } // jphs loop
+	    } // jths loop
+	  } // jph loop
+	} // jth loop
+      } // jxi wavelength loop
+    }
+    fclose(p_outfile);
+  } else {
+    result = -1;
+  }
+  return result;
+}
+
+#ifdef MPI_VERSION
+int InclusionOutputInfo::mpireceive(const mixMPI* mpidata, int pid) {
+  int result = 0;
+  int chk_inpol, chk_iavm, chk_isam, chk_num_theta, chk_num_thetas;
+  int chk_num_phi, chk_num_phis, chk_ndirs, chk_idfc, chk_configs;
+  double chk_exri;
+  MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_iavm, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  result += (chk_inpol == inpol) ? 0 : 1;
+  result += (chk_iavm == iavm) ? 0 : 1;
+  result += (chk_isam == isam) ? 0 : 1;
+  result += (chk_num_theta == _num_theta) ? 0 : 1;
+  result += (chk_num_thetas == _num_thetas) ? 0 : 1;
+  result += (chk_num_phi == _num_phi) ? 0 : 1;
+  result += (chk_num_phis == _num_phis) ? 0 : 1;
+  result += (chk_ndirs == ndirs) ? 0 : 1;
+  result += (chk_exri == exri) ? 0 : 1;
+  result += (chk_idfc == idfc) ? 0 : 1;
+  result += (chk_configs == configurations) ? 0 : 1;
+  if (result == 0) {
+    int xi1, offset, chunk_size;
+    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+    MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    // Receive vectors of single values per scale
+    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    offset = xi1 - _first_xi;
+    MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_raprs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_raprs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+    // Receive vectors whose sizes depend on configurations
+    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    offset = (xi1 - _first_xi) * (configurations + 1);
+    MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+    // Receive vectors whose sizes depend on directions and wavelengths
+    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    offset = (xi1 - _first_xi) * ndirs;
+    MPI_Recv(vec_dir_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_rapr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_rapr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+    MPI_Recv(vec_dir_fl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_fz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqel1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqel2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqer1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqer2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqex1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqex2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqey1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqey2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqez1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqez2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_tqsz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_mull + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(vec_dir_mulllr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  }  else {
+    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+  }
+  return result;
+}
+
+int InclusionOutputInfo::mpisend(const mixMPI *mpidata) {
+  int result = 0;
+  int chunk_size;
+  // Send output metadata for configuration cross-check
+  MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  // Wait for process 0 to cross-check the configuration
+  MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  if (result == 0) {
+    // Process 0 confirmed the consistency of configuration. Send the data.
+    // Send vectors of single values per scale
+    MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_scs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_scs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_abs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_abs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_exs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_exs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_albeds1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_albeds2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_scsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_scsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_absrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_absrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_exsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_exsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_qschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_qschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_pschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_pschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_s0mag1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_s0mag2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_cosav1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_cosav2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_raprs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_raprs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fk1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fk2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fsas11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fsas21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fsas22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_fsas12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+
+    // Send vectors whose sizes depend on configurations
+    chunk_size = xi_block_size * (configurations + 1);
+    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+
+    // Send vectors whose sizes depend on directions and wavelengths
+    chunk_size = xi_block_size * ndirs;
+    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_scs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_scs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_abs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_abs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_exs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_exs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_albeds1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_albeds2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_scsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_scsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_absrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_absrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_exsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_exsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fsas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fsas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fsas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fsas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_qschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_qschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_pschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_pschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_s0mag1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_s0mag2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_cosav1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_cosav2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_rapr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_rapr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_fz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqel1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqel2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqer1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqer2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqex1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqex2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqey1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqey2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqez1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqez2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_tqsz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_mull, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(vec_dir_mulllr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+  }
+  return result;
+}
+#endif // MPI_VERSION
+// >>> END OF InclusionOutputInfo CLASS IMPLEMENTATION <<<
-- 
GitLab