From 9e49e77025c819f96c8b126f292d0180466a4194 Mon Sep 17 00:00:00 2001
From: "Mulas, Giacomo" <gmulas@oa-cagliari.inaf.it>
Date: Wed, 29 May 2024 23:53:41 +0200
Subject: [PATCH] Virtual files are dumped to disk at each iteration in MPI
 processes, drastically reducing memory usage

---
 src/cluster/cluster.cpp | 2120 +++++++++++++++++++--------------------
 src/libnptm/Commons.cpp |    2 +
 src/libnptm/file_io.cpp |  108 +-
 3 files changed, 1080 insertions(+), 1150 deletions(-)

diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 44aff3db..17694a50 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -94,6 +94,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
  *  \param data_file: `string` Name of the input data file.
  *  \param output_path: `string` Directory to write the output files in.
  */
+
 void cluster(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;
@@ -102,6 +103,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
   FILE *timing_file = fopen(timing_name.c_str(), "w");
   Logger *time_logger = new Logger(LOG_DEBG, timing_file);
   Logger *logger = new Logger(LOG_DEBG);
+
+  //===========
+  // Initialise MAGMA
+  //===========
 #ifdef USE_MAGMA
   int device_count;
   cudaGetDeviceCount(&device_count);
@@ -117,11 +122,18 @@ void cluster(const string& config_file, const string& data_file, const string& o
     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 {
@@ -138,6 +150,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
     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);
@@ -152,25 +169,33 @@ void cluster(const string& config_file, const string& data_file, const string& o
       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;
-      //FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
+      // 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);
       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\0");
       p_output->append_line(virtual_line);
       sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n\0",
@@ -228,309 +253,260 @@ void cluster(const string& config_file, const string& data_file, const string& o
       double exri = sqrt(exdc);
       sprintf(virtual_line, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n\0", exri);
       p_output->append_line(virtual_line);
+
+      // Create empty virtual binary file
       VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();
-      VirtualBinaryLine *vbinline = NULL;
-      // fstream *tppoanp = new fstream;
-      // fstream &tppoan = *tppoanp;
       string tppoan_name = output_path + "/c_TPPOAN";
-      // tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
-      // if (tppoan.is_open()) {
 #ifdef USE_MAGMA
-	logger->log("INFO: using MAGMA calls.\n", LOG_INFO);
+      logger->log("INFO: using MAGMA calls.\n", LOG_INFO);
 #elif defined USE_LAPACK
-	logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
+      logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
 #else
-	logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
+      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;
-	// tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(iavm));
-	// tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(isam));
-	// tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(inpol));
-	// tppoan.write(reinterpret_cast<char *>(&nxi), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(nxi));
-	// tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(nth));
-	// tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(nph));
-	// tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
-	vtppoanp->append_line(VirtualBinaryLine(nths));
-	// tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
-	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\0", cid->vk);
-	  p_output->append_line(virtual_line);
-	  sprintf(virtual_line, " \n\0");
-	  p_output->append_line(virtual_line);
-	}
-	// do the first iteration on jxi488 separately, since it seems to be different from the others
-	int jxi488 = 1;
-	chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
+      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\0", cid->vk);
+	p_output->append_line(virtual_line);
+	sprintf(virtual_line, " \n\0");
+	p_output->append_line(virtual_line);
+      }
+
+      // do the first iteration on jxi488 separately, since it seems to be different from the others
+      int jxi488 = 1;
+      chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
-	nvtxRangePush("First iteration");
+      nvtxRangePush("First iteration");
 #endif
-	int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);
+      int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);
 #ifdef USE_NVTX
-	nvtxRangePop();
+      nvtxRangePop();
 #endif
-	chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
-	elapsed = start_iter_1 - t_start;
-	string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";
-	logger->log(message);
-	time_logger->log(message);
-	elapsed = end_iter_1 - start_iter_1;
-	message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
-	logger->log(message);
-	time_logger->log(message);
-	if (jer != 0) {
-	  // First loop failed. Halt the calculation.
-	  // tppoan.close();
-	  fclose(timing_file);
-	  //fclose(output);
-	  delete p_output;
-	  delete p_scattering_angles;
-	  delete cid;
-	  delete logger;
-	  delete time_logger;
-	  delete sconf;
-	  delete gconf;
-	  return;
-	}
-	// 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_OCLU");
-	// reallocate a new one (even if it would be more efficient to emty the existing one
+      chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
+      elapsed = start_iter_1 - t_start;
+      string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";
+      logger->log(message);
+      time_logger->log(message);
+      elapsed = end_iter_1 - start_iter_1;
+      message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
+      logger->log(message);
+      time_logger->log(message);
+      if (jer != 0) {
+	// First loop failed. Halt the calculation.
+	fclose(timing_file);
 	delete p_output;
-	p_output = new VirtualAsciiFile();
-	// now tppoan
-	vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
-	delete vtppoanp;
-	vtppoanp = new VirtualBinaryFile();
+	delete p_scattering_angles;
+	delete cid;
+	delete logger;
+	delete time_logger;
+	delete sconf;
+	delete gconf;
+	return;
+      }
 
-	// here go the calls that send data to be duplicated on other MPI processes from process 0 to others, using MPI broadcasts, but only if MPI is actually used
+      //==================================================
+      // 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_OCLU");
+      delete p_output;
+      vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
+      delete vtppoanp;
+
+      // here go the calls that send data to be duplicated on other MPI processes from process 0 to others, using MPI broadcasts, but only if MPI is actually used
 #ifdef MPI_VERSION
-	if (mpidata->mpirunning) {
-	  gconf->mpibcast(mpidata);
-	  sconf->mpibcast(mpidata);	    
-	  cid->mpibcast(mpidata);
-	  p_scattering_angles->mpibcast(mpidata);
-	}	
+      if (mpidata->mpirunning) {
+	gconf->mpibcast(mpidata);
+	sconf->mpibcast(mpidata);	    
+	cid->mpibcast(mpidata);
+	p_scattering_angles->mpibcast(mpidata);
+      }	
 #endif
-	// Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
-	int ompnumthreads = 1;
-	VirtualAsciiFile **p_outarray = NULL;
-	VirtualBinaryFile **vtppoanarray = NULL;
+      // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
+      int ompnumthreads = 1;
+      // this is for MPI process 0 (or even if we are not using MPI at all)
+      int myjxi488startoffset = 0;
+      int myMPIstride = ompnumthreads;
+      int myMPIblock = ompnumthreads;
+      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access the all later
+      VirtualAsciiFile **p_outarray = NULL;
+      VirtualBinaryFile **vtppoanarray = NULL;
 
 #ifdef USE_NVTX
-	nvtxRangePush("Parallel loop");
+      nvtxRangePush("Parallel loop");
 #endif
+
+      //===========================================
+      // open the OpenMP parallel context, so each thread can initialise its stuff
+      //===========================================
 #pragma omp parallel
-	{
-	  // Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
-	  int myompthread = 0;
+      {
+	// Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
+	int myompthread = 0;
+
 #ifdef _OPENMP
-	  // If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files
-	  myompthread = omp_get_thread_num();
-	  if (myompthread == 0) ompnumthreads = omp_get_num_threads();
+	// If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files
+	myompthread = omp_get_thread_num();
+	if (myompthread == 0) ompnumthreads = omp_get_num_threads();
 #endif
-	  if (myompthread == 0) {
-	    p_outarray = new VirtualAsciiFile*[ompnumthreads];
-	    vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
-	  }
-	  // To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
-	  ClusterIterationData *cid_2 = NULL;
-	  //FILE *output_2 = NULL;
-	  VirtualAsciiFile *p_output_2 = NULL;
-	  VirtualBinaryFile *vtppoanp_2 = NULL;
-	  // fstream *tppoanp_2 = NULL;
-	  // for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
-	  if (myompthread == 0) {
-	    cid_2 = cid;
-	    //output_2 = output;
-	    p_output_2 = p_output;
-	    // tppoanp_2 = tppoanp;
-	    vtppoanp_2 = vtppoanp;
-	  } else {
-	    // this is not thread 0, so do create fresh copies of all local variables
-	    cid_2 = new ClusterIterationData(*cid);
-	    //output_2 = fopen((output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), "w");
-	    p_output_2 = new VirtualAsciiFile();
-	    vtppoanp_2 = new VirtualBinaryFile();
-	    // tppoanp_2 = new fstream;
-	    // tppoanp_2->open((output_path + "/c_TPPOAN_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
+
+	if (myompthread == 0) {
+	  // Initialise some shared variables only on thread 0
+	  p_outarray = new VirtualAsciiFile*[ompnumthreads];
+	  vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
+	  myMPIblock = ompnumthreads;
+	  myMPIstride = myMPIblock;
+	}
+
+#ifdef MPI_VERSION
+	if (myompthread == 0) {
+	  if (mpidata->mpirunning) {
+	    // only go through this if MPI has been actually used
+	    for (int rr=1; rr<mpidata->nprocs; rr++) {
+	      // individually send their respective starting points to other MPI processes: they start immediately after the frequencies computed by previous processes so far
+	      int remotejxi488startoffset = myMPIstride;
+	      MPI_Send(&remotejxi488startoffset, 1, MPI_INT, rr, 3, MPI_COMM_WORLD);
+	      int remoteMPIblock;
+	      MPI_Recv(&remoteMPIblock, 1, MPI_INT, rr, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+	      // update myMPIstride to include the ones due to MPI process rr
+	      myMPIstride += remoteMPIblock;
+	    }
+	    // now I know the total myMPIstride, I can send it to all processes
+	    MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
 	  }
+	}
+#endif
+	// add an omp barrier to make sure that the global variables defined by thread 0 are known to all threads below this
+#pragma omp barrier
+
+	// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
+	ClusterIterationData *cid_2 = NULL;
+	VirtualAsciiFile *p_output_2 = NULL;
+	VirtualBinaryFile *vtppoanp_2 = NULL;
+	// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
+	if (myompthread == 0) {
+	  cid_2 = cid;
+	} else {
+	  // this is not thread 0, so do create fresh copies of all local variables
+	  cid_2 = new ClusterIterationData(*cid);
+	}
+	// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
+	if (myompthread==0) {
+	  logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
+	}
+#pragma omp barrier
+	// ok, now I can actually start the parallel calculations
+	for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {
+	  // the parallel loop over MPI processes covers a different set of indices for each thread
+#pragma omp barrier
+	  int myjxi488 = ixi488+myompthread;
+	  // each thread opens new virtual files and stores their pointers in the shared array
+	  p_output_2 = new VirtualAsciiFile();
+	  vtppoanp_2 = new VirtualBinaryFile();
+	  // each thread puts a copy of the pointers to its virtual files in the shared arrays
 	  p_outarray[myompthread] = p_output_2;
 	  vtppoanarray[myompthread] = vtppoanp_2;
-	  // fstream &tppoan_2 = *tppoanp_2;
-	  // make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
 #pragma omp barrier
-	  if (myompthread==0) {
-	    logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
-	  }
-	  // ok, now I can actually start the parallel calculations
-#pragma omp for
-	  for (jxi488 = cid_2->firstxi; jxi488 <= cid_2->lastxi; jxi488++) {
-	    int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
+
+	  // each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
+	  if (myjxi488<=cid_2->number_of_scales) {
+	    int jer = cluster_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
 	  }
+#pragma omp barrier
 
+#ifdef USE_NVTX
+	  nvtxRangePush("Output concatenation");
+#endif
 #pragma omp barrier
-	  // only threads different from 0 have to free local copies of variables and close local files
-	  if (myompthread != 0) {
-	    delete cid_2;
-	    //fclose(output_2);
-	    // p_output_2->write_to_disk(output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread));
-	    // delete p_output_2;
-	    // tppoanp_2->close();
-	    // delete tppoanp_2;
+	  // threads different from 0 append their virtual files to the one of thread 0, and delete them
+	  if (myompthread == 0) {
+	    for (int ti=1; ti<ompnumthreads; ti++) {
+	      p_outarray[0]->append(*(p_outarray[ti]));
+	      delete p_outarray[ti];
+	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
+	      delete vtppoanarray[ti];
+	    }
 	  }
 #pragma omp barrier
-	  {
-	    string message = "INFO: Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
-	    logger->log(message);
+	  //==============================================
+	  // Collect all virtual files on thread 0 of MPI process 0, and append them to disk
+	  //==============================================
+	  if (myompthread == 0) {
+	    // thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them
+	    p_outarray[0]->append_to_disk(output_path + "/c_OCLU");
+	    delete p_outarray[0];
+	    vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
+	    delete vtppoanarray[0];
+
+#ifdef MPI_VERSION
+	    if (mpidata->mpirunning) {
+	      // only go through this if MPI has been actually used
+	      for (int rr=1; rr<mpidata->nprocs; rr++) {
+		// get the data from process rr, creating a new virtual ascii file
+		VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
+		// append to disk and delete virtual ascii file
+		p_output->append_to_disk(output_path + "/c_OCLU");
+		delete p_output;
+		// get the data from process rr, creating a new virtual binary file
+		VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
+		// append to disk and delete virtual binary file
+		vtppoanp->append_to_disk(output_path + "/c_TPPOAN");
+		delete vtppoanp;
+		int test = MPI_Barrier(MPI_COMM_WORLD);
+	      }
+	    }
+#endif
 	  }
-	} // closes pragma omp parallel
+	  // end block writing to disk
 #ifdef USE_NVTX
-	nvtxRangePop();
-
-	nvtxRangePush("Output concatenation");
+	  nvtxRangePop();
 #endif
 #pragma omp barrier
-	{
-	  // thread 0 already wrote on global files, skip it and take care of appending the others
-	  for (int ti=1; ti<ompnumthreads; ti++) {
-	    p_outarray[0]->append(*(p_outarray[ti]));
-	    delete p_outarray[ti];
-	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
-	    delete vtppoanarray[ti];
-	  }
-	  p_outarray[0]->append_to_disk(output_path + "/c_OCLU");
-	  delete p_outarray[0];
+
+	} // close strided loop running on MPI processes
+	// delete the shared arrays I used to make available to thread 0 the virtual files of other threads
+#pragma omp barrier
+	if (myompthread == 0) {
 	  delete[] p_outarray;
-	  vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
-	  delete vtppoanarray[0];
 	  delete[] vtppoanarray;
-	  // for (int ri = 1; ri < ompnumthreads; ri++) {
-	    // string partial_file_name = output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(ri);
-	    // string message = "Copying ASCII output in MPI process " + to_string(mpidata->rank) + " of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
-	    // logger->log(message, LOG_DEBG);
-	    // FILE *partial_output = fopen(partial_file_name.c_str(), "r");
-	    // // have a look at the getline() or fgets() C functions to read one line at a time, instead of char by char, it would simplify things
-	    // char virtual_line[256];
-	    // int index = 0;
-	    // int c = fgetc(partial_output);
-	    // while (c != EOF) {
-	    //   virtual_line[index++] = c;
-	    //   if (c == '\n') {
-	    // 	virtual_line[index] = '\0';
-	    // 	p_output->append_line(virtual_line);
-	    // 	index = 0;
-	    //   }
-	    //   c = fgetc(partial_output);
-	    // }
-	    // fclose(partial_output);
-	    // remove(partial_file_name.c_str());
-	    // logger->log("done.\n", LOG_DEBG);
-	  //   string partial_file_name = output_path + "/c_TPPOAN_" + to_string(mpidata->rank) + "_" + to_string(ri);
-	  //   string message = "Copying binary output in MPI process " + to_string(mpidata->rank) + " of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
-	  //   logger->log(message, LOG_DEBG);
-	  //   fstream partial_tppoan;
-	  //   partial_tppoan.open(partial_file_name.c_str(), ios::in | ios::binary);
-	  //   partial_tppoan.seekg(0, ios::end);
-	  //   long buffer_size = partial_tppoan.tellg();
-	  //   char *binary_buffer = new char[buffer_size];
-	  //   partial_tppoan.seekg(0, ios::beg);
-	  //   partial_tppoan.read(binary_buffer, buffer_size);
-	  //   // tppoan.write(binary_buffer, buffer_size);
-	  //   partial_tppoan.close();
-	  //   delete[] binary_buffer;
-	  //   remove(partial_file_name.c_str());
-	  //   logger->log("done.\n", LOG_DEBG);
-	  // }
 	}
-	// here go the code to append the files written in MPI processes > 0 to the ones on MPI process 0
-#ifdef MPI_VERSION
-	if (mpidata->mpirunning) {
-	  // only go through this if MPI has been actually used
-	  for (int rr=1; rr<mpidata->nprocs; rr++) {
-	    // get the data from process rr
-	    VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
-	    p_output->append_to_disk(output_path + "/c_OCLU");
-	    delete p_output;
-	    VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
-	    vtppoanp->append_to_disk(output_path + "/c_TPPOAN_bis");
-	    delete vtppoanp;
-	    // // how many openmp threads did process rr use?
-	    // int remotethreads;
-	    // MPI_Recv(&remotethreads, 1, MPI_INT, rr, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	    // for (int ri=0; ri<remotethreads; ri++) {
-	      //   // first get the ASCII local file
-	      //   char *chunk_buffer;
-	      //   int chunk_buffer_size = -1;
-	      //   MPI_Recv(&chunk_buffer_size, 1, MPI_INT, rr, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	      //   while (chunk_buffer_size != 0) {
-	      // 	char *chunk_buffer = new char[chunk_buffer_size];
-	      // 	MPI_Recv(chunk_buffer, chunk_buffer_size, MPI_CHAR, rr, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	      // 	// fputs(chunk_buffer, output);
-	      // 	int index = 0, last_char = 0;
-	      // 	char c;
-	      // 	while (last_char < chunk_buffer_size) {
-	      // 	  c = chunk_buffer[last_char++];
-	      // 	  virtual_line[index++] = c;
-	      // 	  if (c == '\n') {
-	      // 	    virtual_line[index] = '\0';
-	      // 	    index = 0;
-	      // 	    p_output->append_line(virtual_line);
-	      // 	  }
-	      // 	}
-	      // 	delete[] chunk_buffer;
-	      // 	MPI_Recv(&chunk_buffer_size, 1, MPI_INT, rr, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	      //   }
-	      //   // if (ri<remotethreads-1) fprintf(output, "\n");
-
-	    //   // now get the binary local file
-	    //   long buffer_size = 0;
-	    //   // get the size of the buffer
-	    //   MPI_Recv(&buffer_size, 1, MPI_LONG, rr, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	    //   // allocate the bufer
-	    //   char *binary_buffer = new char[buffer_size];
-	    //   // actually receive the buffer
-	    //   MPI_Recv(binary_buffer, buffer_size, MPI_CHAR, rr, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-	    //   // we can write it to disk
-	    //   // tppoan.write(binary_buffer, buffer_size);
-	    //   delete[] binary_buffer;
-	    // }
-	  }
+	{
+	  string message = "INFO: Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
+	  logger->log(message);
 	}
-
-#endif
 #ifdef USE_NVTX
 	nvtxRangePop();
 #endif
-	// tppoanp->close();
-	// delete tppoanp;
-      // } else { // In case TPPOAN could not be opened. Should never happen.
-      // 	logger->err("\nERROR: failed to open TPPOAN file.\n");
-      // }
-      // fclose(output);
-      // p_output->write_to_disk(output_path + "/c_OCLU");
-      // delete p_output;
-      // Clean memory
-      delete cid;
+	delete cid_2;
+      }
       delete p_scattering_angles;
-    } else { // NSPH mismatch between geometry and scatterer configurations.
+    }    
+
+    else { // NSPH mismatch between geometry and scatterer configurations.
       throw UnrecognizedConfigurationException(
 					       "Inconsistent geometry and scatterer configurations."
 					       );
     }
+      
     delete sconf;
     delete gconf;
     chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();
@@ -539,8 +515,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
     logger->log(message);
     logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
     time_logger->log(message);
-  }
-
+  } // end instructions block of MPI process 0
+  
+    //===============================
+    // instruction block for MPI processes different from 0
+    //===============================
 #ifdef MPI_VERSION
   else {
     // here go the code for MPI processes other than 0
@@ -549,17 +528,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
     ScattererConfiguration *sconf = new ScattererConfiguration(mpidata);
     ClusterIterationData *cid = new ClusterIterationData(mpidata);
     ScatteringAngles *p_scattering_angles = new ScatteringAngles(mpidata);
-    // open separate files for other MPI processes
-    // File *output = fopen((output_path + "/c_OCLU_mpi"+ to_string(mpidata->rank)).c_str(), "w");
-    // fstream *tppoanp = new fstream;
-    // fstream &tppoan = *tppoanp;
-    // string tppoan_name = output_path + "/c_TPPOAN_mpi"+ to_string(mpidata->rank);
-    // tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
+    
     // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
     int ompnumthreads = 1;
     VirtualAsciiFile **p_outarray = NULL;
     VirtualBinaryFile **vtppoanarray = NULL;
-
+    int myjxi488startoffset;
+    int myMPIstride = ompnumthreads;
+    int myMPIblock = ompnumthreads;
+      
 #pragma omp parallel
     {
       // Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
@@ -570,965 +547,855 @@ void cluster(const string& config_file, const string& data_file, const string& o
       if (myompthread == 0) ompnumthreads = omp_get_num_threads();
 #endif
       if (myompthread == 0) {
+	// receive the start parameter from MPI process 0
+	MPI_Recv(&myjxi488startoffset, 1, MPI_INT, 0, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+	// send my number of omp threads to process 0
+	MPI_Send(&ompnumthreads, 1, MPI_INT, 0, 3, MPI_COMM_WORLD);
+	// receive myMPIstride sent by MPI process 0 to all processes
+	MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
+	// allocate virtualfiles for each thread
 	p_outarray = new VirtualAsciiFile*[ompnumthreads];
 	vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
       }
+#pragma omp barrier
       // To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
       ClusterIterationData *cid_2 = NULL;
-      //FILE *output_2 = NULL;
       VirtualAsciiFile *p_output_2 = NULL;
       VirtualBinaryFile *vtppoanp_2 = NULL;
       // fstream *tppoanp_2 = NULL;
       // for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
       if (myompthread == 0) {
 	cid_2 = cid;
-	// output_2 = output;
-	// tppoanp_2 = tppoanp;
       } else {
 	// this is not thread 0, so do create fresh copies of all local variables
 	cid_2 = new ClusterIterationData(*cid);
       }
-      // output_2 = fopen((output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), "w");
-      p_output_2 = new VirtualAsciiFile();
-      p_outarray[myompthread] = p_output_2;
-      vtppoanp_2 = new VirtualBinaryFile();
-      vtppoanarray[myompthread] = vtppoanp_2;
-      // tppoanp_2 = new fstream;
-      // tppoanp_2->open((output_path + "/c_TPPOAN_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
-      // fstream &tppoan_2 = *tppoanp_2;
       // make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
 #pragma omp barrier
-      if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
       // ok, now I can actually start the parallel calculations
-#pragma omp for
-      for (int jxi488 = cid_2->firstxi; jxi488 <= cid_2->lastxi; jxi488++) {
-	int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
-      }
-
+      for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {
+	// the parallel loop over MPI processes covers a different set of indices for each thread
 #pragma omp barrier
-      // only threads different from 0 have to free local copies of variables
-      if (myompthread != 0) {
-	delete cid_2;
-      }
-      // fclose(output_2);
-      // p_output_2->write_to_disk(output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread));
-      // delete p_output_2;
-      // tppoanp_2->close();
-      // delete tppoanp_2;
+	int myjxi488 = ixi488+myjxi488startoffset+myompthread;
+	// each thread opens new virtual files and stores their pointers in the shared array
+	p_output_2 = new VirtualAsciiFile();
+	vtppoanp_2 = new VirtualBinaryFile();
+	// each thread puts a copy of the pointers to its virtual files in the shared arrays
+	p_outarray[myompthread] = p_output_2;
+	vtppoanarray[myompthread] = vtppoanp_2;
 #pragma omp barrier
-      {
-	string message = "INFO: Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
-	logger->log(message);
-      }
-    } // closes pragma omp parallel
+	if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
+	// ok, now I can actually start the parallel calculations
+	// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
+	if (myjxi488<=cid_2->number_of_scales) {
+	  int jer = cluster_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
+	} // close the OMP parallel for loop
+
 #pragma omp barrier
-    {
-      for (int ti=1; ti<ompnumthreads; ti++) {
-	p_outarray[0]->append(*(p_outarray[ti]));
-	delete p_outarray[ti];
-	vtppoanarray[0]->append(*(vtppoanarray[ti]));
-	delete vtppoanarray[ti];
+	// threads different from 0 append their virtual files to the one of thread 0, and delete them
+	if (myompthread == 0) {
+	  for (int ti=1; ti<ompnumthreads; ti++) {
+	    p_outarray[0]->append(*(p_outarray[ti]));
+	    delete p_outarray[ti];
+	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
+	    delete vtppoanarray[ti];
+	  }
+	  // thread 0 sends the collected virtualfiles to thread 0 of MPI process 0, then deletes them
+	  for (int rr=1; rr<mpidata->nprocs; rr++) {
+	    if (rr == mpidata->rank) {
+	      p_outarray[0]->mpisend(mpidata);
+	      delete p_outarray[0];
+	      vtppoanarray[0]->mpisend(mpidata);
+	      delete vtppoanarray[0];
+	    }
+	    int test = MPI_Barrier(MPI_COMM_WORLD);
+	  }
+	}
+      } // close strided loop running on MPI processes
+      
+	// Clean memory
+#pragma omp barrier
+      if (myompthread == 0) {
+	delete[] p_outarray;
+	delete[] vtppoanarray;
       }
-      p_outarray[0]->mpisend(mpidata);
-      delete p_outarray[0];
-      delete[] p_outarray;
-      vtppoanarray[0]->mpisend(mpidata);
-      delete vtppoanarray[0];
-      delete[] vtppoanarray;
-      // // tell MPI process 0 how many threads we have on this process (not necessarily the same across all processes)
-      // MPI_Send(&ompnumthreads, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
-      // // reopen local files, send them all to MPI process 0
-      // for (int ri = 0; ri < ompnumthreads; ri++) {
-	// 	string partial_file_name = output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(ri);
-	// 	string message = "Copying ASCII output in MPI process " + to_string(mpidata->rank) + " of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
-	// 	logger->log(message, LOG_DEBG);
-	// 	fstream partial_output;
-	// 	partial_output.open(partial_file_name.c_str(), ios::in);
-	// 	partial_output.seekg(0, ios::end);
-	// 	const long partial_output_size = partial_output.tellg();
-	// 	partial_output.close();
-	// 	partial_output.open(partial_file_name.c_str(), ios::in);
-	// 	int chunk_buffer_size = 25165824; // Length of char array  with 24Mb size
-	// 	char *chunk_buffer = new char[chunk_buffer_size]();
-	// 	int full_chunks = (int)(partial_output_size / chunk_buffer_size);
-	// 	for (int fi = 0; fi < full_chunks; fi++) {
-	// 	  partial_output.read(chunk_buffer, chunk_buffer_size);
-	// 	  // If EOF is reached, do not send EOF character.
-	// 	  long ptr_position = partial_output.tellg();
-	// 	  if (ptr_position == partial_output_size) {
-	// 	    chunk_buffer[chunk_buffer_size - 1] = '\0';
-	// 	  }
-	// 	  // Send the size of the buffer that is being transmitted (Node-0 does not know whether it is full or not)
-	// 	  MPI_Send(&chunk_buffer_size, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
-	// 	  // Actually send the file contents to Node-0
-	// 	  MPI_Send(chunk_buffer, chunk_buffer_size, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
-	// 	}
-	// 	long ptr_position = partial_output.tellg();
-	// 	if (ptr_position < partial_output_size) {
-	// 	  // Send the last partial buffer
-	// 	  chunk_buffer_size = partial_output_size - ptr_position;
-	// 	  delete[] chunk_buffer;
-	// 	  chunk_buffer = new char[chunk_buffer_size];
-	// 	  partial_output.read(chunk_buffer, chunk_buffer_size);
-	// 	  // chunk_buffer[chunk_buffer_size - 1] = '\0';
-	// 	  // Send the size of the buffer that is being transmitted (Node-0 does not know whether it is full or not)
-	// 	  MPI_Send(&chunk_buffer_size, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
-	// 	  // Actually send the file contents to Node-0
-	// 	  MPI_Send(chunk_buffer, chunk_buffer_size, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
-	// 	}
-	// 	// Send a size 0 flag to inform Node-0 that the transmission is complete
-	// 	chunk_buffer_size = 0;
-	// 	MPI_Send(&chunk_buffer_size, 1, MPI_INT, 0, 1, MPI_COMM_WORLD);
-	// 	partial_output.close();
-	// 	delete[] chunk_buffer;
-	// 	remove(partial_file_name.c_str());
-	// 	logger->log("done.\n", LOG_DEBG);
-	
-      // 	string partial_file_name = output_path + "/c_TPPOAN_" + to_string(mpidata->rank) + "_" + to_string(ri);
-      // 	string message = "Copying binary output in MPI process " + to_string(mpidata->rank) + " of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
-      // 	logger->log(message, LOG_DEBG);
-      // 	fstream partial_tppoan;
-      // 	partial_tppoan.open(partial_file_name.c_str(), ios::in | ios::binary);
-      // 	partial_tppoan.seekg(0, ios::end);
-      // 	long buffer_size = partial_tppoan.tellg();
-      // 	char *binary_buffer = new char[buffer_size];
-      // 	partial_tppoan.seekg(0, ios::beg);
-      // 	partial_tppoan.read(binary_buffer, buffer_size);
-      // 	// tell MPI process 0 how large is the buffer
-      // 	MPI_Send(&buffer_size, 1, MPI_LONG, 0, 1, MPI_COMM_WORLD);
-      // 	// actually send the buffer
-      // 	MPI_Send(binary_buffer, buffer_size, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
-      // 	// // tppoan.write(binary_buffer, buffer_size);
-      // 	partial_tppoan.close();
-      // 	delete[] binary_buffer;
-      // 	remove(partial_file_name.c_str());
-      // 	logger->log("done.\n", LOG_DEBG);
-      // }
+      delete cid_2;
+
     }
-    // Clean memory
-    delete cid;
     delete p_scattering_angles;
     delete sconf;
     delete gconf;
-
-  }
 #endif
 #ifdef USE_MAGMA
-  logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");
-  magma_finalize();
+    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");
+    magma_finalize();
 #endif
-  fclose(timing_file);
-  delete time_logger;
-  delete logger;
+    fclose(timing_file);
+    delete time_logger;
+    delete logger;
+  }
 }
 
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp)
-{
-  int nxi = sconf->number_of_scales;
-  char virtual_line[256];
-  string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
-  Logger *logger = new Logger(LOG_DEBG);
-  logger->log(message);
-  chrono::duration<double> elapsed;
-  chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end;
-  int jer = 0;
-  int lcalc = 0;
-  int jaw = 1;
-  int li = gconf->li;
-  int le = gconf->le;
-  int lm = 0;
-  if (le > lm) lm = le;
-  if (li > lm) lm = li;
-  int nsph = gconf->number_of_spheres;
-  np_int mxndm = gconf->mxndm;
-  int iavm = gconf->iavm;
-  int inpol = gconf->in_pol;
-  int npnt = gconf->npnt;
-  int npntts = gconf->npntts;
-  int isam = gconf->iavm;
-  int jwtm = gconf->jwtm;
-  np_int ndit = 2 * nsph * cid->c4->nlim;
-  int isq, ibf;
+  int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp)
+  {
+    int nxi = sconf->number_of_scales;
+    char virtual_line[256];
+    string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
+    Logger *logger = new Logger(LOG_DEBG);
+    logger->log(message);
+    chrono::duration<double> elapsed;
+    chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end;
+    int jer = 0;
+    int lcalc = 0;
+    int jaw = 1;
+    int li = gconf->li;
+    int le = gconf->le;
+    int lm = 0;
+    if (le > lm) lm = le;
+    if (li > lm) lm = li;
+    int nsph = gconf->number_of_spheres;
+    np_int mxndm = gconf->mxndm;
+    int iavm = gconf->iavm;
+    int inpol = gconf->in_pol;
+    int npnt = gconf->npnt;
+    int npntts = gconf->npntts;
+    int isam = gconf->iavm;
+    int jwtm = gconf->jwtm;
+    np_int ndit = 2 * nsph * cid->c4->nlim;
+    int isq, ibf;
 
 #ifdef USE_NVTX
-  nvtxRangePush("Prepare matrix calculation");
+    nvtxRangePush("Prepare matrix calculation");
 #endif
-  sprintf(virtual_line, "========== JXI =%3d ====================\n\0", jxi488);
-  output->append_line(virtual_line);
-  double xi = sconf->get_scale(jxi488 - 1);
-  double exdc = sconf->exdc;
-  double exri = sqrt(exdc);
-  int idfc = (int)sconf->idfc;
-  double vkarg = 0.0;
-  if (idfc >= 0) {
-    cid->vk = xi * cid->wn;
-    vkarg = cid->vk;
-    sprintf(virtual_line, "  VK=%15.7lE, XI=%15.7lE\n\0", cid->vk, xi);
-    output->append_line(virtual_line);
-  } else {
-    vkarg = xi * cid->vk;
-    cid->sqsfi = 1.0 / (xi * xi);
-    sprintf(virtual_line, "  XI=%15.7lE\n\0", xi);
-    output->append_line(virtual_line);
-  }
-  hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4);
-  if (jer != 0) {
-    sprintf(virtual_line, "  STOP IN HJV\n\0");
+    sprintf(virtual_line, "========== JXI =%3d ====================\n\0", jxi488);
     output->append_line(virtual_line);
-    return jer;
-    // break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
-  }
-  for (int i132 = 1; i132 <= nsph; i132++) {
-    int iogi = cid->c1->iog[i132 - 1];
-    if (iogi != i132) {
-      for (int l123 = 1; l123 <= li; l123++) {
-	cid->c1->rmi[l123 - 1][i132 - 1] = cid->c1->rmi[l123 - 1][iogi - 1];
-	cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1];
-      } // l123 loop
+    double xi = sconf->get_scale(jxi488 - 1);
+    double exdc = sconf->exdc;
+    double exri = sqrt(exdc);
+    int idfc = (int)sconf->idfc;
+    double vkarg = 0.0;
+    if (idfc >= 0) {
+      cid->vk = xi * cid->wn;
+      vkarg = cid->vk;
+      sprintf(virtual_line, "  VK=%15.7lE, XI=%15.7lE\n\0", cid->vk, xi);
+      output->append_line(virtual_line);
     } else {
-      int nsh = cid->c1->nshl[i132 - 1];
-      int ici = (nsh + 1) / 2;
-      if (idfc == 0) {
-	for (int ic = 0; ic < ici; ic++)
-	  cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
+      vkarg = xi * cid->vk;
+      cid->sqsfi = 1.0 / (xi * xi);
+      sprintf(virtual_line, "  XI=%15.7lE\n\0", xi);
+      output->append_line(virtual_line);
+    }
+    hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4);
+    if (jer != 0) {
+      sprintf(virtual_line, "  STOP IN HJV\n\0");
+      output->append_line(virtual_line);
+      return jer;
+      // break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
+    }
+    for (int i132 = 1; i132 <= nsph; i132++) {
+      int iogi = cid->c1->iog[i132 - 1];
+      if (iogi != i132) {
+	for (int l123 = 1; l123 <= li; l123++) {
+	  cid->c1->rmi[l123 - 1][i132 - 1] = cid->c1->rmi[l123 - 1][iogi - 1];
+	  cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1];
+	} // l123 loop
       } else {
-	if (jxi488 == 1) {
+	int nsh = cid->c1->nshl[i132 - 1];
+	int ici = (nsh + 1) / 2;
+	if (idfc == 0) {
 	  for (int ic = 0; ic < ici; ic++)
-	    cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
+	    cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
+	} else {
+	  if (jxi488 == 1) {
+	    for (int ic = 0; ic < ici; ic++)
+	      cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
+	  }
+	}
+	if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;
+	dme(
+	    cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri,
+	    cid->c1, cid->c2, jer, lcalc, cid->arg
+	    );
+	if (jer != 0) {
+	  sprintf(virtual_line, "  STOP IN DME\n\0");
+	  output->append_line(virtual_line);
+	  return jer;
+	  //break;
 	}
       }
-      if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;
-      dme(
-	  cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri,
-	  cid->c1, cid->c2, jer, lcalc, cid->arg
-	  );
       if (jer != 0) {
-	sprintf(virtual_line, "  STOP IN DME\n\0");
-	output->append_line(virtual_line);
 	return jer;
 	//break;
       }
-    }
-    if (jer != 0) {
-      return jer;
-      //break;
-    }
-  } // i132 loop
+    } // i132 loop
 #ifdef USE_NVTX
-  nvtxRangePop();
+    nvtxRangePop();
 #endif
-  interval_start = chrono::high_resolution_clock::now();
+    interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
-  nvtxRangePush("Calculate inverted matrix");
+    nvtxRangePush("Calculate inverted matrix");
 #endif
-  cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
+    cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
 #ifdef USE_NVTX
-  nvtxRangePop();
+    nvtxRangePop();
 #endif
-  interval_end = chrono::high_resolution_clock::now();
-  elapsed = interval_end - interval_start;
-  message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
-  logger->log(message);
-  interval_start = chrono::high_resolution_clock::now();
+    interval_end = chrono::high_resolution_clock::now();
+    elapsed = interval_end - interval_start;
+    message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
+    logger->log(message);
+    interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
-  nvtxRangePush("Invert the matrix");
+    nvtxRangePush("Invert the matrix");
 #endif
-  invert_matrix(cid->am, ndit, jer, mxndm);
+    invert_matrix(cid->am, ndit, jer, mxndm);
 #ifdef USE_NVTX
-  nvtxRangePop();
+    nvtxRangePop();
 #endif
-  interval_end = chrono::high_resolution_clock::now();
-  elapsed = interval_end - interval_start;
-  message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
-  logger->log(message);
-  if (jer != 0) {
-    message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n";
-    logger->err(message);
-    return jer;
-    // break; // jxi488 loop: goes to memory clean
-  }
-  interval_start = chrono::high_resolution_clock::now();
+    interval_end = chrono::high_resolution_clock::now();
+    elapsed = interval_end - interval_start;
+    message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
+    logger->log(message);
+    if (jer != 0) {
+      message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n";
+      logger->err(message);
+      return jer;
+      // break; // jxi488 loop: goes to memory clean
+    }
+    interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
-  nvtxRangePush("Average calculation");
+    nvtxRangePush("Average calculation");
 #endif
-  ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
-  if (idfc >= 0) {
-    if (jxi488 == jwtm) {
-      int nlemt = 2 * cid->c4->nlem;
-      string ttms_name = output_path + "/c_TTMS.hd5";
-      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5");
-      ttms_name = output_path + "/c_TTMS";
-      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m);
-    }
-  }
-  // label 156: continue from here
-  if (inpol == 0) {
-    sprintf(virtual_line, "   LIN\n\0");
-    output->append_line(virtual_line);
-  } else { // label 158
-    sprintf(virtual_line, "  CIRC\n\0");
-    output->append_line(virtual_line);
-  }
-  // label 160
-  double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
-  double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
-  dcomplex s0 = 0.0 + 0.0 * I;
-  scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4);
-  double sqk = cid->vk * cid->vk * exdc;
-  aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps);
-  rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);
-  if (cid->c4->li != cid->c4->le) {
-    sprintf(virtual_line, "     SPHERES; LMX=LI\n\0");
-    output->append_line(virtual_line);
-  }
-  for (int i170 = 1; i170 <= nsph; i170++) {
-    if (cid->c1->iog[i170 - 1] >= i170) {
-      int i = i170 - 1;
-      double albeds = cid->c1->sscs[i] / cid->c1->sexs[i];
-      cid->c1->sqscs[i] *= cid->sqsfi;
-      cid->c1->sqabs[i] *= cid->sqsfi;
-      cid->c1->sqexs[i] *= cid->sqsfi;
-      sprintf(virtual_line, "     SPHERE %2d\n\0", i170);
-      output->append_line(virtual_line);
-      if (cid->c1->nshl[i] != 1) {
-	sprintf(virtual_line, "  SIZE=%15.7lE\n\0", cid->c2->vsz[i]);
-	output->append_line(virtual_line);
-      } else { // label 162
-	sprintf(virtual_line, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n\0", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i]));
-	output->append_line(virtual_line);
+    ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
+    if (idfc >= 0) {
+      if (jxi488 == jwtm) {
+	int nlemt = 2 * cid->c4->nlem;
+	string ttms_name = output_path + "/c_TTMS.hd5";
+	TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5");
+	ttms_name = output_path + "/c_TTMS";
+	TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m);
       }
-      // label 164
-      sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n\0");
-      output->append_line(virtual_line);
-      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i], albeds);
-      output->append_line(virtual_line);
-      sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n\0");
-      output->append_line(virtual_line);
-      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]);
-      output->append_line(virtual_line);
-      sprintf(virtual_line, "  FSAS=%15.7lE%15.7lE\n\0", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i]));
-      output->append_line(virtual_line);
-      double alamb = 2.0 * 3.141592653589793 / cid->vk;
-      sprintf(virtual_line, "INSERTION: CS_SPHERE  %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i]);
-      output->append_line(virtual_line);
-      csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcsv[i];
-      s0 = cid->c1->fsas[i] * exri;
-      qschu = imag(s0) * csch;
-      pschu = real(s0) * csch;
-      s0mag = cabs(s0) * cs0;
-      sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag);
-      output->append_line(virtual_line);
-      double rapr = cid->c1->sexs[i] - cid->gaps[i];
-      double cosav = cid->gaps[i] / cid->c1->sscs[i];
-      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
+    }
+    // label 156: continue from here
+    if (inpol == 0) {
+      sprintf(virtual_line, "   LIN\n\0");
       output->append_line(virtual_line);
-      sprintf(virtual_line, "  IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[0][i], cid->tqss[0][i]);
+    } else { // label 158
+      sprintf(virtual_line, "  CIRC\n\0");
       output->append_line(virtual_line);
-      sprintf(virtual_line, "  IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[1][i], cid->tqss[1][i]);
+    }
+    // label 160
+    double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
+    double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
+    dcomplex s0 = 0.0 + 0.0 * I;
+    scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4);
+    double sqk = cid->vk * cid->vk * exdc;
+    aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps);
+    rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);
+    if (cid->c4->li != cid->c4->le) {
+      sprintf(virtual_line, "     SPHERES; LMX=LI\n\0");
       output->append_line(virtual_line);
     }
-  } // i170 loop
-  sprintf(virtual_line, "  FSAT=%15.7lE%15.7lE\n\0", real(cid->c3->tfsas), imag(cid->c3->tfsas));
-  output->append_line(virtual_line);
-  csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs;
-  s0 = cid->c3->tfsas * exri;
-  qschu = imag(s0) * csch;
-  pschu = real(s0) * csch;
-  s0mag = cabs(s0) * cs0;
-  sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag);
-  output->append_line(virtual_line);
-  // tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double));
-  vtppoanp->append_line(VirtualBinaryLine(cid->vk));
-  pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4);
-  apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm);
+    for (int i170 = 1; i170 <= nsph; i170++) {
+      if (cid->c1->iog[i170 - 1] >= i170) {
+	int i = i170 - 1;
+	double albeds = cid->c1->sscs[i] / cid->c1->sexs[i];
+	cid->c1->sqscs[i] *= cid->sqsfi;
+	cid->c1->sqabs[i] *= cid->sqsfi;
+	cid->c1->sqexs[i] *= cid->sqsfi;
+	sprintf(virtual_line, "     SPHERE %2d\n\0", i170);
+	output->append_line(virtual_line);
+	if (cid->c1->nshl[i] != 1) {
+	  sprintf(virtual_line, "  SIZE=%15.7lE\n\0", cid->c2->vsz[i]);
+	  output->append_line(virtual_line);
+	} else { // label 162
+	  sprintf(virtual_line, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n\0", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i]));
+	  output->append_line(virtual_line);
+	}
+	// label 164
+	sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n\0");
+	output->append_line(virtual_line);
+	sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i], albeds);
+	output->append_line(virtual_line);
+	sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n\0");
+	output->append_line(virtual_line);
+	sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]);
+	output->append_line(virtual_line);
+	sprintf(virtual_line, "  FSAS=%15.7lE%15.7lE\n\0", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i]));
+	output->append_line(virtual_line);
+	double alamb = 2.0 * 3.141592653589793 / cid->vk;
+	sprintf(virtual_line, "INSERTION: CS_SPHERE  %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i]);
+	output->append_line(virtual_line);
+	csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcsv[i];
+	s0 = cid->c1->fsas[i] * exri;
+	qschu = imag(s0) * csch;
+	pschu = real(s0) * csch;
+	s0mag = cabs(s0) * cs0;
+	sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag);
+	output->append_line(virtual_line);
+	double rapr = cid->c1->sexs[i] - cid->gaps[i];
+	double cosav = cid->gaps[i] / cid->c1->sscs[i];
+	sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
+	output->append_line(virtual_line);
+	sprintf(virtual_line, "  IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[0][i], cid->tqss[0][i]);
+	output->append_line(virtual_line);
+	sprintf(virtual_line, "  IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[1][i], cid->tqss[1][i]);
+	output->append_line(virtual_line);
+      }
+    } // i170 loop
+    sprintf(virtual_line, "  FSAT=%15.7lE%15.7lE\n\0", real(cid->c3->tfsas), imag(cid->c3->tfsas));
+    output->append_line(virtual_line);
+    csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs;
+    s0 = cid->c3->tfsas * exri;
+    qschu = imag(s0) * csch;
+    pschu = real(s0) * csch;
+    s0mag = cabs(s0) * cs0;
+    sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag);
+    output->append_line(virtual_line);
+    // tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double));
+    vtppoanp->append_line(VirtualBinaryLine(cid->vk));
+    pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4);
+    apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm);
 #ifdef USE_NVTX
-  nvtxRangePop();
+    nvtxRangePop();
 #endif
-  interval_end = chrono::high_resolution_clock::now();
-  elapsed = interval_end - interval_start;
-  message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n\0";
-  logger->log(message);
-  interval_start = chrono::high_resolution_clock::now();
+    interval_end = chrono::high_resolution_clock::now();
+    elapsed = interval_end - interval_start;
+    message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n\0";
+    logger->log(message);
+    interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
-  nvtxRangePush("Angle loop");
+    nvtxRangePush("Angle loop");
 #endif
-  double th = sa->th;
-  for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
-    double ph = sa->ph;
-    double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
-    for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
-      int jw = 0;
-      if (sa->nk != 1 || jxi488 <= 1) {
-	upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
-	if (isam >= 0) {
-	  wmamp(
-		0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0,
-		nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1
-		);
-	  // label 182
-	  apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
-	  raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
-	  jw = 1;
-	}
-      } else { // label 180, NK == 1 AND JXI488 == 1
-	if (isam >= 0) {
-	  // label 182
-	  apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
-	  raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
-	  jw = 1;
-	}
-      }
-      // label 184
-      double thsl = sa->ths;
-      double phsph = 0.0;
-      for (int jths = 1; jths <= sa->nths; jths++) {
-	double ths = thsl;
-	int icspnv = 0;
-	if (isam > 1) ths += sa->thsca;
-	if (isam >= 1) {
-	  phsph = 0.0;
-	  if (ths < 0.0 || ths > 180.0) phsph = 180.0;
-	  if (ths < 0.0) ths *= -1.0;
-	  if (ths > 180.0) ths = 360.0 - ths;
-	  if (phsph != 0.0) icspnv = 1;
-	}
-	// label 186
-	double phs = sa->phs;
-	for (int jphs = 1; jphs <= sa->nphs; jphs++) {
-	  double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
-	  if (isam >= 1) {
-	    phs = sa->ph + phsph;
-	    if (phs > 360.0) phs -= 360.0;
-	  }
-	  // label 188
-	  bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
-	  if (!goto190) {
-	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
-	    if (isam >= 0)
-	      wmamp(
-		    2, costs, sints, cosps, sinps, inpol, cid->c4->le,
-		    0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1
-		    );
-	  }
-	  // label 190
-	  if (sa->nkks != 1 || jxi488 <= 1) {
-	    upvsp(
-		  cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,
-		  cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp
+    double th = sa->th;
+    for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
+      double ph = sa->ph;
+      double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
+      for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
+	int jw = 0;
+	if (sa->nk != 1 || jxi488 <= 1) {
+	  upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
+	  if (isam >= 0) {
+	    wmamp(
+		  0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0,
+		  nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1
 		  );
-	    if (isam < 0) {
-	      wmasp(
-		    cost, sint, cosp, sinp, costs, sints, cosps, sinps,
-		    cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le,
-		    0, nsph, cid->argi, cid->args, cid->c1
-		    );
-	    } else { // label 192
-	      for (int i193 = 0; i193 < 3; i193++) {
-		cid->up[i193] = cid->upmp[i193];
-		cid->un[i193] = cid->unmp[i193];
-		cid->ups[i193] = cid->upsmp[i193];
-		cid->uns[i193] = cid->unsmp[i193];
-	      }
-	    }
+	    // label 182
+	    apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
+	    raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
+	    jw = 1;
 	  }
-	  // label 194
-	  if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->c6);
-	  if (isam < 0) {
+	} else { // label 180, NK == 1 AND JXI488 == 1
+	  if (isam >= 0) {
+	    // label 182
 	    apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
 	    raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
 	    jw = 1;
 	  }
-	  // label 196
-	  // tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
-	  vtppoanp->append_line(VirtualBinaryLine(th));
-	  // tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
-	  vtppoanp->append_line(VirtualBinaryLine(ph));
-	  // tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
-	  vtppoanp->append_line(VirtualBinaryLine(ths));
-	  // tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
-	  vtppoanp->append_line(VirtualBinaryLine(phs));
-	  // tppoan.write(reinterpret_cast<char *>(&(cid->scan)), sizeof(double));
-	  vtppoanp->append_line(VirtualBinaryLine(cid->scan));
-	  if (jaw != 0) {
-	    jaw = 0;
-	    mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->cext);
-	    // We now have some implicit loops writing to binary
-	    for (int i = 0; i < 4; i++) {
-	      for (int j = 0; j < 4; j++) {
-		double value = cid->cext[i][j];
-		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-		vtppoanp->append_line(VirtualBinaryLine(value));
+	}
+	// label 184
+	double thsl = sa->ths;
+	double phsph = 0.0;
+	for (int jths = 1; jths <= sa->nths; jths++) {
+	  double ths = thsl;
+	  int icspnv = 0;
+	  if (isam > 1) ths += sa->thsca;
+	  if (isam >= 1) {
+	    phsph = 0.0;
+	    if (ths < 0.0 || ths > 180.0) phsph = 180.0;
+	    if (ths < 0.0) ths *= -1.0;
+	    if (ths > 180.0) ths = 360.0 - ths;
+	    if (phsph != 0.0) icspnv = 1;
+	  }
+	  // label 186
+	  double phs = sa->phs;
+	  for (int jphs = 1; jphs <= sa->nphs; jphs++) {
+	    double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
+	    if (isam >= 1) {
+	      phs = sa->ph + phsph;
+	      if (phs > 360.0) phs -= 360.0;
+	    }
+	    // label 188
+	    bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
+	    if (!goto190) {
+	      upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
+	      if (isam >= 0)
+		wmamp(
+		      2, costs, sints, cosps, sinps, inpol, cid->c4->le,
+		      0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1
+		      );
+	    }
+	    // label 190
+	    if (sa->nkks != 1 || jxi488 <= 1) {
+	      upvsp(
+		    cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,
+		    cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp
+		    );
+	      if (isam < 0) {
+		wmasp(
+		      cost, sint, cosp, sinp, costs, sints, cosps, sinps,
+		      cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le,
+		      0, nsph, cid->argi, cid->args, cid->c1
+		      );
+	      } else { // label 192
+		for (int i193 = 0; i193 < 3; i193++) {
+		  cid->up[i193] = cid->upmp[i193];
+		  cid->un[i193] = cid->unmp[i193];
+		  cid->ups[i193] = cid->upsmp[i193];
+		  cid->uns[i193] = cid->unsmp[i193];
+		}
 	      }
 	    }
-	    for (int i = 0; i < 2; i++) {
-	      double value = cid->c1ao->scscm[i];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = real(cid->c1ao->scscpm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1ao->scscpm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = cid->c1ao->ecscm[i];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = real(cid->c1ao->ecscpm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1ao->ecscpm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
+	    // label 194
+	    if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->c6);
+	    if (isam < 0) {
+	      apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
+	      raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
+	      jw = 1;
 	    }
-	    for (int i = 0; i < 3; i++) {
-	      for (int j = 0; j < 2; j++) {
-		double value = cid->gapm[i][j];
+	    // label 196
+	    // tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
+	    vtppoanp->append_line(VirtualBinaryLine(th));
+	    // tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
+	    vtppoanp->append_line(VirtualBinaryLine(ph));
+	    // tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
+	    vtppoanp->append_line(VirtualBinaryLine(ths));
+	    // tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
+	    vtppoanp->append_line(VirtualBinaryLine(phs));
+	    // tppoan.write(reinterpret_cast<char *>(&(cid->scan)), sizeof(double));
+	    vtppoanp->append_line(VirtualBinaryLine(cid->scan));
+	    if (jaw != 0) {
+	      jaw = 0;
+	      mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->cext);
+	      // We now have some implicit loops writing to binary
+	      for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+		  double value = cid->cext[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      for (int i = 0; i < 2; i++) {
+		double value = cid->c1ao->scscm[i];
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+		value = real(cid->c1ao->scscpm[i]);
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+		value = imag(cid->c1ao->scscpm[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = real(cid->gappm[i][j]);
+		value = cid->c1ao->ecscm[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->gappm[i][j]);
+		value = real(cid->c1ao->ecscpm[i]);
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+		value = imag(cid->c1ao->ecscpm[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
 	      }
-	    }
-	    sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm);
-	    output->append_line(virtual_line);
-	    int jlr = 2;
-	    for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
-	      int ipol = (ilr210 % 2 == 0) ? 1 : -1;
-	      if (ilr210 == 2) jlr = 1;
-	      double extsm = cid->c1ao->ecscm[ilr210 - 1];
-	      double qextm = extsm * cid->sqsfi / cid->c3->gcs;
-	      double extrm = extsm / cid->c3->ecs;
-	      double scasm = cid->c1ao->scscm[ilr210 - 1];
-	      double albdm = scasm / extsm;
-	      double qscam = scasm * cid->sqsfi / cid->c3->gcs;
-	      double scarm = scasm / cid->c3->scs;
-	      double abssm = extsm - scasm;
-	      double qabsm = abssm * cid->sqsfi / cid->c3->gcs;
-	      double absrm = abssm / cid->c3->acs;
-	      double acsecs = cid->c3->acs / cid->c3->ecs;
-	      if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
-	      dcomplex s0m = cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
-	      double qschum = imag(s0m) * csch;
-	      double pschum = real(s0m) * csch;
-	      double s0magm = cabs(s0m) * cs0;
-	      double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas);
-	      double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas);
-	      if (inpol == 0) {
-		sprintf(virtual_line, "   LIN %2d\n\0", ipol);
-		output->append_line(virtual_line);
-	      } else { // label 206
-		sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
-		output->append_line(virtual_line);
+	      for (int i = 0; i < 3; i++) {
+		for (int j = 0; j < 2; j++) {
+		  double value = cid->gapm[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = real(cid->gappm[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = imag(cid->gappm[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
 	      }
-	      // label 208
-	      sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0");
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", scasm, abssm, extsm, albdm);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0");
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", qscam, qabsm, qextm);
+	      sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm);
 	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0");
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", scarm, absrm, extrm);
-	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
-		      ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
-		      imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
-		      real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1ao->fsacm[jlr - 1][ilr210 - 1])
-		      );
-	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0",
-		      ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
-		      );
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschum, pschum, s0magm);
-	      output->append_line(virtual_line);
-	      double rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
-	      double cosav = cid->gapm[2][ilr210 - 1] / cid->c1ao->scscm[ilr210 - 1];
-	      double fz = rapr;
-	      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fk=%15.7lE\n\0", fz);
+	      int jlr = 2;
+	      for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
+		int ipol = (ilr210 % 2 == 0) ? 1 : -1;
+		if (ilr210 == 2) jlr = 1;
+		double extsm = cid->c1ao->ecscm[ilr210 - 1];
+		double qextm = extsm * cid->sqsfi / cid->c3->gcs;
+		double extrm = extsm / cid->c3->ecs;
+		double scasm = cid->c1ao->scscm[ilr210 - 1];
+		double albdm = scasm / extsm;
+		double qscam = scasm * cid->sqsfi / cid->c3->gcs;
+		double scarm = scasm / cid->c3->scs;
+		double abssm = extsm - scasm;
+		double qabsm = abssm * cid->sqsfi / cid->c3->gcs;
+		double absrm = abssm / cid->c3->acs;
+		double acsecs = cid->c3->acs / cid->c3->ecs;
+		if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
+		dcomplex s0m = cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
+		double qschum = imag(s0m) * csch;
+		double pschum = real(s0m) * csch;
+		double s0magm = cabs(s0m) * cs0;
+		double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas);
+		double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas);
+		if (inpol == 0) {
+		  sprintf(virtual_line, "   LIN %2d\n\0", ipol);
+		  output->append_line(virtual_line);
+		} else { // label 206
+		  sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
+		  output->append_line(virtual_line);
+		}
+		// label 208
+		sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0");
+		output->append_line(virtual_line);
+		sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", scasm, abssm, extsm, albdm);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0");
+		output->append_line(virtual_line);
+		sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", qscam, qabsm, qextm);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0");
+		output->append_line(virtual_line);
+		sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", scarm, absrm, extrm);
+		output->append_line(virtual_line);
+		sprintf(
+			virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
+			ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
+			imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
+			real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1ao->fsacm[jlr - 1][ilr210 - 1])
+			);
+		output->append_line(virtual_line);
+		sprintf(
+			virtual_line, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0",
+			ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
+			);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschum, pschum, s0magm);
+		output->append_line(virtual_line);
+		double rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
+		double cosav = cid->gapm[2][ilr210 - 1] / cid->c1ao->scscm[ilr210 - 1];
+		double fz = rapr;
+		sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, "  Fk=%15.7lE\n\0", fz);
+		output->append_line(virtual_line);
+		double alamb = 2.0 * 3.141592653589793 / cid->vk;
+		sprintf(virtual_line, "INSERTION: CSM_CLUSTER  %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, scasm, abssm, extsm);
+		output->append_line(virtual_line);
+	      } // ilr210 loop
+	      double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]);
+	      double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->c1ao->fsacm[0][0]);
+	      sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rmbrif);
 	      output->append_line(virtual_line);
-	      double alamb = 2.0 * 3.141592653589793 / cid->vk;
-	      sprintf(virtual_line, "INSERTION: CSM_CLUSTER  %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, scasm, abssm, extsm);
+	      sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rmdchr);
 	      output->append_line(virtual_line);
-	    } // ilr210 loop
-	    double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]);
-	    double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->c1ao->fsacm[0][0]);
-	    sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rmbrif);
-	    output->append_line(virtual_line);
-	    sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rmdchr);
-	    output->append_line(virtual_line);
-	  }
-	  // label 212
-	  sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n\0", jth486, jph484, jths, jphs);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n\0", th, ph, ths, phs);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  SCAND=%10.3lE\n\0", cid->scan);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  CFMP=%15.7lE, SFMP=%15.7lE\n\0", cid->cfmp, cid->sfmp);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  CFSP=%15.7lE, SFSP=%15.7lE\n\0", cid->cfsp, cid->sfsp);
-	  output->append_line(virtual_line);
-	  if (isam >= 0) {
-	    sprintf(virtual_line, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->un[0], cid->un[1], cid->un[2]);
-	    output->append_line(virtual_line);
-	    sprintf(virtual_line, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->uns[0], cid->uns[1], cid->uns[2]);
+	    }
+	    // label 212
+	    sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n\0", jth486, jph484, jths, jphs);
 	    output->append_line(virtual_line);
-	  } else { // label 214
-	    sprintf(virtual_line, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n\0", cid->un[0], cid->un[1], cid->un[2]);
+	    sprintf(virtual_line, "  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n\0", th, ph, ths, phs);
 	    output->append_line(virtual_line);
-	  }
-	  // label 220
-	  if (inpol == 0) {
-	    sprintf(virtual_line, "   LIN\n\0");
+	    sprintf(virtual_line, "  SCAND=%10.3lE\n\0", cid->scan);
 	    output->append_line(virtual_line);
-	  } else { // label 222
-	    sprintf(virtual_line, "  CIRC\n\0");
+	    sprintf(virtual_line, "  CFMP=%15.7lE, SFMP=%15.7lE\n\0", cid->cfmp, cid->sfmp);
 	    output->append_line(virtual_line);
-	  }
-	  // label 224
-	  scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4);
-	  if (cid->c4->li != cid->c4->le) {
-	    sprintf(virtual_line, "     SPHERES; LMX=MIN0(LI,LE)\n\0");
+	    sprintf(virtual_line, "  CFSP=%15.7lE, SFSP=%15.7lE\n\0", cid->cfsp, cid->sfsp);
 	    output->append_line(virtual_line);
-	  }
-	  for (int i226 = 1; i226 <= nsph; i226++) {
-	    if (cid->c1->iog[i226 - 1] >= i226) {
-	      sprintf(virtual_line, "     SPHERE %2d\n\0", i226);
+	    if (isam >= 0) {
+	      sprintf(virtual_line, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->un[0], cid->un[1], cid->un[2]);
 	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n\0",
-		      real(cid->c1->sas[i226 - 1][0][0]), imag(cid->c1->sas[i226 - 1][0][0]),
-		      real(cid->c1->sas[i226 - 1][1][0]), imag(cid->c1->sas[i226 - 1][1][0])
-		      );
+	      sprintf(virtual_line, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->uns[0], cid->uns[1], cid->uns[2]);
 	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n\0",
-		      real(cid->c1->sas[i226 - 1][0][1]), imag(cid->c1->sas[i226 - 1][0][1]),
-		      real(cid->c1->sas[i226 - 1][1][1]), imag(cid->c1->sas[i226 - 1][1][1])
-		      );
+	    } else { // label 214
+	      sprintf(virtual_line, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n\0", cid->un[0], cid->un[1], cid->un[2]);
+	      output->append_line(virtual_line);
+	    }
+	    // label 220
+	    if (inpol == 0) {
+	      sprintf(virtual_line, "   LIN\n\0");
 	      output->append_line(virtual_line);
-	      for (int j225 = 0; j225 < 16; j225++) {
-		cid->c1->vint[j225] = cid->c1->vints[i226 - 1][j225];
-	      } // j225 loop
-	      mmulc(cid->c1->vint, cid->cmullr, cid->cmul);
-	      sprintf(virtual_line, "  MULS\n\0");
+	    } else { // label 222
+	      sprintf(virtual_line, "  CIRC\n\0");
 	      output->append_line(virtual_line);
-	      for (int i1 = 0; i1 < 4; i1++) {
+	    }
+	    // label 224
+	    scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4);
+	    if (cid->c4->li != cid->c4->le) {
+	      sprintf(virtual_line, "     SPHERES; LMX=MIN0(LI,LE)\n\0");
+	      output->append_line(virtual_line);
+	    }
+	    for (int i226 = 1; i226 <= nsph; i226++) {
+	      if (cid->c1->iog[i226 - 1] >= i226) {
+		sprintf(virtual_line, "     SPHERE %2d\n\0", i226);
+		output->append_line(virtual_line);
 		sprintf(
-			virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
-			cid->cmul[i1][0], cid->cmul[i1][1], cid->cmul[i1][2], cid->cmul[i1][3]
+			virtual_line, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n\0",
+			real(cid->c1->sas[i226 - 1][0][0]), imag(cid->c1->sas[i226 - 1][0][0]),
+			real(cid->c1->sas[i226 - 1][1][0]), imag(cid->c1->sas[i226 - 1][1][0])
 			);
 		output->append_line(virtual_line);
-	      } // i1 loop
-	      sprintf(virtual_line, "  MULSLR\n\0");
-	      output->append_line(virtual_line);
-	      for (int i1 = 0; i1 < 4; i1++) {
 		sprintf(
-			virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
-			cid->cmullr[i1][0], cid->cmullr[i1][1], cid->cmullr[i1][2], cid->cmullr[i1][3]
+			virtual_line, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n\0",
+			real(cid->c1->sas[i226 - 1][0][1]), imag(cid->c1->sas[i226 - 1][0][1]),
+			real(cid->c1->sas[i226 - 1][1][1]), imag(cid->c1->sas[i226 - 1][1][1])
 			);
 		output->append_line(virtual_line);
-	      } // i1 loop
-	    }
-	  } // i226 loop
-	  sprintf(
-		  virtual_line, "  SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n\0",
-		  real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]),
-		  real(cid->c3->tsas[1][0]), imag(cid->c3->tsas[1][0])
-		  );
-	  output->append_line(virtual_line);
-	  sprintf(
-		  virtual_line, "  SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n\0",
-		  real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]),
-		  real(cid->c3->tsas[1][1]), imag(cid->c3->tsas[1][1])
-		  );
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "     CLUSTER\n\0");
-	  output->append_line(virtual_line);
-	  pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4);
-	  mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext);
-	  mmulc(cid->c1->vint, cid->cmullr, cid->cmul);
-	  if (jw != 0) {
-	    jw = 0;
-	    // Some implicit loops writing to binary.
-	    for (int i = 0; i < 4; i++) {
-	      for (int j = 0; j < 4; j++) {
-		double value = cid->cext[i][j];
-		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-		vtppoanp->append_line(VirtualBinaryLine(value));
+		for (int j225 = 0; j225 < 16; j225++) {
+		  cid->c1->vint[j225] = cid->c1->vints[i226 - 1][j225];
+		} // j225 loop
+		mmulc(cid->c1->vint, cid->cmullr, cid->cmul);
+		sprintf(virtual_line, "  MULS\n\0");
+		output->append_line(virtual_line);
+		for (int i1 = 0; i1 < 4; i1++) {
+		  sprintf(
+			  virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
+			  cid->cmul[i1][0], cid->cmul[i1][1], cid->cmul[i1][2], cid->cmul[i1][3]
+			  );
+		  output->append_line(virtual_line);
+		} // i1 loop
+		sprintf(virtual_line, "  MULSLR\n\0");
+		output->append_line(virtual_line);
+		for (int i1 = 0; i1 < 4; i1++) {
+		  sprintf(
+			  virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
+			  cid->cmullr[i1][0], cid->cmullr[i1][1], cid->cmullr[i1][2], cid->cmullr[i1][3]
+			  );
+		  output->append_line(virtual_line);
+		} // i1 loop
 	      }
-	    }
-	    for (int i = 0; i < 2; i++) {
-	      double value = cid->c1ao->scsc[i];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = real(cid->c1ao->scscp[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1ao->scscp[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = cid->c1ao->ecsc[i];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = real(cid->c1ao->ecscp[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1ao->ecscp[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	    }
-	    for (int i = 0; i < 3; i++) {
-	      for (int j = 0; j < 2; j++) {
-		double value = cid->gap[i][j];
+	    } // i226 loop
+	    sprintf(
+		    virtual_line, "  SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n\0",
+		    real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]),
+		    real(cid->c3->tsas[1][0]), imag(cid->c3->tsas[1][0])
+		    );
+	    output->append_line(virtual_line);
+	    sprintf(
+		    virtual_line, "  SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n\0",
+		    real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]),
+		    real(cid->c3->tsas[1][1]), imag(cid->c3->tsas[1][1])
+		    );
+	    output->append_line(virtual_line);
+	    sprintf(virtual_line, "     CLUSTER\n\0");
+	    output->append_line(virtual_line);
+	    pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4);
+	    mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext);
+	    mmulc(cid->c1->vint, cid->cmullr, cid->cmul);
+	    if (jw != 0) {
+	      jw = 0;
+	      // Some implicit loops writing to binary.
+	      for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+		  double value = cid->cext[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      for (int i = 0; i < 2; i++) {
+		double value = cid->c1ao->scsc[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = real(cid->gapp[i][j]);
+		value = real(cid->c1ao->scscp[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->gapp[i][j]);
+		value = imag(cid->c1ao->scscp[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-	      }
-	    }
-	    for (int i = 0; i < 2; i++) {
-	      for (int j = 0; j < 3; j++) {
-		double value = cid->tqce[i][j];
+		value = cid->c1ao->ecsc[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = real(cid->tqcpe[i][j]);
+		value = real(cid->c1ao->ecscp[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->tqcpe[i][j]);
+		value = imag(cid->c1ao->ecscp[i]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
 	      }
-	    }
-	    for (int i = 0; i < 2; i++) {
-	      for (int j = 0; j < 3; j++) {
-		double value = cid->tqcs[i][j];
+	      for (int i = 0; i < 3; i++) {
+		for (int j = 0; j < 2; j++) {
+		  double value = cid->gap[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = real(cid->gapp[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = imag(cid->gapp[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      for (int i = 0; i < 2; i++) {
+		for (int j = 0; j < 3; j++) {
+		  double value = cid->tqce[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = real(cid->tqcpe[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = imag(cid->tqcpe[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      for (int i = 0; i < 2; i++) {
+		for (int j = 0; j < 3; j++) {
+		  double value = cid->tqcs[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = real(cid->tqcps[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		  value = imag(cid->tqcps[i][j]);
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      for (int i = 0; i < 3; i++) {
+		double value = cid->u[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = real(cid->tqcps[i][j]);
+		value = cid->up[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->tqcps[i][j]);
+		value = cid->un[i];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
 	      }
 	    }
-	    for (int i = 0; i < 3; i++) {
-	      double value = cid->u[i];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = cid->up[i];
+	    // label 254
+	    for (int i = 0; i < 16; i++) {
+	      double value = real(cid->c1->vint[i]);
 	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = cid->un[i];
+	      value = imag(cid->c1->vint[i]);
 	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	      vtppoanp->append_line(VirtualBinaryLine(value));
 	    }
-	  }
-	  // label 254
-	  for (int i = 0; i < 16; i++) {
-	    double value = real(cid->c1->vint[i]);
-	    // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	    vtppoanp->append_line(VirtualBinaryLine(value));
-	    value = imag(cid->c1->vint[i]);
-	    // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	    vtppoanp->append_line(VirtualBinaryLine(value));
-	  }
-	  for (int i = 0; i < 4; i++) {
-	    for (int j = 0; j < 4; j++) {
-	      double value = cid->cmul[i][j];
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
+	    for (int i = 0; i < 4; i++) {
+	      for (int j = 0; j < 4; j++) {
+		double value = cid->cmul[i][j];
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+	      }
 	    }
-	  }
-	  int jlr = 2;
-	  for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
-	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;
-	    if (ilr290 == 2) jlr = 1;
-	    double extsec = cid->c1ao->ecsc[ilr290 - 1];
-	    double qext = extsec * cid->sqsfi / cid->c3->gcs;
-	    double extrat = extsec / cid->c3->ecs;
-	    double scasec = cid->c1ao->scsc[ilr290 - 1];
-	    double albedc = scasec / extsec;
-	    double qsca = scasec * cid->sqsfi / cid->c3->gcs;
-	    double scarat = scasec / cid->c3->scs;
-	    double abssec = extsec - scasec;
-	    double qabs = abssec * cid->sqsfi / cid->c3->gcs;
-	    double absrat = 1.0;
-	    double ratio = cid->c3->acs / cid->c3->ecs;
-	    if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs;
-	    s0 = cid->c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
-	    double qschu = imag(s0) * csch;
-	    double pschu = real(s0) * csch;
-	    s0mag = cabs(s0) * cs0;
-	    double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);
-	    double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);
-	    if (inpol == 0) {
-	      sprintf(virtual_line, "   LIN %2d\n\0", ipol);
+	    int jlr = 2;
+	    for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
+	      int ipol = (ilr290 % 2 == 0) ? 1 : -1;
+	      if (ilr290 == 2) jlr = 1;
+	      double extsec = cid->c1ao->ecsc[ilr290 - 1];
+	      double qext = extsec * cid->sqsfi / cid->c3->gcs;
+	      double extrat = extsec / cid->c3->ecs;
+	      double scasec = cid->c1ao->scsc[ilr290 - 1];
+	      double albedc = scasec / extsec;
+	      double qsca = scasec * cid->sqsfi / cid->c3->gcs;
+	      double scarat = scasec / cid->c3->scs;
+	      double abssec = extsec - scasec;
+	      double qabs = abssec * cid->sqsfi / cid->c3->gcs;
+	      double absrat = 1.0;
+	      double ratio = cid->c3->acs / cid->c3->ecs;
+	      if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs;
+	      s0 = cid->c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
+	      double qschu = imag(s0) * csch;
+	      double pschu = real(s0) * csch;
+	      s0mag = cabs(s0) * cs0;
+	      double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);
+	      double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);
+	      if (inpol == 0) {
+		sprintf(virtual_line, "   LIN %2d\n\0", ipol);
+		output->append_line(virtual_line);
+	      } else { // label 273
+		sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
+		output->append_line(virtual_line);
+	      }
+	      // label 275
+	      sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0");
 	      output->append_line(virtual_line);
-	    } else { // label 273
-	      sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
+	      sprintf(
+		      virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0",
+		      scasec, abssec, extsec, albedc
+		      );
 	      output->append_line(virtual_line);
-	    }
-	    // label 275
-	    sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0");
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0",
-		    scasec, abssec, extsec, albedc
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0");
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, " %14.7lE%15.7lE%15.7lE\n\0",
-		    qsca, qabs, qext
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0");
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, " %14.7lE%15.7lE%15.7lE\n\0",
-		    scarat, absrat, extrat
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
-		    ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]),
-		    jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->fsac[jlr - 1][ilr290 - 1])
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, "   SAC(%1d,%1d)=%15.7lE%15.7lE    SAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
-		    ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]),
-		    jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->sac[jlr - 1][ilr290 - 1])
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0",
-		    ilr290, ilr290, refinr, ilr290, ilr290, extcor
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0",
-		    qschu, pschu, s0mag
-		    );
-	    output->append_line(virtual_line);
-	    double alamb = 2.0 * 3.141592653589793 / cid->vk;
-	    if (ilr290 == 1) {
-	      sprintf(virtual_line, "INSERTION: CS1_CLUSTER  %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec);
-	    } else if (ilr290 == 2) {
-	      sprintf(virtual_line, "INSERTION: CS2_CLUSTER  %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec);
-	    }
-	    output->append_line(virtual_line);
-	    bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
-	    if (!goto190) {
-	      cid->gapv[0] = cid->gap[0][ilr290 - 1];
-	      cid->gapv[1] = cid->gap[1][ilr290 - 1];
-	      cid->gapv[2] = cid->gap[2][ilr290 - 1];
-	      double extins = cid->c1ao->ecsc[ilr290 - 1];
-	      double scatts = cid->c1ao->scsc[ilr290 - 1];
-	      double rapr, cosav, fp, fn, fk, fx, fy, fz;
-	      rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
-	      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
+	      sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0");
 	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n\0", fp, fn, fk);
+	      sprintf(
+		      virtual_line, " %14.7lE%15.7lE%15.7lE\n\0",
+		      qsca, qabs, qext
+		      );
 	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n\0", fx, fy, fz);
+	      sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0");
 	      output->append_line(virtual_line);
-	      cid->tqev[0] = cid->tqce[ilr290 - 1][0];
-	      cid->tqev[1] = cid->tqce[ilr290 - 1][1];
-	      cid->tqev[2] = cid->tqce[ilr290 - 1][2];
-	      cid->tqsv[0] = cid->tqcs[ilr290 - 1][0];
-	      cid->tqsv[1] = cid->tqcs[ilr290 - 1][1];
-	      cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
-	      double tep, ten, tek, tsp, tsn, tsk;
-	      tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk);
-	      sprintf(virtual_line, "   TQEl=%15.7lE,  TQEr=%15.7lE,  TQEk=%15.7lE\n\0", tep, ten, tek);
+	      sprintf(
+		      virtual_line, " %14.7lE%15.7lE%15.7lE\n\0",
+		      scarat, absrat, extrat
+		      );
 	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n\0", tsp, tsn, tsk);
+	      sprintf(
+		      virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
+		      ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]),
+		      jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->fsac[jlr - 1][ilr290 - 1])
+		      );
 	      output->append_line(virtual_line);
 	      sprintf(
-		      virtual_line, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n\0",
-		      cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2]
+		      virtual_line, "   SAC(%1d,%1d)=%15.7lE%15.7lE    SAC(%1d,%1d)=%15.7lE%15.7lE\n\0",
+		      ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]),
+		      jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->sac[jlr - 1][ilr290 - 1])
 		      );
 	      output->append_line(virtual_line);
 	      sprintf(
-		      virtual_line, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n\0",
-		      cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2]
+		      virtual_line, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0",
+		      ilr290, ilr290, refinr, ilr290, ilr290, extcor
 		      );
 	      output->append_line(virtual_line);
-	    }
-	  } //ilr290 loop
-	  double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]);
-	  double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->c1ao->fsac[0][0]);
-	  sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rbirif);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rdichr);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  MULC\n\0");
-	  output->append_line(virtual_line);
-	  for (int i = 0; i < 4; i++) {
-	    sprintf(
-		    virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
-		    cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
-		    );
-	    output->append_line(virtual_line);
-	  }
-	  sprintf(virtual_line, "  MULCLR\n\0");
-	  output->append_line(virtual_line);
-	  for (int i = 0; i < 4; i++) {
-	    sprintf(
-		    virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
-		    cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
-		    );
-	    output->append_line(virtual_line);
-	  }
-	  if (iavm != 0) {
-	    mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul);
-	    // Some implicit loops writing to binary.
-	    for (int i = 0; i < 16; i++) {
-	      double value;
-	      value = real(cid->c1ao->vintm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1ao->vintm[i]);
-	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-	      vtppoanp->append_line(VirtualBinaryLine(value));
-	    }
-	    for (int i = 0; i < 4; i++) {
-	      for (int j = 0; j < 4; j++) {
-		double value = cid->cmul[i][j];
-		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
-		vtppoanp->append_line(VirtualBinaryLine(value));
-	      }
-	    }
-	    sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm);
-	    output->append_line(virtual_line);
-	    if (inpol == 0) {
-	      sprintf(virtual_line, "   LIN\n\0");
+	      sprintf(
+		      virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0",
+		      qschu, pschu, s0mag
+		      );
 	      output->append_line(virtual_line);
-	    } else { // label 316
-	      sprintf(virtual_line, "  CIRC\n\0");
+	      double alamb = 2.0 * 3.141592653589793 / cid->vk;
+	      if (ilr290 == 1) {
+		sprintf(virtual_line, "INSERTION: CS1_CLUSTER  %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec);
+	      } else if (ilr290 == 2) {
+		sprintf(virtual_line, "INSERTION: CS2_CLUSTER  %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec);
+	      }
 	      output->append_line(virtual_line);
-	    }
-	    // label 318
+	      bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
+	      if (!goto190) {
+		cid->gapv[0] = cid->gap[0][ilr290 - 1];
+		cid->gapv[1] = cid->gap[1][ilr290 - 1];
+		cid->gapv[2] = cid->gap[2][ilr290 - 1];
+		double extins = cid->c1ao->ecsc[ilr290 - 1];
+		double scatts = cid->c1ao->scsc[ilr290 - 1];
+		double rapr, cosav, fp, fn, fk, fx, fy, fz;
+		rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
+		sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n\0", fp, fn, fk);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n\0", fx, fy, fz);
+		output->append_line(virtual_line);
+		cid->tqev[0] = cid->tqce[ilr290 - 1][0];
+		cid->tqev[1] = cid->tqce[ilr290 - 1][1];
+		cid->tqev[2] = cid->tqce[ilr290 - 1][2];
+		cid->tqsv[0] = cid->tqcs[ilr290 - 1][0];
+		cid->tqsv[1] = cid->tqcs[ilr290 - 1][1];
+		cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
+		double tep, ten, tek, tsp, tsn, tsk;
+		tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk);
+		sprintf(virtual_line, "   TQEl=%15.7lE,  TQEr=%15.7lE,  TQEk=%15.7lE\n\0", tep, ten, tek);
+		output->append_line(virtual_line);
+		sprintf(virtual_line, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n\0", tsp, tsn, tsk);
+		output->append_line(virtual_line);
+		sprintf(
+			virtual_line, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n\0",
+			cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2]
+			);
+		output->append_line(virtual_line);
+		sprintf(
+			virtual_line, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n\0",
+			cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2]
+			);
+		output->append_line(virtual_line);
+	      }
+	    } //ilr290 loop
+	    double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]);
+	    double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->c1ao->fsac[0][0]);
+	    sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rbirif);
+	    output->append_line(virtual_line);
+	    sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rdichr);
+	    output->append_line(virtual_line);
 	    sprintf(virtual_line, "  MULC\n\0");
 	    output->append_line(virtual_line);
 	    for (int i = 0; i < 4; i++) {
@@ -1547,27 +1414,74 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 		      );
 	      output->append_line(virtual_line);
 	    }
-	  }
-	  // label 420, continues jphs loop
-	  if (isam < 1) phs += sa->phsstp;
-	} // jphs loop, labeled 480
-	if (isam <= 1) thsl += sa->thsstp;
-      } // jths loop, labeled 482
-      ph += sa->phstp;
-    } // jph484 loop
-    th += sa->thstp;
-  } // jth486 loop
+	    if (iavm != 0) {
+	      mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul);
+	      // Some implicit loops writing to binary.
+	      for (int i = 0; i < 16; i++) {
+		double value;
+		value = real(cid->c1ao->vintm[i]);
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+		value = imag(cid->c1ao->vintm[i]);
+		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		vtppoanp->append_line(VirtualBinaryLine(value));
+	      }
+	      for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+		  double value = cid->cmul[i][j];
+		  // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+		  vtppoanp->append_line(VirtualBinaryLine(value));
+		}
+	      }
+	      sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm);
+	      output->append_line(virtual_line);
+	      if (inpol == 0) {
+		sprintf(virtual_line, "   LIN\n\0");
+		output->append_line(virtual_line);
+	      } else { // label 316
+		sprintf(virtual_line, "  CIRC\n\0");
+		output->append_line(virtual_line);
+	      }
+	      // label 318
+	      sprintf(virtual_line, "  MULC\n\0");
+	      output->append_line(virtual_line);
+	      for (int i = 0; i < 4; i++) {
+		sprintf(
+			virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
+			cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
+			);
+		output->append_line(virtual_line);
+	      }
+	      sprintf(virtual_line, "  MULCLR\n\0");
+	      output->append_line(virtual_line);
+	      for (int i = 0; i < 4; i++) {
+		sprintf(
+			virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n\0",
+			cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
+			);
+		output->append_line(virtual_line);
+	      }
+	    }
+	    // label 420, continues jphs loop
+	    if (isam < 1) phs += sa->phsstp;
+	  } // jphs loop, labeled 480
+	  if (isam <= 1) thsl += sa->thsstp;
+	} // jths loop, labeled 482
+	ph += sa->phstp;
+      } // jph484 loop
+      th += sa->thstp;
+    } // jth486 loop
 #ifdef USE_NVTX
-  nvtxRangePop();
+    nvtxRangePop();
 #endif
-  interval_end = chrono::high_resolution_clock::now();
-  elapsed = interval_end - interval_start;
-  message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
-  logger->log(message);
+    interval_end = chrono::high_resolution_clock::now();
+    elapsed = interval_end - interval_start;
+    message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
+    logger->log(message);
   
-  logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
+    logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
 
-  delete logger;
+    delete logger;
 
-  return jer;
-}
+    return jer;
+  }
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index f6a1c46a..1410bbea 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1184,6 +1184,8 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
   vk = rhs.vk;
   firstxi = rhs.firstxi;
   lastxi = rhs.lastxi;
+  xiblock = rhs.xiblock;
+  number_of_scales = rhs.number_of_scales;
 }
 
 #ifdef MPI_VERSION
diff --git a/src/libnptm/file_io.cpp b/src/libnptm/file_io.cpp
index 00afb199..28efaba4 100644
--- a/src/libnptm/file_io.cpp
+++ b/src/libnptm/file_io.cpp
@@ -277,9 +277,10 @@ VirtualAsciiFile::VirtualAsciiFile(const mixMPI *mpidata, int rr) {
 VirtualAsciiFile::~VirtualAsciiFile() {
   // is it necessary to pop them out one by one? isn't there the dedicated method of std::vector to clean the vector?
   // besides, shouldn't this be done anyway by the destructor of std:vector?
-  while (!_file_lines->size() > 0) {
-    _file_lines->pop_back();
-  }
+  _file_lines->clear();
+  // while (_file_lines->size() > 0) {
+  //   _file_lines->pop_back();
+  // }
   if (_file_lines != NULL) delete _file_lines;
 }
 
@@ -300,15 +301,17 @@ void VirtualAsciiFile::append_line(const string& line) {
 int VirtualAsciiFile::append_to_disk(const std::string& file_name) {
   // dump to disk the contents of the virtualasciifile, appending at the end of the given file_name
   int result = 0;
-  fstream output_file;
-  output_file.open(file_name, ios::app);
-  if (output_file.is_open()) {
-    for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-      output_file << *it;
+  if (_file_lines->size() > 0) {
+    fstream output_file;
+    output_file.open(file_name, ios::app);
+    if (output_file.is_open()) {
+      for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+	output_file << *it;
+      }
+      output_file.close();
+    } else {
+      result = 1;
     }
-    output_file.close();
-  } else {
-    result = 1;
   }
   return result;
 }
@@ -334,15 +337,17 @@ int VirtualAsciiFile::insert(int32_t position, VirtualAsciiFile& rhs, int32_t st
 int VirtualAsciiFile::write_to_disk(const std::string& file_name) {
   // dump to disk the contents of the virtualasciifile, replacing the given file_name
   int result = 0;
-  fstream output_file;
-  output_file.open(file_name, ios::out);
-  if (output_file.is_open()) {
-    for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-      output_file << *it;
+  if (_file_lines->size() > 0) {
+    fstream output_file;
+    output_file.open(file_name, ios::out);
+    if (output_file.is_open()) {
+      for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+	output_file << *it;
+      }
+      output_file.close();
+    } else {
+      result = 1;
     }
-    output_file.close();
-  } else {
-    result = 1;
   }
   return result;
 }
@@ -354,13 +359,15 @@ void VirtualAsciiFile::mpisend(const mixMPI *mpidata) {
   int32_t num_lines =  _file_lines->size();
   MPI_Send(&num_lines, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
   // now loop over data to send
-  int32_t mysize;
-  for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-    // send the length of each string
-    mysize = (int32_t) it->size();
-    MPI_Send(&mysize, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    // send the string itself
-    MPI_Send(it->c_str(), mysize+1, MPI_CHAR, 0, 10, MPI_COMM_WORLD);
+  if (num_lines>0) {
+    int32_t mysize;
+    for (vector<string>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+      // send the length of each string
+      mysize = (int32_t) it->size();
+      MPI_Send(&mysize, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      // send the string itself
+      MPI_Send(it->c_str(), mysize+1, MPI_CHAR, 0, 10, MPI_COMM_WORLD);
+    }
   }
 }
 #endif
@@ -483,9 +490,10 @@ VirtualBinaryFile::~VirtualBinaryFile() {
   // for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
   //   delete it;
   // }
-  while (!_file_lines->size() > 0) {
-    _file_lines->pop_back();
-  }
+  _file_lines->clear();
+  // while (_file_lines->size() > 0) {
+  //   _file_lines->pop_back();
+  // }
   if (_file_lines != NULL) delete _file_lines;
 }
 
@@ -506,15 +514,17 @@ void VirtualBinaryFile::append_line(const VirtualBinaryLine& line) {
 int VirtualBinaryFile::append_to_disk(const std::string& file_name) {
   // dump to disk the contents of the virtualasciifile, appending at the end of the given file_name
   int result = 0;
-  fstream output_file;
-  output_file.open(file_name, ios::app | ios::binary);
-  if (output_file.is_open()) {
-    for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-      output_file.write(it->_data_pointer, it->_data_size);
+  if (_file_lines->size() > 0) {
+    fstream output_file;
+    output_file.open(file_name, ios::app | ios::binary);
+    if (output_file.is_open()) {
+      for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+	output_file.write(it->_data_pointer, it->_data_size);
+      }
+      output_file.close();
+    } else {
+      result = 1;
     }
-    output_file.close();
-  } else {
-    result = 1;
   }
   return result;
 }
@@ -541,15 +551,17 @@ int VirtualBinaryFile::append_to_disk(const std::string& file_name) {
 int VirtualBinaryFile::write_to_disk(const std::string& file_name) {
   // dump to disk the contents of the virtualasciifile, replacing the given file_name
   int result = 0;
-  fstream output_file;
-  output_file.open(file_name, ios::out | ios::binary);
-  if (output_file.is_open()) {
-    for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-      output_file.write(it->_data_pointer, it->_data_size);
+  if (_file_lines->size() > 0) {
+    fstream output_file;
+    output_file.open(file_name, ios::out | ios::binary);
+    if (output_file.is_open()) {
+      for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+	output_file.write(it->_data_pointer, it->_data_size);
+      }
+      output_file.close();
+    } else {
+      result = 1;
     }
-    output_file.close();
-  } else {
-    result = 1;
   }
   return result;
 }
@@ -561,8 +573,10 @@ void VirtualBinaryFile::mpisend(const mixMPI *mpidata) {
   int32_t num_lines =  _file_lines->size();
   MPI_Send(&num_lines, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
   // now loop over data to send
-  for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
-    it->mpisend(mpidata);
+  if (num_lines>0) {
+    for (vector<VirtualBinaryLine>::iterator it = _file_lines->begin(); it != _file_lines->end(); ++it) {
+      it->mpisend(mpidata);
+    }
   }
 }
 #endif
-- 
GitLab