diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 3b722c4051c10e25d34e2ee44f603af6edb7d08c..1d28b0a593efa00f77ab2ee27627d08aa84eb4e9 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -91,6 +91,10 @@
 #include "../include/utils.h"
 #endif
 
+#ifndef INCLUDE_OUTPUTS_H_
+#include "../include/outputs.h"
+#endif
+
 using namespace std;
 
 // >>> InclusionIterationData header <<< //
@@ -816,11 +820,11 @@ InclusionIterationData::~InclusionIterationData() {
  *  \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
  *  \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object.
  *  \param cid: `InclusionIterationData *` Pointer to an `InclusionIterationData` object.
- *  \param output: `VirtualAsciiFile *` Pointer to a `VirtualAsciiFile` object.
+ *  \param oi: `InclusionOutputInfo *` Pointer to an `InclusionOutputInfo` object.
  *  \param output_path: `const string &` Path to the output directory.
  *  \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object.
  */
-int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp);
+int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp);
 
 /*! \brief C++ implementation of INCLU
  *
@@ -929,81 +933,17 @@ void inclusion(const string& config_file, const string& data_file, const string&
       ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
       double wp = sconf->wp;
       // Open empty virtual ascii file for output
-      VirtualAsciiFile *p_output = new VirtualAsciiFile();
-      char virtual_line[256];
+      InclusionOutputInfo *p_output = new InclusionOutputInfo(sconf, gconf, mpidata);
       InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count);
       const np_int ndi = cid->c1->nsph * cid->c1->nlim;
       const np_int ndit = 2 * ndi;
       logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
       time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
       
-      //==========================
-      // Write a block of info to the ascii output file
-      //==========================
-      sprintf(virtual_line, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
-      p_output->append_line(virtual_line);
-#ifdef USE_ILP64
-      sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
-	      nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
-	      gconf->npntts, gconf->iavm, gconf->iavm
-	      );
-#else
-      sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
-	      nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
-	      gconf->npntts, gconf->iavm, gconf->iavm
-	      );
-#endif // USE_ILP64
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
-      p_output->append_line(virtual_line);
-      for (int ri = 0; ri < nsph; ri++) {
-	sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n",
-		gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
-		);
-	p_output->append_line(virtual_line);
-      }
-      sprintf(virtual_line, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
-      p_output->append_line(virtual_line);
-      sprintf(
-	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
-	      p_scattering_angles->th, p_scattering_angles->thstp,
-	      p_scattering_angles->thlst, p_scattering_angles->ths,
-	      p_scattering_angles->thsstp, p_scattering_angles->thslst
-	      );
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
-      p_output->append_line(virtual_line);
-      sprintf(
-	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
-	      p_scattering_angles->ph, p_scattering_angles->phstp,
-	      p_scattering_angles->phlst, p_scattering_angles->phs,
-	      p_scattering_angles->phsstp, p_scattering_angles->phslst
-	      );
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, " READ(IR,*)JWTM\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, " %5d\n", gconf->jwtm);
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)NSPHT\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)(XIV(I),I=1,NXI)\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)NSHL(I),ROS(I)\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
-      p_output->append_line(virtual_line);
-      sprintf(virtual_line, " \n");
-      p_output->append_line(virtual_line);
       instr(sconf, cid->c1);
       thdps(cid->c1->lm, cid->zpv);
       double exdc = sconf->exdc;
       double exri = sqrt(exdc);
-      sprintf(virtual_line, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
-      p_output->append_line(virtual_line);
 
       // Create an empty bynary file
       VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();
@@ -1037,10 +977,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
       vtppoanp->append_line(VirtualBinaryLine(nphs));
       if (sconf->idfc < 0) {
 	cid->vk = cid->xip * cid->wn;
-	sprintf(virtual_line, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);
-	p_output->append_line(virtual_line);
-	sprintf(virtual_line, " \n");
-	p_output->append_line(virtual_line);
+	p_output->vec_vk[0] = cid->vk;
       }
       
       // do the first iteration on jxi488 separately, since it seems to be different from the others
@@ -1092,8 +1029,6 @@ void inclusion(const string& config_file, const string& data_file, const string&
       //==================================================
       // do the first outputs here, so that I open here the new files, afterwards I only append
       //==================================================
-      p_output->write_to_disk(output_path + "/c_OINCLU");
-      delete p_output;
       vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
       delete vtppoanp;
       
@@ -1113,7 +1048,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
       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 them all later
-      VirtualAsciiFile **p_outarray = NULL;
+      InclusionOutputInfo **p_outarray = NULL;
       VirtualBinaryFile **vtppoanarray = NULL;
 
 #ifdef USE_NVTX
@@ -1136,7 +1071,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
 
 	if (myompthread == 0) {
 	  // Initialise some shared variables only on thread 0
-	  p_outarray = new VirtualAsciiFile*[ompnumthreads];
+	  p_outarray = new InclusionOutputInfo*[ompnumthreads];
 	  vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
 	  myMPIblock = ompnumthreads;
 	  myMPIstride = myMPIblock;
@@ -1165,11 +1100,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
 
 	// 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
 	InclusionIterationData *cid_2 = NULL;
-	VirtualAsciiFile *p_output_2 = NULL;
+	InclusionOutputInfo *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;
+	  // OMP thread 0 of MPI process 0 holds the pointer to the full output structure
+	  p_output_2 = p_output;
+	  p_outarray[0] = p_output_2;
 	} else {
 	  // this is not thread 0, so do create fresh copies of all local variables
 	  cid_2 = new InclusionIterationData(*cid);
@@ -1185,16 +1123,24 @@ void inclusion(const string& config_file, const string& data_file, const string&
 #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;
 #pragma omp barrier
 
 	  // 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) {
+	    if (myompthread > 0) {
+	      // UPDATE: non-0 threads need to allocate memory for one scale at a time.
+	      p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
+	      p_outarray[myompthread] = p_output_2;
+	    }
 	    int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
+	  } else {
+	    if (myompthread > 0) {
+	      // If there is no input for this thread, set output pointer to NULL.
+	      p_outarray[myompthread] = NULL;
+	    }
 	  }
 #pragma omp barrier
 
@@ -1205,8 +1151,11 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  // 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];
+	      if (p_outarray[ti] != NULL) {
+		p_outarray[0]->insert(*(p_outarray[ti]));
+		delete p_outarray[ti];
+		p_outarray[ti] = NULL;
+	      }
 	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	      delete vtppoanarray[ti];
 	    }
@@ -1217,8 +1166,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  //==============================================
 	  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_OINCLU");
-	    delete p_outarray[0];
+	    // p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");
+	    // delete p_outarray[0];
 	    vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
 	    delete vtppoanarray[0];
 
@@ -1226,11 +1175,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	    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 by receiving it in total memory structure
+		p_outarray[0]->mpireceive(mpidata, rr);
 		// get the data from process rr, creating a new virtual ascii file
-		VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
+		// VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
 		// append to disk and delete virtual ascii file
-		p_output->append_to_disk(output_path + "/c_OINCLU");
-		delete p_output;
+		// p_output->append_to_disk(output_path + "/c_OINCLU");
+		// 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
@@ -1264,12 +1216,15 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	delete cid_2;
       }
       delete p_scattering_angles;
-    }    
+      p_output->write(output_path + "/c_OINCLU.hd5", "HDF5");
+      p_output->write(output_path + "/c_OINCLU", "LEGACY");
+      delete p_output;
+    } // closes s_nsph == nsph check
 	  
     else { // Sphere number inconsistency error.
       throw UnrecognizedConfigurationException(
-					       "Inconsistent geometry and scatterer configurations."
-					       );
+	"Inconsistent geometry and scatterer configurations."
+      );
     }
 
     delete sconf;
@@ -1302,7 +1257,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
 
     // 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;
+    InclusionOutputInfo **p_outarray = NULL;
     VirtualBinaryFile **vtppoanarray = NULL;
     int myjxi488startoffset;
     int myMPIstride = ompnumthreads;
@@ -1325,13 +1280,13 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	// receive myMPIstride sent by MPI process 0 to all processes
 	MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
 	// allocate virtual files for each thread
-	p_outarray = new VirtualAsciiFile*[ompnumthreads];
+	p_outarray = new InclusionOutputInfo*[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
       InclusionIterationData *cid_2 = NULL;
-      VirtualAsciiFile *p_output_2 = NULL;
+      InclusionOutputInfo *p_output_2 = NULL;
       VirtualBinaryFile *vtppoanp_2 = NULL;
       // PLACEHOLDER
       // 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
@@ -1349,25 +1304,41 @@ void inclusion(const string& config_file, const string& data_file, const string&
 #pragma omp barrier
 	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
 	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) {
+	  if (myompthread > 0) {
+	    // UPDATE: non-0 threads need to allocate memory for one scale at a time.
+	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
+	    p_outarray[myompthread] = p_output_2;
+	  } else {
+	    // Thread 0 of non-zero MPI processes needs to allocate memory for the
+	    // output of all threads.
+	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
+	    p_outarray[0] = p_output_2;
+	  }
 	  int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
-	} // close the OMP parallel for loop
+	} else {
+	  if (myompthread > 0) {
+	    // If there is no input for this thread, set the output pointer to NULL.
+	    p_outarray[myompthread] = NULL;
+	  }	  
+	}
 
 #pragma omp barrier
 	// 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];
+	    if (p_outarray[ti] != NULL) {
+	      p_outarray[0]->insert(*(p_outarray[ti]));
+	      delete p_outarray[ti];
+	      p_outarray[ti] = NULL;
+	    }
 	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	    delete vtppoanarray[ti];
 	  }
@@ -1407,7 +1378,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
 #endif
 }
 	
-int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp) {
+int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp) {
+  const dcomplex cc0 = 0.0 + I * 0.0;
   int nxi = sconf->number_of_scales;
   char virtual_line[256];
   string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
@@ -1434,12 +1406,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
   int last_configuration;
   dcomplex ent, entn;
   double enti;
+  int num_configs = sconf->configurations;
+  int ndirs = sa->nkks;
+  int oindex = -1;
+  int jindex = jxi488 - output->first_xi + 1;
 
 #ifdef USE_NVTX
   nvtxRangePush("Prepare matrix calculation");
 #endif
-  sprintf(virtual_line, "========== JXI =%3d ====================\n", jxi488);
-  output->append_line(virtual_line);
   double xi = sconf->get_scale(jxi488 - 1);
   double exdc = sconf->exdc;
   double exri = sqrt(exdc);
@@ -1448,14 +1422,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
   if (idfc >= 0) {
     cid->vk = xi * cid->wn;
     vkarg = cid->vk;
-    sprintf(virtual_line, "  VK=%15.7lE, XI=%15.7lE\n", cid->vk, xi);
-    output->append_line(virtual_line);
+    output->vec_vk[jindex - 1] = cid->vk;
+    output->vec_xi[jindex - 1] = xi;
     // goes to 120
   } else { // label 119
     vkarg = xi * cid->vk;
     cid->sqsfi = 1.0 / (xi * xi);
-    sprintf(virtual_line, "  XI=%15.7lE\n", xi);
-    output->append_line(virtual_line);
+    output->vec_vk[jindex - 1] = cid->vk;
+    output->vec_xi[jindex - 1] = xi;
   }
   // label 120
   double sze = vkarg * cid->extr;
@@ -1494,8 +1468,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
       }
       indme(i133, npnt, npntts, vkarg, ent, enti, entn, jer, lcalc, cid->arg, cid->c1);
       if (jer != 0) {
-	sprintf(virtual_line, "  STOP IN INDME\n");
-	output->append_line(virtual_line);
+	output->vec_ier[jindex - 1] = 1;
+	output->vec_jxi[jindex - 1] = -jxi488;
 	message = "ERROR: indme failed with error code " + to_string(jer) + ".\n";
 	logger->log(message, LOG_ERRO);
 	delete logger;
@@ -1506,8 +1480,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
   } // i133 loop
   ospv(cid->c1, vkarg, sze, exri, entn, enti, jer, lcalc, cid->arg);
   if (jer != 0) {
-    sprintf(virtual_line, "  STOP IN OSPV\n");
-    output->append_line(virtual_line);
+    output->vec_ier[jindex - 1] = 2;
+    output->vec_jxi[jindex - 1] = -jxi488;
     message = "ERROR: ospv failed with error code " + to_string(jer) + ".\n";
     logger->log(message, LOG_ERRO);
     delete logger;
@@ -1610,19 +1584,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
     }
   }
   // label 156: continue from here
+  last_configuration = 0;
   for (int i168 = 1; i168 <= nsph; i168++) {
     if (cid->c1->iog[i168 - 1] >= i168) {
+      int i = i168 - 1;
+      last_configuration++;
+      oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
       if (cid->c1->nshl[i168 - 1] != 1) {
-	sprintf(virtual_line, "  SPHERE N.%2d: SIZE=%15.7lE\n", i168, cid->c1->vsz[i168 - 1]);
-	output->append_line(virtual_line);
+	output->vec_sphere_ref_indices[oindex] = cc0;
+	output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
       } else {
-	sprintf(virtual_line, "  SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", i168, cid->c1->vsz[i168 - 1], real(cid->c1->vkt[i168 - 1]), imag(cid->c1->vkt[i168 - 1]));
-	output->append_line(virtual_line);
+	output->vec_sphere_ref_indices[oindex] = cid->c1->vkt[i];
+	output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
       }
     } 
   } // i168 loop
-  sprintf(virtual_line, "  EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", sze, real(entn), imag(entn));
-  output->append_line(virtual_line);
+  last_configuration++;
+  oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
+  output->vec_sphere_sizes[oindex] = sze;
+  output->vec_sphere_ref_indices[oindex] = entn;
   // label 160
   double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
   double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;
@@ -1642,6 +1622,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
   nvtxRangePush("Angle loop");
 #endif
   double th = sa->th;
+  int done_dirs = 0;
   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;
@@ -1767,8 +1748,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 		vtppoanp->append_line(VirtualBinaryLine(value));
 	      }
 	    }
-	    sprintf(virtual_line, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);
-	    output->append_line(virtual_line);
 	    int jlr = 2;
 	    for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
 	      int ipol = (ilr210 % 2 == 0) ? 1 : -1;
@@ -1784,68 +1763,72 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 	      double qschum = imag(s0m) * csch;
 	      double pschum = real(s0m) * csch;
 	      double s0magm = cabs(s0m) * cs0;
-	      if (inpol == 0) {
-		sprintf(virtual_line, "   LIN %2d\n", ipol);
-		output->append_line(virtual_line);
-	      } else { // label 206
-		sprintf(virtual_line, "  CIRC %2d\n", ipol);
-		output->append_line(virtual_line);
-	      }
+	      // if (inpol == 0) {
+	      // sprintf(virtual_line, "   LIN %2d\n", ipol);
+	      // output->append_line(virtual_line);
+	      // } else { // label 206
+	      // sprintf(virtual_line, "  CIRC %2d\n", ipol);
+	      // output->append_line(virtual_line);
+	      // }
 	      // label 208
-	      sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
-	      output->append_line(virtual_line);
-	      double alamb = 2.0 * 3.141592653589793238 / cid->vk;
-	      sprintf(virtual_line, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n", ipol, alamb, scasm, abssm, extsm);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
-	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "  FSAS(%1d,%1d)=%15.7lE%15.7lE   FSAS(%1d,%1d)=%15.7lE%15.7lE\n",
-		      ilr210, ilr210, real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]),
-		      imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
-		      real(cid->c1->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1->fsacm[jlr - 1][ilr210 - 1])
-		      );
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
-	      output->append_line(virtual_line);
 	      double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
 	      double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1];
 	      double fz = rapr;
-	      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fk=%15.7lE\n", fz);
-	      output->append_line(virtual_line);
+	      if (ipol == -1) {
+		output->vec_scs1[jindex - 1] = scasm;
+		output->vec_abs1[jindex - 1] = abssm;
+		output->vec_exs1[jindex - 1] = extsm;
+		output->vec_albeds1[jindex - 1] = albdm;
+		output->vec_scsrt1[jindex - 1] = qscam;
+		output->vec_absrt1[jindex - 1] = qabsm;
+		output->vec_exsrt1[jindex - 1] = qextm;
+		output->vec_fsas11[jindex - 1] = cid->c1->fsacm[0][0];
+		output->vec_fsas21[jindex - 1] = cid->c1->fsacm[1][0];
+		output->vec_qschu1[jindex -1] = qschum;
+		output->vec_pschu1[jindex -1] = pschum;
+		output->vec_s0mag1[jindex -1] = s0magm;
+		output->vec_cosav1[jindex -1] = cosav;
+		output->vec_raprs1[jindex -1] = rapr;
+		output->vec_fk1[jindex -1] = fz;
+	      } else if (ipol == 1) {
+		output->vec_scs2[jindex - 1] = scasm;
+		output->vec_abs2[jindex - 1] = abssm;
+		output->vec_exs2[jindex - 1] = extsm;
+		output->vec_albeds2[jindex - 1] = albdm;
+		output->vec_scsrt2[jindex - 1] = qscam;
+		output->vec_absrt2[jindex - 1] = qabsm;
+		output->vec_exsrt2[jindex - 1] = qextm;
+		output->vec_fsas22[jindex - 1] = cid->c1->fsacm[1][1];
+		output->vec_fsas12[jindex - 1] = cid->c1->fsacm[0][1];
+		output->vec_qschu2[jindex -1] = qschum;
+		output->vec_pschu2[jindex -1] = pschum;
+		output->vec_s0mag2[jindex -1] = s0magm;
+		output->vec_cosav2[jindex -1] = cosav;
+		output->vec_raprs2[jindex -1] = rapr;
+		output->vec_fk2[jindex -1] = fz;
+	      }
 	    } // ilr210 loop
-	    double rmbrif = (real(cid->c1->fsacm[0][0]) - real(cid->c1->fsacm[1][1])) / real(cid->c1->fsacm[0][0]);
-	    double rmdchr = (imag(cid->c1->fsacm[0][0]) - imag(cid->c1->fsacm[1][1])) / imag(cid->c1->fsacm[0][0]);
-	    sprintf(virtual_line, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rmbrif);
-	    output->append_line(virtual_line);
-	    sprintf(virtual_line, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rmdchr);
-	    output->append_line(virtual_line);
 	  }
 	  // label 212
-	  sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", 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", th, ph, ths, phs);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  SCAND=%10.3lE\n", cid->scan);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  CFMP=%15.7lE, SFMP=%15.7lE\n", cid->cfmp, cid->sfmp);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  CFSP=%15.7lE, SFSP=%15.7lE\n", cid->cfsp, cid->sfsp);
-	  output->append_line(virtual_line);
+	  output->vec_dir_scand[done_dirs] = cid->scan;
+	  output->vec_dir_cfmp[done_dirs] = cid->cfmp;
+	  output->vec_dir_sfmp[done_dirs] = cid->sfmp;
+	  output->vec_dir_cfsp[done_dirs] = cid->cfsp;
+	  output->vec_dir_sfsp[done_dirs] = cid->sfsp;
 	  if (isam >= 0) {
-	    sprintf(virtual_line, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n", 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", cid->uns[0], cid->uns[1], cid->uns[2]);
-	    output->append_line(virtual_line);
+	    output->vec_dir_un[3 * done_dirs] = cid->un[0];
+	    output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
+	    output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
+	    output->vec_dir_uns[3 * done_dirs] = cid->uns[0];
+	    output->vec_dir_uns[3 * done_dirs + 1] = cid->uns[1];
+	    output->vec_dir_uns[3 * done_dirs + 2] = cid->uns[2];
 	  } else { // label 214
-	    sprintf(virtual_line, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", cid->un[0], cid->un[1], cid->un[2]);
-	    output->append_line(virtual_line);
+	    output->vec_dir_un[3 * done_dirs] = cid->un[0];
+	    output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
+	    output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
+	    output->vec_dir_uns[3 * done_dirs] = 0.0;
+	    output->vec_dir_uns[3 * done_dirs + 1] = 0.0;
+	    output->vec_dir_uns[3 * done_dirs + 2] = 0.0;
 	  }
 	  // label 220
 	  pcros(cid->vk, exri, cid->c1);
@@ -1926,8 +1909,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 	      vtppoanp->append_line(VirtualBinaryLine(value));
 	    }
 	  }
-	  sprintf(virtual_line, "     SINGLE SCATTERER\n");
-	  output->append_line(virtual_line);
+	  oindex = (jindex - 1) * ndirs + done_dirs;
 	  int jlr = 2;
 	  for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
 	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;
@@ -1943,53 +1925,38 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 	    double qschu = imag(s0) * csch;
 	    double pschu = real(s0) * csch;
 	    double s0mag = cabs(s0) * cs0;
-	    if (inpol == 0) {
-	      sprintf(virtual_line, "   LIN %2d\n", ipol);
-	      output->append_line(virtual_line);
-	    } else { // label 273
-	      sprintf(virtual_line, "  CIRC %2d\n", ipol);
-	      output->append_line(virtual_line);
-	    }
 	    // label 275
-	    sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
-		    scasec, abssec, extsec, albedc
-		    );
-	    output->append_line(virtual_line);
-	    double alam = 2.0 * 3.141592653589793238 / cid->vk;
-	    sprintf(virtual_line, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n", ipol, alam, scasec, abssec, extsec);
-	    sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, " %14.7lE%15.7lE%15.7lE\n",
-		    qsca, qabs, qext
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line,
-		    "  FSAS(%1d,%1d)=%15.7lE%15.7lE   FSAS(%1d,%1d)=%15.7lE%15.7lE\n",
-		    ilr290, ilr290, real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]),
-		    imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,
-		    real(cid->c1->fsac[jlr - 1][ilr290 - 1]),
-		    imag(cid->c1->fsac[jlr - 1][ilr290 - 1])
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line,
-		    "   SAS(%1d,%1d)=%15.7lE%15.7lE    SAS(%1d,%1d)=%15.7lE%15.7lE\n",
-		    ilr290, ilr290, real(cid->c1->sac[ilr290 - 1][ilr290 - 1]),
-		    imag(cid->c1->sac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,
-		    real(cid->c1->sac[jlr - 1][ilr290 - 1]),
-		    imag(cid->c1->sac[jlr - 1][ilr290 - 1])
-		    );
-	    output->append_line(virtual_line);
-	    sprintf(
-		    virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
-		    qschu, pschu, s0mag
-		    );
-	    output->append_line(virtual_line);
+	    if (ipol == -1) {
+	      output->vec_dir_scs1[oindex] = scasec;
+	      output->vec_dir_abs1[oindex] = abssec;
+	      output->vec_dir_exs1[oindex] = extsec;
+	      output->vec_dir_albeds1[oindex] = albedc;
+	      output->vec_dir_scsrt1[oindex] = qsca;
+	      output->vec_dir_absrt1[oindex] = qabs;
+	      output->vec_dir_exsrt1[oindex] = qext;
+	      output->vec_dir_fsas11[oindex] = cid->c1->fsac[0][0];
+	      output->vec_dir_fsas21[oindex] = cid->c1->fsac[1][0];
+	      output->vec_dir_sas11[oindex] = cid->c1->sac[0][0];
+	      output->vec_dir_sas21[oindex] = cid->c1->sac[1][0];
+	      output->vec_dir_qschu1[oindex] = qschu;
+	      output->vec_dir_pschu1[oindex] = pschu;
+	      output->vec_dir_s0mag1[oindex] = s0mag;
+	    } else if (ipol == 1) {
+	      output->vec_dir_scs2[oindex] = scasec;
+	      output->vec_dir_abs2[oindex] = abssec;
+	      output->vec_dir_exs2[oindex] = extsec;
+	      output->vec_dir_albeds2[oindex] = albedc;
+	      output->vec_dir_scsrt2[oindex] = qsca;
+	      output->vec_dir_absrt2[oindex] = qabs;
+	      output->vec_dir_exsrt2[oindex] = qext;
+	      output->vec_dir_fsas22[oindex] = cid->c1->fsac[1][1];
+	      output->vec_dir_fsas12[oindex] = cid->c1->fsac[0][1];
+	      output->vec_dir_sas22[oindex] = cid->c1->sac[1][1];
+	      output->vec_dir_sas12[oindex] = cid->c1->sac[0][1];
+	      output->vec_dir_qschu2[oindex] = qschu;
+	      output->vec_dir_pschu2[oindex] = pschu;
+	      output->vec_dir_s0mag2[oindex] = s0mag;
+	    }
 	    bool goto290 = isam >= 0 && (jths > 1 || jphs > 1);
 	    if (!goto290) {
 	      cid->gapv[0] = cid->gap[0][ilr290 - 1];
@@ -1999,12 +1966,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 	      double scatts = cid->c1->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", cosav, rapr);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz);
-	      output->append_line(virtual_line);
+	      if (ipol == -1) {
+		output->vec_dir_cosav1[oindex] = cosav;
+		output->vec_dir_rapr1[oindex] = rapr;
+		output->vec_dir_fl1[oindex] = fp;
+		output->vec_dir_fr1[oindex] = fn;
+		output->vec_dir_fk1[oindex] = fk;
+		output->vec_dir_fx1[oindex] = fx;
+		output->vec_dir_fy1[oindex] = fy;
+		output->vec_dir_fz1[oindex] = fz;
+	      } else if (ipol == 1) {
+		output->vec_dir_cosav2[oindex] = cosav;
+		output->vec_dir_rapr2[oindex] = rapr;
+		output->vec_dir_fl2[oindex] = fp;
+		output->vec_dir_fr2[oindex] = fn;
+		output->vec_dir_fk2[oindex] = fk;
+		output->vec_dir_fx2[oindex] = fx;
+		output->vec_dir_fy2[oindex] = fy;
+		output->vec_dir_fz2[oindex] = fz;
+	      }
 	      cid->tqev[0] = cid->tqce[ilr290 - 1][0];
 	      cid->tqev[1] = cid->tqce[ilr290 - 1][1];
 	      cid->tqev[2] = cid->tqce[ilr290 - 1][2];
@@ -2013,45 +1993,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 	      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", tep, ten, tek);
-	      output->append_line(virtual_line);
-	      sprintf(virtual_line, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n", tsp, tsn, tsk);
-	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n",
-		      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",
-		      cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2]
-		      );
-	      output->append_line(virtual_line);
+	      if (ipol == -1) {
+		output->vec_dir_tqel1[oindex] = tep;
+		output->vec_dir_tqer1[oindex] = ten;
+		output->vec_dir_tqek1[oindex] = tek;
+		output->vec_dir_tqsl1[oindex] = tsp;
+		output->vec_dir_tqsr1[oindex] = tsn;
+		output->vec_dir_tqsk1[oindex] = tsk;
+		output->vec_dir_tqex1[oindex] = cid->tqce[0][0];
+		output->vec_dir_tqey1[oindex] = cid->tqce[0][1];
+		output->vec_dir_tqez1[oindex] = cid->tqce[0][2];
+		output->vec_dir_tqsx1[oindex] = cid->tqcs[0][0];
+		output->vec_dir_tqsy1[oindex] = cid->tqcs[0][1];
+		output->vec_dir_tqsz1[oindex] = cid->tqcs[0][2];
+	      } else if (ipol == 1) {
+		output->vec_dir_tqel2[oindex] = tep;
+		output->vec_dir_tqer2[oindex] = ten;
+		output->vec_dir_tqek2[oindex] = tek;
+		output->vec_dir_tqsl2[oindex] = tsp;
+		output->vec_dir_tqsr2[oindex] = tsn;
+		output->vec_dir_tqsk2[oindex] = tsk;
+		output->vec_dir_tqex2[oindex] = cid->tqce[1][0];
+		output->vec_dir_tqey2[oindex] = cid->tqce[1][1];
+		output->vec_dir_tqez2[oindex] = cid->tqce[1][2];
+		output->vec_dir_tqsx2[oindex] = cid->tqcs[1][0];
+		output->vec_dir_tqsy2[oindex] = cid->tqcs[1][1];
+		output->vec_dir_tqsz2[oindex] = cid->tqcs[1][2];
+	      }
 	    }
 	  } //ilr290 loop
-	  double rbirif = (real(cid->c1->fsac[0][0]) - real(cid->c1->fsac[1][1])) / real(cid->c1->fsac[0][0]);
-	  double rdichr = (imag(cid->c1->fsac[0][0]) - imag(cid->c1->fsac[1][1])) / imag(cid->c1->fsac[0][0]);
-	  sprintf(virtual_line, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rbirif);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rdichr);
-	  output->append_line(virtual_line);
-	  sprintf(virtual_line, "  MULL\n");
-	  output->append_line(virtual_line);
 	  for (int i = 0; i < 4; i++) {
-	    sprintf(
-		    virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
-		    cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
-		    );
-	    output->append_line(virtual_line);
+	    oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
+	    output->vec_dir_mull[oindex] = cid->cmul[i][0];
+	    output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
+	    output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
+	    output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
 	  }
-	  sprintf(virtual_line, "  MULLLR\n");
-	  output->append_line(virtual_line);
 	  for (int i = 0; i < 4; i++) {
-	    sprintf(
-		    virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
-		    cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
-		    );
-	    output->append_line(virtual_line);
+	    oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
+	    output->vec_dir_mulllr[oindex] = cid->cmullr[i][0];
+	    output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
+	    output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
+	    output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
 	  }
 	  if (iavm != 0) {
 	    mmulc(cid->c1->vintm, cid->cmullr, cid->cmul);
@@ -2069,37 +2052,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
 		vtppoanp->append_line(VirtualBinaryLine(value));
 	      }
 	    }
-	    sprintf(virtual_line, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);
-	    output->append_line(virtual_line);
-	    if (inpol == 0) {
-	      sprintf(virtual_line, "   LIN\n");
-	      output->append_line(virtual_line);
-	    } else { // label 316
-	      sprintf(virtual_line, "  CIRC\n");
-	      output->append_line(virtual_line);
-	    }
 	    // label 318
-	    sprintf(virtual_line, "  MULL\n");
-	    output->append_line(virtual_line);
 	    for (int i = 0; i < 4; i++) {
-	      sprintf(
-		      virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
-		      cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
-		      );
-	      output->append_line(virtual_line);
+	      oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
+	      output->vec_dir_mull[oindex] = cid->cmul[i][0];
+	      output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
+	      output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
+	      output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
 	    }
-	    sprintf(virtual_line, "  MULLLR\n");
-	    output->append_line(virtual_line);
 	    for (int i = 0; i < 4; i++) {
-	      sprintf(
-		      virtual_line, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
-		      cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
-		      );
-	      output->append_line(virtual_line);
+	      oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
+	      output->vec_dir_mulllr[oindex] = cid->cmul[i][0];
+	      output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
+	      output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
+	      output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
 	    }
 	  }
 	  // label 420, continues jphs loop
 	  if (isam < 1) phs += sa->phsstp;
+	  done_dirs++;
 	} // jphs loop, labeled 480
 	if (isam <= 1) thsl += sa->thsstp;
       } // jths loop, labeled 482