Commit f32264ab authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement SphereOutputInfo class and use it in np_sphere

parent 098577ec

src/include/outputs.h

+107 −21
Original line number Diff line number Diff line
@@ -467,9 +467,9 @@ public: 


  /*! \brief `ClusterOutputInfo` constructor from HDF5 input.

   *

   * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.

   * \param hdf5_name: `const string &` Path to the HDF5 file to be read.

   */   

  ClusterOutputInfo(const std::string &hdf5_file_name);

  ClusterOutputInfo(const std::string &hdf5_name);



  /*! \brief `ClusterOutputInfo` instance destroyer.

   */
@@ -479,13 +479,13 @@ public: 
   *

   * \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 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 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

    int first_xi = 1, int xi_length = 0

  );

  

  /*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
@@ -647,7 +647,7 @@ public: 
  int jwtm;

  //! \brief Vector of scale (wavelength) indices.

  int *vec_jxi;

  //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).

  //! \brief Vector of error severities (0 - success, 1 - INDME, 2 - OSPV).

  short *vec_ier;

  //! \brief Vector of vacuum wave numbers.

  double *vec_vk;
@@ -897,9 +897,9 @@ public: 


  /*! \brief `InclusionOutputInfo` constructor from HDF5 input.

   *

   * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.

   * \param hdf5_name: `const string &` Path to the HDF5 file to be read.

   */   

  InclusionOutputInfo(const std::string &hdf5_file_name);

  InclusionOutputInfo(const std::string &hdf5_name);



  /*! \brief `InclusionOutputInfo` instance destroyer.

   */
@@ -909,13 +909,13 @@ public: 
   *

   * \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 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 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

    int first_xi = 1, int xi_length = 0

  );

  

  /*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
@@ -1004,18 +1004,16 @@ protected: 
public:

  //! \brief Read-only view on the ID of the first scale

  const int &first_xi = _first_xi;

  //! \brief Number of spheres.

  int nsph;

  //! \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.
@@ -1054,18 +1052,106 @@ public: 
  int xi_block_size;

  //! \brief Index of the wavelength for T-matrix output.

  int jwtm;

  //! \brief Number of sphere types.

  int configurations;

  //! \brief Highest expansion order achieved in calculations.

  int lcalc;

  //! \brief Harmonic functions argument.

  dcomplex arg;

  //! \brief Vector of scale (wavelength) indices.

  int *vec_jxi;

  //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).

  //! \brief Vector of error severities (0 - success, 1 - DME).

  short *vec_ier;

  //! \brief Vector of vacuum wave numbers.

  double *vec_vk;

  //! \brief Vector of computed scales.

  double *vec_xi;

  //! \brief Vector of sphere sizes (one for every scale).

  //! \brief Vector of sphere sizes (one for every configuration and scale).

  double *vec_sphere_sizes;

  //! \brief Vector of sphere refractive indices  (one for every scale).

  //! \brief Vector of sphere refractive indices (one for every configuration and scale).

  dcomplex *vec_sphere_ref_indices;

  //! \brief Vector of sphere scattering cross-sections.

  double *vec_scs;

  //! \brief Vector of sphere absorption cross-sections.

  double *vec_abs;

  //! \brief Vector of sphere extinction cross-sections.

  double *vec_exs;

  //! \brief Vector of sphere albedos.

  double *vec_albeds;

  //! \brief Vector of sphere scattering-to-geometric cross-sections.

  double *vec_scsrt;

  //! \brief Vector of sphere absorption-to-geometric cross-sections.

  double *vec_absrt;

  //! \brief Vector of sphere extinction-to-geometric cross-sections.

  double *vec_exsrt;

  //! \brief Vector of sphere forward scattering amplitudes.

  dcomplex *vec_fsas;

  //! \brief Vector of sphere QSCHU.

  double *vec_qschu;

  //! \brief Vector of sphere PSCHU.

  double *vec_pschu;

  //! \brief Vector of sphere S0MAG.

  double *vec_s0mag;

  //! \brief Vector of sphere average asymmetry parameter.

  double *vec_cosav;

  //! \brief Vector of sphere average radiation pressure force (N).

  double *vec_raprs;

  //! \brief Vector of sphere average extinction torque along incidence direction (parallel polarization).

  double *vec_tqek1;

  //! \brief Vector of sphere average extinction torque along incidence direction (perpendicular polarization).

  double *vec_tqek2;

  //! \brief Vector of sphere average scattering torque along incidence direction (parallel polarization).

  double *vec_tqsk1;

  //! \brief Vector of sphere average scattering torque along incidence direction (perpendicular polarization).

  double *vec_tqsk2;

  //| \brief Vector of total forward scattering amplitudes.

  dcomplex *vec_fsat;

  //! \brief Vector of total QSCHU.

  double *vec_qschut;

  //! \brief Vector of total PSCHU.

  double *vec_pschut;

  //! \brief Vector of total S0MAG.

  double *vec_s0magt;

  //! \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 sphere differential scattering amplitude with polarization parallel to parallel incidence field.

  dcomplex *vec_dir_sas11;

  //! \brief Vector of sphere differential scattering amplitude with polarization perpendicular to the parallel incidence field.

  dcomplex *vec_dir_sas21;

  //! \brief Vector of sphere differential scattering amplitude with polarization perpendicular to perpendicular incidence field.

  dcomplex *vec_dir_sas12;

  //! \brief Vector of sphere differential scattering amplitude with polarization parallel the perpendicular incidence field.

  dcomplex *vec_dir_sas22;

  //! \brief Vector of differential radiation pressure force components along the X axis.

  double *vec_dir_fx;

  //! \brief Vector of differential radiation pressure force components along the Y axis.

  double *vec_dir_fy;

  //! \brief Vector of differential radiation pressure force components along the Z axis.

  double *vec_dir_fz;

  //! \brief Vector of sphere Mueller transormation matrices referred to meridional plane.

  double *vec_dir_muls;

  //! \brief Vector of sphere Mueller transormation matrices referred to scattering plane.

  double *vec_dir_mulslr;

  

  /*! \brief `SphereOutputInfo` default instance constructor.

   *
@@ -1082,9 +1168,9 @@ public: 


  /*! \brief `SphereOutputInfo` constructor from HDF5 input.

   *

   * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.

   * \param hdf5_name: `const string &` Path to the HDF5 file to be read.

   */   

  SphereOutputInfo(const std::string &hdf5_file_name);

  SphereOutputInfo(const std::string &hdf5_name);



  /*! \brief `InclusionOutputInfo` instance destroyer.

   */
@@ -1094,13 +1180,13 @@ public: 
   *

   * \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 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 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

    int first_xi = 1, int xi_length = 0

  );

  

  /*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.

src/libnptm/outputs.cpp

+1161 −36

File changed.

Preview size limit exceeded, changes collapsed.

src/sphere/sphere.cpp

+85 −133
Original line number Diff line number Diff line
@@ -51,6 +51,10 @@ 
#include "../include/TransitionMatrix.h"

#endif



#ifndef INCLUDE_OUTPUTS_H_

#include "../include/outputs.h"

#endif



using namespace std;



/*! \brief C++ implementation of SPH
@@ -90,6 +94,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
  }

  int s_nsph = sconf->number_of_spheres;

  int nsph = gconf->number_of_spheres;

  int configurations = sconf->configurations;

  if (s_nsph == nsph) {

    int isq, ibf;

    double *argi, *args, *gaps;
@@ -127,6 +132,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
    }

    double frx = 0.0, fry = 0.0, frz = 0.0;

    double cfmp, cfsp, sfmp, sfsp;

    const dcomplex cc0 = 0.0 * I + 0.0;

    int jw;

    int l_max = gconf->l_max;

    ParticleDescriptor *c1 = new ParticleDescriptorSphere(gconf, sconf);
@@ -150,48 +156,9 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
    argi = new double[1];

    args = new double[1];

    gaps = new double[2];

    FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w");

    fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");

    fprintf(

	    output,

	    " %5d%5d%5d%5d%5d%5d\n",

	    nsph,

	    l_max,

	    in_pol,

	    npnt,

	    npntts,

	    meridional_type

	    );

    fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");

    fprintf(

	    output,

	    "  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",

	    in_theta_start,

	    in_theta_step,

	    in_theta_end,

	    sc_theta_start,

	    sc_theta_step,

	    sc_theta_end

	    );

    fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");

    fprintf(

	    output,

	    "  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",

	    in_phi_start,

	    in_phi_step,

	    in_phi_end,

	    sc_phi_start,

	    sc_phi_step,

	    sc_phi_end

	    );

    fprintf(output, " READ(IR,*)JWTM\n");

    fprintf(output, " %5d\n", jwtm);

    fprintf(output, "  READ(ITIN)NSPHT\n");

    fprintf(output, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");

    fprintf(output, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");

    fprintf(output, "  READ(ITIN)(XIV(I),I=1,NXI)\n");

    fprintf(output, "  READ(ITIN)NSHL(I),ROS(I)\n");

    fprintf(output, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n \n");

    // FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w");

    mixMPI *mpidata = new mixMPI(); // Fake MPI configuration

    SphereOutputInfo *p_output = new SphereOutputInfo(sconf, gconf, mpidata);

    double sml = 1.0e-3;

    int nth = 0, nph = 0;

    if (in_theta_step != 0.0)
@@ -253,7 +220,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
    thdps(l_max, zpv);

    double exdc = sconf->exdc;

    double exri = sqrt(exdc);

    fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);

    fstream tppoan;

    string tppoan_name = output_path + "/c_TPPOAN";

    tppoan.open(tppoan_name.c_str(), ios::binary|ios::out);
@@ -271,9 +237,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou 


      for (int nsi = 0; nsi < nsph; nsi++)

	tppoan.write(reinterpret_cast<char *>(&(c1->iog[nsi])), sizeof(int));

      if (in_pol == 0) fprintf(output, "   LIN\n");

      else fprintf(output, "  CIRC\n");

      fprintf(output, " \n");

      double wp = sconf->wp;

      double xip = sconf->xip;

      double wn = wp / 3.0e8;
@@ -283,22 +246,24 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
      int nxi = sconf->number_of_scales;

      if (idfc < 0) {

	vk = xip * wn;

	fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);

	fprintf(output, " \n");

	p_output->vec_vk[0] = vk;

      }

      int ndirs = nkks;

      for (int jxi488 = 0; jxi488 < nxi; jxi488++) {

	int oindex = 0;

	int jxi = jxi488 + 1;

	logger->log("INFO: running scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");

	fprintf(output, "========== JXI =%3d ====================\n", jxi);

	double xi = sconf->get_scale(jxi488);

	if (idfc >= 0) {

	  vk = xi * wn;

	  vkarg = vk;

	  fprintf(output, "  VK=%15.7lE, XI=%15.7lE\n", xi, vk);

	  p_output->vec_vk[jxi488] = vk;

	  p_output->vec_xi[jxi488] = xi;

	} else { // IDFC < 0

	  vkarg = xi * vk;

	  sqsfi = 1.0 / (xi * xi);

	  fprintf(output, "  XI=%15.7lE\n", xi);

	  p_output->vec_vk[jxi488] = vk;

	  p_output->vec_xi[jxi488] = xi;

	}

	tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));

	for (int i132 = 0; i132 < nsph; i132++) {
@@ -329,10 +294,10 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
	  int lcalc = 0;

	  dme(l_max, i, npnt, npntts, vkarg, exdc, exri, c1, jer, lcalc, arg);

	  if (jer != 0) {

	    fprintf(output, "  STOP IN DME\n");

	    fprintf(output, "  AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, real(arg), imag(arg));

	    p_output->vec_ier[jxi] = 1;

	    p_output->lcalc = lcalc;

	    p_output->arg = arg;

	    tppoan.close();

	    fclose(output);

	    delete sconf;

	    delete gconf;

	    delete c1;
@@ -398,55 +363,47 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
	double sqk = vk * vk * exdc;

	aps(zpv, l_max, nsph, c1, sqk, gaps);

	rabas(in_pol, l_max, nsph, c1, tqse, tqspe, tqss, tqsps);

	int last_configuration = 0;

	for (int i170 = 0; i170 < nsph; i170++) {

	  int i = i170 + 1;

	  if (c1->iog[i170] >= i) {

	    last_configuration++;

	    oindex = jxi488 * configurations + last_configuration - 1;

	    double albeds = c1->sscs[i170] / c1->sexs[i170];

	    c1->sqscs[i170] *= sqsfi;

	    c1->sqabs[i170] *= sqsfi;

	    c1->sqexs[i170] *= sqsfi;

	    fprintf(output, "     SPHERE %2d\n", i);

	    if (c1->nshl[i170] != 1) {

	      fprintf(output, "  SIZE=%15.7lE\n", c1->vsz[i170]);

	      p_output->vec_sphere_ref_indices[oindex] = cc0;

	      p_output->vec_sphere_sizes[oindex] = c1->vsz[i170];

	    } else {

	      fprintf(

		      output,

		      "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",

		      c1->vsz[i170],

		      real(c1->vkt[i170]),

		      imag(c1->vkt[i170])

		      );

	    }

	    fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");

	    fprintf(

		    output,

		    " %14.7lE%15.7lE%15.7lE%15.7lE\n",

		    c1->sscs[i170], c1->sabs[i170],

		    c1->sexs[i170], albeds

		    );

	    fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");

	    fprintf(

		    output,

		    " %14.7lE%15.7lE%15.7lE\n",

		    c1->sqscs[i170], c1->sqabs[i170],

		    c1->sqexs[i170]

		    );

	    fprintf(output, "  FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i170]), imag(c1->fsas[i170]));

	      p_output->vec_sphere_ref_indices[oindex] = c1->vkt[i170];

	      p_output->vec_sphere_sizes[oindex] = c1->vsz[i170];

	    }

	    p_output->vec_scs[oindex] = c1->sscs[i170];

	    p_output->vec_abs[oindex] = c1->sabs[i170];

	    p_output->vec_exs[oindex] = c1->sexs[i170];

	    p_output->vec_albeds[oindex] = albeds;

	    p_output->vec_scsrt[oindex] = c1->sqscs[i170];

	    p_output->vec_absrt[oindex] = c1->sqabs[i170];

	    p_output->vec_exsrt[oindex] = c1->sqexs[i170];

	    p_output->vec_fsas[oindex] = c1->fsas[i170];

	    double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];

	    s0 = c1->fsas[i170] * exri;

	    double qschu = csch * imag(s0);

	    double pschu = csch * real(s0);

	    double s0mag = cs0 * cabs(s0);

	    fprintf(

		    output,

		    "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",

		    qschu, pschu, s0mag

		    );

	    p_output->vec_qschu[oindex] = qschu;

	    p_output->vec_pschu[oindex] = pschu;

	    p_output->vec_s0mag[oindex] = s0mag;

	    double rapr = c1->sexs[i170] - gaps[i170];

	    double cosav = gaps[i170] / c1->sscs[i170];

	    fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);

	    fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170], tqss[0][i170]);

	    fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);

	    p_output->vec_cosav[oindex] = cosav;

	    p_output->vec_raprs[oindex] = rapr;

	    p_output->vec_tqek1[oindex] = tqse[0][i170];

	    p_output->vec_tqsk1[oindex] = tqss[0][i170];

	    p_output->vec_tqek2[oindex] = tqse[1][i170];

	    p_output->vec_tqsk2[oindex] = tqss[1][i170];

	    tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));

	    tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));

	    double val = real(tqspe[0][i170]);
@@ -470,19 +427,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
	  } // End if iog[i170] >= i

	} // i170 loop

	if (nsph != 1) {

	  fprintf(output, "  FSAT=(%15.7lE,%15.7lE)\n", real(tfsas), imag(tfsas));

	  p_output->vec_fsat[jxi488] = tfsas;

	  double csch = 2.0 * vk * sqsfi / gcs;

	  s0 = tfsas * exri;

	  double qschu = csch * imag(s0);

	  double pschu = csch * real(s0);

	  double s0mag = cs0 * cabs(s0);

	  fprintf(

		  output,

		  "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",

		  qschu, pschu, s0mag

		  );

	  p_output->vec_qschut[jxi488] = qschu;

	  p_output->vec_pschut[jxi488] = pschu;

	  p_output->vec_s0magt[jxi488] = s0mag;

	}

	th = th1;

	int done_dirs = 0;

	for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section

	  int jth = jth486 + 1;

	  ph = ph1;
@@ -556,61 +512,51 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
		  tppoan.write(reinterpret_cast<char *>(&(u[1])), sizeof(double));

		  tppoan.write(reinterpret_cast<char *>(&(u[2])), sizeof(double));

		}

		fprintf(

			output,

			"********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",

			jth, jph, jths, jphs

			);

		fprintf(

			output,

			"  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",

			th, ph, ths, phs

			);

		fprintf(output, "  SCAND=%10.3lE\n", scan);

		fprintf(output, "  CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);

		fprintf(output, "  CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);

		p_output->vec_dir_scand[done_dirs] = scan;

		p_output->vec_dir_cfmp[done_dirs] = cfmp;

		p_output->vec_dir_cfsp[done_dirs] = cfsp;

		p_output->vec_dir_sfmp[done_dirs] = sfmp;

		p_output->vec_dir_sfsp[done_dirs] = sfsp;

		if (meridional_type >= 0) {

		  fprintf(output, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);

		  fprintf(output, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);

		  p_output->vec_dir_un[3 * done_dirs] = un[0];

		  p_output->vec_dir_un[3 * done_dirs + 1] = un[1];

		  p_output->vec_dir_un[3 * done_dirs + 2] = un[2];

		  p_output->vec_dir_uns[3 * done_dirs] = uns[0];

		  p_output->vec_dir_uns[3 * done_dirs + 1] = uns[1];

		  p_output->vec_dir_uns[3 * done_dirs + 2] = uns[2];

		} else {

		  fprintf(output, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);

		  p_output->vec_dir_un[3 * done_dirs] = un[0];

		  p_output->vec_dir_un[3 * done_dirs + 1] = un[1];

		  p_output->vec_dir_un[3 * done_dirs + 2] = un[2];

		  p_output->vec_dir_uns[3 * done_dirs] = 0.0;

		  p_output->vec_dir_uns[3 * done_dirs + 1] = 0.0;

		  p_output->vec_dir_uns[3 * done_dirs + 2] = 0.0;

		}

		sscr2(nsph, l_max, vk, exri, c1);

		last_configuration = 0;

		for (int ns226 = 0; ns226 < nsph; ns226++) {

		  int ns = ns226 + 1;

		  fprintf(output, "     SPHERE %2d\n", ns);

		  fprintf(

			  output, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",

			  real(c1->sas[ns226][0][0]), imag(c1->sas[ns226][0][0]),

			  real(c1->sas[ns226][1][0]), imag(c1->sas[ns226][1][0])

			  );

		  fprintf(

			  output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",

			  real(c1->sas[ns226][0][1]), imag(c1->sas[ns226][0][1]),

			  real(c1->sas[ns226][1][1]), imag(c1->sas[ns226][1][1])

			  );

		  if (jths == 1 && jphs == 1)

		    fprintf(

			    output, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",

			    frx, fry, frz

			    );

		  oindex = jxi488 * nsph * ndirs + nsph * done_dirs + ns226;

		  p_output->vec_dir_sas11[oindex] = c1->sas[ns226][0][0];

		  p_output->vec_dir_sas21[oindex] = c1->sas[ns226][1][0];

		  p_output->vec_dir_sas12[oindex] = c1->sas[ns226][0][1];

		  p_output->vec_dir_sas22[oindex] = c1->sas[ns226][1][1];

		  p_output->vec_dir_fx[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;

		  p_output->vec_dir_fy[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;

		  p_output->vec_dir_fz[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;

		  for (int i225 = 0; i225 < 16; i225++) c1->vint[i225] = c1->vints[ns226][i225];

		  mmulc(c1->vint, cmullr, cmul);

		  fprintf(output, "  MULS\n        ");

		  for (int imul = 0; imul < 4; imul++) {

		    int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

		    for (int jmul = 0; jmul < 4; jmul++) {

		      fprintf(output, "%15.7lE", cmul[imul][jmul]);

		      p_output->vec_dir_muls[muls_index + jmul] = cmul[imul][jmul];

		    }

		    if (imul < 3) fprintf(output, "\n        ");

		    else fprintf(output, "\n");

		  }

		  fprintf(output, "  MULSLR\n        ");

		  for (int imul = 0; imul < 4; imul++) {

		    int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

		    for (int jmul = 0; jmul < 4; jmul++) {

		      fprintf(output, "%15.7lE", cmullr[imul][jmul]);

		      p_output->vec_dir_mulslr[muls_index + jmul] = cmullr[imul][jmul];

		    }

		    if (imul < 3) fprintf(output, "\n        ");

		    else fprintf(output, "\n");

		  }

		  for (int vi = 0; vi < 16; vi++) {

		    double value = real(c1->vint[vi]);
@@ -625,6 +571,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
		  }

		} // ns226 loop

		if (meridional_type < 1) phs += sc_phi_step;

		done_dirs++;

	      } // jphs480 loop

	      if (meridional_type <= 1) thsl += sc_theta_step;

	    } // jths482 loop
@@ -632,13 +579,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou 
	  } // jph484 loop on elevation

	  th += in_theta_step;

	} // jth486 loop on azimuth

	p_output->vec_jxi[jxi488] = jxi488 + 1;

	logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");

      } //jxi488 loop on scales

      tppoan.close();

    } else { // In case TPPOAN could not be opened. Should never happen.

      logger->err("ERROR: failed to open TPPOAN file.\n");

    }

    fclose(output);

    // fclose(output);

    p_output->write(output_path + "/c_OSPH.hd5", "HDF5");

    p_output->write(output_path + "/c_OSPH", "LEGACY");

    delete p_output;

    delete mpidata;

    delete c1;

    for (int zi = l_max - 1; zi > -1; zi--) {

      for (int zj = 0; zj < 3; zj++) {