Create instructions to build C++ version of inclusion

lib_LTLIBRARIES=libnptm/libnptm.la

libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

if BUILDFORTRAN

PROGS=cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS

PROGS=cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu inclusion/np_inclusion sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS

bin_PROGRAMS=$(PROGS)

EDFBCLUSOURCES=../src/cluster/edfb_clu.f

CLUSOURCES=../src/cluster/clu.f

trapping_frfme_SOURCES=$(FRFMESOURCES)

trapping_lffft_SOURCES=$(LFFFTSOURCES)

cluster_np_cluster_SOURCES=../src/cluster/np_cluster.cpp ../src/cluster/cluster.cpp

inclusion_np_inclusion_SOURCES=../src/inclusion/np_inclusion.cpp ../src/inclusion/inclusion.cpp

sphere_np_sphere_SOURCES=../src/sphere/np_sphere.cpp ../src/sphere/sphere.cpp

trapping_np_trapping_SOURCES=../src/trapping/np_trapping.cpp ../src/trapping/cfrfme.cpp ../src/trapping/clffft.cpp

testing_test_TEDF_SOURCES=../src/testing/test_TEDF.cpp

testing_test_TTMS_SOURCES=../src/testing/test_TTMS.cpp

else

PROGS=cluster/np_cluster sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS

PROGS=cluster/np_cluster inclusion/np_inclusion sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS

bin_PROGRAMS=$(PROGS)

cluster_np_cluster_SOURCES=../src/cluster/np_cluster.cpp ../src/cluster/cluster.cpp

inclusion_np_inclusion_SOURCES=../src/inclusion/np_inclusion.cpp ../src/inclusion/inclusion.cpp

sphere_np_sphere_SOURCES=../src/sphere/np_sphere.cpp ../src/sphere/sphere.cpp

trapping_np_trapping_SOURCES=../src/trapping/np_trapping.cpp ../src/trapping/cfrfme.cpp ../src/trapping/clffft.cpp

testing_test_TEDF_SOURCES=../src/testing/test_TEDF.cpp

CONFIG_CLEAN_FILES =

CONFIG_CLEAN_VPATH_FILES =

@BUILDFORTRAN_FALSE@am__EXEEXT_1 = cluster/np_cluster$(EXEEXT) \

@BUILDFORTRAN_FALSE@	inclusion/np_inclusion$(EXEEXT) \

@BUILDFORTRAN_FALSE@	sphere/np_sphere$(EXEEXT) \

@BUILDFORTRAN_FALSE@	trapping/np_trapping$(EXEEXT) \

@BUILDFORTRAN_FALSE@	testing/test_TEDF$(EXEEXT) \

@BUILDFORTRAN_TRUE@	cluster/np_cluster$(EXEEXT) \

@BUILDFORTRAN_TRUE@	inclusion/edfb_inclu$(EXEEXT) \

@BUILDFORTRAN_TRUE@	inclusion/inclu$(EXEEXT) \

@BUILDFORTRAN_TRUE@	inclusion/np_inclusion$(EXEEXT) \

@BUILDFORTRAN_TRUE@	sphere/edfb_sph$(EXEEXT) \

@BUILDFORTRAN_TRUE@	sphere/sph$(EXEEXT) \

@BUILDFORTRAN_TRUE@	sphere/np_sphere$(EXEEXT) \

inclusion_inclu_DEPENDENCIES = libnptm/libnptm.la \

	$(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \

	$(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1)

am__inclusion_np_inclusion_SOURCES_DIST =  \

	../src/inclusion/np_inclusion.cpp \

	../src/inclusion/inclusion.cpp

@BUILDFORTRAN_FALSE@am_inclusion_np_inclusion_OBJECTS =  \

@BUILDFORTRAN_FALSE@	../src/inclusion/np_inclusion.$(OBJEXT) \

@BUILDFORTRAN_FALSE@	../src/inclusion/inclusion.$(OBJEXT)

@BUILDFORTRAN_TRUE@am_inclusion_np_inclusion_OBJECTS =  \

@BUILDFORTRAN_TRUE@	../src/inclusion/np_inclusion.$(OBJEXT) \

@BUILDFORTRAN_TRUE@	../src/inclusion/inclusion.$(OBJEXT)

inclusion_np_inclusion_OBJECTS = $(am_inclusion_np_inclusion_OBJECTS)

inclusion_np_inclusion_LDADD = $(LDADD)

inclusion_np_inclusion_DEPENDENCIES = libnptm/libnptm.la \

	$(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \

	$(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1)

am__sphere_edfb_sph_SOURCES_DIST = ../src/sphere/edfb_sph.f

@BUILDFORTRAN_TRUE@am__objects_5 = ../src/sphere/edfb_sph.$(OBJEXT)

@BUILDFORTRAN_TRUE@am_sphere_edfb_sph_OBJECTS = $(am__objects_5)

am__maybe_remake_depfiles = depfiles

am__depfiles_remade = ../src/cluster/$(DEPDIR)/cluster.Po \

	../src/cluster/$(DEPDIR)/np_cluster.Po \

	../src/inclusion/$(DEPDIR)/inclusion.Po \

	../src/inclusion/$(DEPDIR)/np_inclusion.Po \

	../src/libnptm/$(DEPDIR)/Commons.Plo \

	../src/libnptm/$(DEPDIR)/Configuration.Plo \

	../src/libnptm/$(DEPDIR)/Parsers.Plo \

SOURCES = $(libnptm_libnptm_la_SOURCES) $(cluster_clu_SOURCES) \

	$(cluster_edfb_clu_SOURCES) $(cluster_np_cluster_SOURCES) \

	$(inclusion_edfb_inclu_SOURCES) $(inclusion_inclu_SOURCES) \

	$(sphere_edfb_sph_SOURCES) $(sphere_np_sphere_SOURCES) \

	$(sphere_sph_SOURCES) $(testing_test_TEDF_SOURCES) \

	$(testing_test_TTMS_SOURCES) $(trapping_frfme_SOURCES) \

	$(trapping_lffft_SOURCES) $(trapping_np_trapping_SOURCES)

	$(inclusion_np_inclusion_SOURCES) $(sphere_edfb_sph_SOURCES) \

	$(sphere_np_sphere_SOURCES) $(sphere_sph_SOURCES) \

	$(testing_test_TEDF_SOURCES) $(testing_test_TTMS_SOURCES) \

	$(trapping_frfme_SOURCES) $(trapping_lffft_SOURCES) \

	$(trapping_np_trapping_SOURCES)

DIST_SOURCES = $(libnptm_libnptm_la_SOURCES) \

	$(am__cluster_clu_SOURCES_DIST) \

	$(am__cluster_edfb_clu_SOURCES_DIST) \

	$(am__cluster_np_cluster_SOURCES_DIST) \

	$(am__inclusion_edfb_inclu_SOURCES_DIST) \

	$(am__inclusion_inclu_SOURCES_DIST) \

	$(am__inclusion_np_inclusion_SOURCES_DIST) \

	$(am__sphere_edfb_sph_SOURCES_DIST) \

	$(am__sphere_np_sphere_SOURCES_DIST) \

	$(am__sphere_sph_SOURCES_DIST) \

LDADD = libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}

lib_LTLIBRARIES = libnptm/libnptm.la

libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

@BUILDFORTRAN_FALSE@PROGS = cluster/np_cluster sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_TRUE@PROGS = cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_FALSE@PROGS = cluster/np_cluster inclusion/np_inclusion sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_TRUE@PROGS = cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu inclusion/np_inclusion sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_TRUE@EDFBCLUSOURCES = ../src/cluster/edfb_clu.f

@BUILDFORTRAN_TRUE@CLUSOURCES = ../src/cluster/clu.f

@BUILDFORTRAN_TRUE@EDFBINCLUSOURCES = ../src/inclusion/edfb_inclu.f

@BUILDFORTRAN_TRUE@trapping_lffft_SOURCES = $(LFFFTSOURCES)

@BUILDFORTRAN_FALSE@cluster_np_cluster_SOURCES = ../src/cluster/np_cluster.cpp ../src/cluster/cluster.cpp

@BUILDFORTRAN_TRUE@cluster_np_cluster_SOURCES = ../src/cluster/np_cluster.cpp ../src/cluster/cluster.cpp

@BUILDFORTRAN_FALSE@inclusion_np_inclusion_SOURCES = ../src/inclusion/np_inclusion.cpp ../src/inclusion/inclusion.cpp

@BUILDFORTRAN_TRUE@inclusion_np_inclusion_SOURCES = ../src/inclusion/np_inclusion.cpp ../src/inclusion/inclusion.cpp

@BUILDFORTRAN_FALSE@sphere_np_sphere_SOURCES = ../src/sphere/np_sphere.cpp ../src/sphere/sphere.cpp

@BUILDFORTRAN_TRUE@sphere_np_sphere_SOURCES = ../src/sphere/np_sphere.cpp ../src/sphere/sphere.cpp

@BUILDFORTRAN_FALSE@trapping_np_trapping_SOURCES = ../src/trapping/np_trapping.cpp ../src/trapping/cfrfme.cpp ../src/trapping/clffft.cpp

inclusion/inclu$(EXEEXT): $(inclusion_inclu_OBJECTS) $(inclusion_inclu_DEPENDENCIES) $(EXTRA_inclusion_inclu_DEPENDENCIES) inclusion/$(am__dirstamp)

	@rm -f inclusion/inclu$(EXEEXT)

	$(AM_V_F77LD)$(F77LINK) $(inclusion_inclu_OBJECTS) $(inclusion_inclu_LDADD) $(LIBS)

../src/inclusion/np_inclusion.$(OBJEXT):  \

	../src/inclusion/$(am__dirstamp) \

	../src/inclusion/$(DEPDIR)/$(am__dirstamp)

../src/inclusion/inclusion.$(OBJEXT):  \

	../src/inclusion/$(am__dirstamp) \

	../src/inclusion/$(DEPDIR)/$(am__dirstamp)

inclusion/np_inclusion$(EXEEXT): $(inclusion_np_inclusion_OBJECTS) $(inclusion_np_inclusion_DEPENDENCIES) $(EXTRA_inclusion_np_inclusion_DEPENDENCIES) inclusion/$(am__dirstamp)

	@rm -f inclusion/np_inclusion$(EXEEXT)

	$(AM_V_CXXLD)$(CXXLINK) $(inclusion_np_inclusion_OBJECTS) $(inclusion_np_inclusion_LDADD) $(LIBS)

../src/sphere/$(am__dirstamp):

	@$(MKDIR_P) ../src/sphere

	@: > ../src/sphere/$(am__dirstamp)

@AMDEP_TRUE@@am__include@ @am__quote@../src/cluster/$(DEPDIR)/cluster.Po@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/cluster/$(DEPDIR)/np_cluster.Po@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/inclusion/$(DEPDIR)/inclusion.Po@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/inclusion/$(DEPDIR)/np_inclusion.Po@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/Commons.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/Configuration.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/Parsers.Plo@am__quote@ # am--include-marker

	-rm -f $(am__CONFIG_DISTCLEAN_FILES)

		-rm -f ../src/cluster/$(DEPDIR)/cluster.Po

	-rm -f ../src/cluster/$(DEPDIR)/np_cluster.Po

	-rm -f ../src/inclusion/$(DEPDIR)/inclusion.Po

	-rm -f ../src/inclusion/$(DEPDIR)/np_inclusion.Po

	-rm -f ../src/libnptm/$(DEPDIR)/Commons.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/Configuration.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/Parsers.Plo

	-rm -rf $(top_srcdir)/autom4te.cache

		-rm -f ../src/cluster/$(DEPDIR)/cluster.Po

	-rm -f ../src/cluster/$(DEPDIR)/np_cluster.Po

	-rm -f ../src/inclusion/$(DEPDIR)/inclusion.Po

	-rm -f ../src/inclusion/$(DEPDIR)/np_inclusion.Po

	-rm -f ../src/libnptm/$(DEPDIR)/Commons.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/Configuration.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/Parsers.Plo

/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari

   This program is free software: you can redistribute it and/or modify

   it under the terms of the GNU General Public License as published by

   the Free Software Foundation, either version 3 of the License, or

   (at your option) any later version.

   This program is distributed in the hope that it will be useful,

   but WITHOUT ANY WARRANTY; without even the implied warranty of

   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

   GNU General Public License for more details.

   A copy of the GNU General Public License is distributed along with

   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.

 */

/*! \file inclusion.cpp

 *

 * \brief Implementation of the calculation for a sphere with inclusions.

 */

#include <chrono>

#include <cstdio>

#include <exception>

#include <fstream>

#include <hdf5.h>

#include <string>

#ifdef _OPENMP

#include <omp.h>

#endif

#ifdef USE_MPI

#ifndef MPI_VERSION

#include <mpi.h>

#endif

#endif

#ifdef USE_NVTX

#include <nvtx3/nvToolsExt.h>

#endif

#ifdef USE_MAGMA

#include <cuda_runtime.h>

#endif

#ifndef INCLUDE_TYPES_H_

#include "../include/types.h"

#endif

#ifndef INCLUDE_ERRORS_H_

#include "../include/errors.h"

#endif

#ifndef INCLUDE_LOGGING_H_

#include "../include/logging.h"

#endif

#ifndef INCLUDE_CONFIGURATION_H_

#include "../include/Configuration.h"

#endif

#ifndef INCLUDE_COMMONS_H_

#include "../include/Commons.h"

#endif

#ifndef INCLUDE_SPH_SUBS_H_

#include "../include/sph_subs.h"

#endif

#ifndef INCLUDE_CLU_SUBS_H_

#include "../include/clu_subs.h"

#endif

#ifndef INCLUDE_TRANSITIONMATRIX_H_

#include "../include/TransitionMatrix.h"

#endif

#ifndef INCLUDE_ALGEBRAIC_H_

#include "../include/algebraic.h"

#endif

#ifndef INCLUDE_LIST_H_

#include "../include/List.h"

#endif

#ifndef INCLUDE_FILE_IO_H_

#include "../include/file_io.h"

#endif

#ifndef INCLUDE_UTILS_H_

#include "../include/utils.h"

#endif

using namespace std;

/*! \brief C++ implementation of INCLU

 *

 * \param config_file: `string` Name of the configuration file.

 * \param data_file: `string` Name of the input data file.

 * \param output_path: `string` Directory to write the output files in.

 * \param mpidata: `mixMPI *` Pointer to an instance of MPI data settings.

 */

void inclusion(const string& config_file, const string& data_file, const string& output_path, const mixMPI *mpidata) {

  chrono::time_point<chrono::high_resolution_clock> t_start = chrono::high_resolution_clock::now();

  chrono::duration<double> elapsed;

  string message;

  string timing_name;

  FILE *timing_file;

  Logger *time_logger;

  if (mpidata->rank == 0) {

    timing_name = output_path + "/c_timing_mpi"+ to_string(mpidata->rank) +".log";

    timing_file = fopen(timing_name.c_str(), "w");

    time_logger = new Logger(LOG_DEBG, timing_file);

  }

  Logger *logger = new Logger(LOG_DEBG);

  int device_count = 0;

  //===========

  // Initialise MAGMA

  //===========

#ifdef USE_MAGMA

  cudaGetDeviceCount(&device_count);

  logger->log("DEBUG: Proc-" + to_string(mpidata->rank) + " found " + to_string(device_count) + " CUDA devices.\n", LOG_DEBG);

  logger->log("INFO: Process " + to_string(mpidata->rank) + " initializes MAGMA.\n");

  magma_int_t magma_result = magma_init();

  if (magma_result != MAGMA_SUCCESS) {

    logger->err("ERROR: Process " + to_string(mpidata->rank) + " failed to initilize MAGMA.\n");

    logger->err("PROC-" + to_string(mpidata->rank) + ": MAGMA error code " + to_string(magma_result) + "\n");

    if (mpidata->rank == 0) {

      fclose(timing_file);

      delete time_logger;

    }

    delete logger;

    return;

  }

#endif

  // end MAGMA initialisation

  //===========================

  // the following only happens on MPI process 0

  //===========================

  if (mpidata->rank == 0) {

#ifdef USE_NVTX

    nvtxRangePush("Set up");

#endif

    //=======================

    // Initialise sconf from configuration file

    //=======================

    logger->log("INFO: making legacy configuration...", LOG_INFO);

    ScattererConfiguration *sconf = NULL;

    try {

      sconf = ScattererConfiguration::from_dedfb(config_file);

    } catch(const OpenConfigurationFileException &ex) {

      logger->err("\nERROR: failed to open scatterer configuration file.\n");

      string message = "FILE: " + string(ex.what()) + "\n";

      logger->err(message);

      fclose(timing_file);

      delete time_logger;

      delete logger;

      return;

    }

    sconf->write_formatted(output_path + "/c_OEDFB");

    sconf->write_binary(output_path + "/c_TEDF");

    sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");

    // end logger initialisation

    //========================

    // Initialise gconf from configuration files

    //========================

    GeometryConfiguration *gconf = NULL;

    try {

      gconf = GeometryConfiguration::from_legacy(data_file);

    } catch (const OpenConfigurationFileException &ex) {

      logger->err("\nERROR: failed to open geometry configuration file.\n");

      string message = "FILE: " + string(ex.what()) + "\n";

      logger->err(message);

      if (sconf) delete sconf;

      fclose(timing_file);

      delete time_logger;

      delete logger;

      return;

    }

    logger->log(" done.\n", LOG_INFO);

    //end gconf initialisation

#ifdef USE_NVTX

    nvtxRangePop();

#endif

    int s_nsph = sconf->number_of_spheres;

    int nsph = gconf->number_of_spheres;

    // Sanity check on number of sphere consistency, should always be verified

    if (s_nsph == nsph) {

      // Shortcuts to variables stored in configuration objects

      ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);

      double wp = sconf->wp;

      // Open empty virtual ascii file for output

      VirtualAsciiFile *p_output = new VirtualAsciiFile();

      // for the time being, this is ok. When we can, add some logic in the sprintf calls that checks if a longer buffer would be needed, and in case expand it

      // in any case, replace all sprintf() with snprintf(), to avoid in any case writing more than the available buffer size

      char virtual_line[256];

      // Create and initialise pristine cid for MPI proc 0 and thread 0

      ClusterIterationData *cid = new ClusterIterationData(gconf, sconf, mpidata, device_count);

      const int ndi = cid->c4->nsph * cid->c4->nlim;

      np_int ndit = 2 * ndi;

      logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      //==========================

      // Write a block of info to the ascii output file

      //==========================

      sprintf(virtual_line, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");

      p_output->append_line(virtual_line);

#ifdef USE_ILP64

      sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#else

      sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#endif      

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");

      p_output->append_line(virtual_line);

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

	sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n",

		gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)

		);

	p_output->append_line(virtual_line);

      }

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

      p_output->append_line(virtual_line);

      sprintf(

	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",

	      p_scattering_angles->th, p_scattering_angles->thstp,

	      p_scattering_angles->thlst, p_scattering_angles->ths,

	      p_scattering_angles->thsstp, p_scattering_angles->thslst

	      );

      p_output->append_line(virtual_line);

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

      p_output->append_line(virtual_line);

      sprintf(

	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",

	      p_scattering_angles->ph, p_scattering_angles->phstp,

	      p_scattering_angles->phlst, p_scattering_angles->phs,

	      p_scattering_angles->phsstp, p_scattering_angles->phslst

	      );

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " READ(IR,*)JWTM\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " %5d\n", gconf->jwtm);

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)NSPHT\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(XIV(I),I=1,NXI)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)NSHL(I),ROS(I)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");

      p_output->append_line(virtual_line);

      sprintf(virtual_line, " \n");

      p_output->append_line(virtual_line);

      str(sconf, cid->c1, cid->c1ao, cid->c3, cid->c4, cid->c6);

      thdps(cid->c4->lm, cid->zpv);

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);

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

      p_output->append_line(virtual_line);

      // Create empty virtual binary file

      VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();

      string tppoan_name = output_path + "/c_TPPOAN";

#ifdef USE_MAGMA

      logger->log("INFO: using MAGMA calls.\n", LOG_INFO);

#elif defined USE_LAPACK

      logger->log("INFO: using LAPACK calls.\n", LOG_INFO);

#else

      logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);

#endif

      int iavm = gconf->iavm;

      int isam = gconf->isam;

      int inpol = gconf->in_pol;

      int nxi = sconf->number_of_scales;

      int nth = p_scattering_angles->nth;

      int nths = p_scattering_angles->nths;

      int nph = p_scattering_angles->nph;

      int nphs = p_scattering_angles->nphs;

      //========================

      // write a block of info to virtual binary file

      //========================

      vtppoanp->append_line(VirtualBinaryLine(iavm));

      vtppoanp->append_line(VirtualBinaryLine(isam));

      vtppoanp->append_line(VirtualBinaryLine(inpol));

      vtppoanp->append_line(VirtualBinaryLine(nxi));

      vtppoanp->append_line(VirtualBinaryLine(nth));

      vtppoanp->append_line(VirtualBinaryLine(nph));

      vtppoanp->append_line(VirtualBinaryLine(nths));

      vtppoanp->append_line(VirtualBinaryLine(nphs));

      if (sconf->idfc < 0) {

	cid->vk = cid->xip * cid->wn;

	sprintf(virtual_line, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);

	p_output->append_line(virtual_line);

	sprintf(virtual_line, " \n");

	p_output->append_line(virtual_line);

      }

      //==================================================

      // do the first outputs here, so that I open here the new files, afterwards I only append

      //==================================================

      p_output->write_to_disk(output_path + "/c_OINCLU");

      delete p_output;

      vtppoanp->write_to_disk(output_path + "/c_TPPOAN");

      delete vtppoanp;

    }

    delete sconf;

    delete gconf;

#ifdef USE_MAGMA

    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");

    magma_finalize();

#endif

    chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();

    elapsed = t_end - t_start;

    string message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";

    logger->log(message);

    logger->log("Finished: output written to " + output_path + "/c_OINCLU\n");

    time_logger->log(message);

    fclose(timing_file);

    delete time_logger;

    delete logger;

  } // end instructions block of MPI process 0

}

using namespace std;

extern void inclusion(const string& config_file, const string& data_file, const string& output_path, const mixMPI& mpidata);

extern void inclusion(const string& config_file, const string& data_file, const string& output_path, const mixMPI *mpidata);

/*! \brief Main program entry point.

 *