diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 06410a2b71fd03a834806a6421f2bfb39ad97025..a00425f17549acabe349860a19c4877928416a4f 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -21,6 +21,7 @@
 #include <cstdio>
 #include <exception>
 #include <fstream>
+#include <hdf5.h>
 #include <string>
 #ifdef _OPENMP
 #include <omp.h>
@@ -73,11 +74,19 @@
 #include "../include/algebraic.h"
 #endif
 
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
 using namespace std;
 
 // I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
 
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, FILE *output, const string& output_path, fstream& tppoan);
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, VirtualAsciiFile *output, const string& output_path, fstream& tppoan);
 
 /*! \brief C++ implementation of CLU
  *
@@ -152,49 +161,71 @@ void cluster(const string& config_file, const string& data_file, const string& o
       // 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");
+      //FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
+      VirtualAsciiFile *p_output = new VirtualAsciiFile();
+      char virtual_line[256];
       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");
-      fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
-      fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
+      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",
 	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
 	      gconf->npntts, gconf->iavm, gconf->iavm
 	      );
-      fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
-      for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
-						gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
-						);
-      fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
-      fprintf(
-	      output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n\0");
+      p_output->append_line(virtual_line);
+      for (int ri = 0; ri < nsph; ri++) {
+	sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n\0",
+		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\0");
+      p_output->append_line(virtual_line);
+      sprintf(
+	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n\0",
 	      p_scattering_angles->th, p_scattering_angles->thstp,
 	      p_scattering_angles->thlst, p_scattering_angles->ths,
 	      p_scattering_angles->thsstp, p_scattering_angles->thslst
 	      );
-      fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
-      fprintf(
-	      output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(
+	      virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n\0",
 	      p_scattering_angles->ph, p_scattering_angles->phstp,
 	      p_scattering_angles->phlst, p_scattering_angles->phs,
 	      p_scattering_angles->phsstp, p_scattering_angles->phslst
 	      );
-      fprintf(output, " READ(IR,*)JWTM\n");
-      fprintf(output, " %5d\n", gconf->jwtm);
-      fprintf(output, "  READ(ITIN)NSPHT\n");
-      fprintf(output, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");
-      fprintf(output, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
-      fprintf(output, "  READ(ITIN)(XIV(I),I=1,NXI)\n");
-      fprintf(output, "  READ(ITIN)NSHL(I),ROS(I)\n");
-      fprintf(output, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
-      fprintf(output, " \n");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, " READ(IR,*)JWTM\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, " %5d\n\0", gconf->jwtm);
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)NSPHT\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)(IOG(I),I=1,NSPH)\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)(XIV(I),I=1,NXI)\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)NSHL(I),ROS(I)\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n\0");
+      p_output->append_line(virtual_line);
+      sprintf(virtual_line, " \n\0");
+      p_output->append_line(virtual_line);
       str(sconf, cid->c1, cid->c1ao, cid->c3, cid->c4, cid->c6);
       thdps(cid->c4->lm, cid->zpv);
       double exdc = sconf->exdc;
       double exri = sqrt(exdc);
-      fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+      sprintf(virtual_line, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n\0", exri);
+      p_output->append_line(virtual_line);
       fstream *tppoanp = new fstream;
       fstream &tppoan = *tppoanp;
       string tppoan_name = output_path + "/c_TPPOAN";
@@ -225,8 +256,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
 	if (sconf->idfc < 0) {
 	  cid->vk = cid->xip * cid->wn;
-	  fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);
-	  fprintf(output, " \n");
+	  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;
@@ -234,7 +267,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 #ifdef USE_NVTX
 	nvtxRangePush("First iteration");
 #endif
-	int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, output, output_path, tppoan);
+	int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, tppoan);
 #ifdef USE_NVTX
 	nvtxRangePop();
 #endif
@@ -251,7 +284,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	  // First loop failed. Halt the calculation.
 	  tppoan.close();
 	  fclose(timing_file);
-	  fclose(output);
+	  //fclose(output);
+	  delete p_output;
 	  delete p_scattering_angles;
 	  delete cid;
 	  delete logger;
@@ -287,17 +321,20 @@ void cluster(const string& config_file, const string& data_file, const string& o
 #endif
 	  // 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;
+	  //FILE *output_2 = NULL;
+	  VirtualAsciiFile *p_output_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;
+	    //output_2 = output;
+	    p_output_2 = p_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");
+	    //output_2 = fopen((output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), "w");
+	    p_output_2 = new VirtualAsciiFile();
 	    tppoanp_2 = new fstream;
 	    tppoanp_2->open((output_path + "/c_TPPOAN_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
 	  }
@@ -310,14 +347,16 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	  // 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, output_2, output_path, *tppoanp_2);
+	    int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, *tppoanp_2);
 	  }
 
 #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);
+	    //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;
 	  }
@@ -341,9 +380,16 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	    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");
+	    char virtual_line[256];
+	    int index = 0;
 	    int c = fgetc(partial_output);
 	    while (c != EOF) {
-	      fputc(c, output);
+	      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);
@@ -384,7 +430,18 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	      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);
+		// 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);
 	      }
@@ -413,7 +470,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
       } else { // In case TPPOAN could not be opened. Should never happen.
 	logger->err("\nERROR: failed to open TPPOAN file.\n");
       }
-      fclose(output);
+      // fclose(output);
+      p_output->write_to_disk(output_path + "/c_OCLU");
+      delete p_output;
       // Clean memory
       delete cid;
       delete p_scattering_angles;
@@ -460,7 +519,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
 #endif
       // 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;
+      //FILE *output_2 = NULL;
+      VirtualAsciiFile *p_output_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) {
@@ -471,7 +531,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	// 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");
+      // output_2 = fopen((output_path + "/c_OCLU_" + to_string(mpidata->rank) + "_" + to_string(myompthread)).c_str(), "w");
+      p_output_2 = new VirtualAsciiFile();
       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;
@@ -481,7 +542,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
       // 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, output_2, output_path, *tppoanp_2);
+	int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, *tppoanp_2);
       }
 
 #pragma omp barrier
@@ -489,7 +550,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
       if (myompthread != 0) {
 	delete cid_2;
       }
-      fclose(output_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;
 #pragma omp barrier
@@ -587,9 +650,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
   delete logger;
 }
 
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, FILE *output, const string& output_path, fstream& tppoan)
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, VirtualAsciiFile *output, const string& output_path, fstream& tppoan)
 {
   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);
@@ -617,7 +681,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 #ifdef USE_NVTX
   nvtxRangePush("Prepare matrix calculation");
 #endif
-  fprintf(output, "========== JXI =%3d ====================\n", jxi488);
+  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);
@@ -626,15 +691,18 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   if (idfc >= 0) {
     cid->vk = xi * cid->wn;
     vkarg = cid->vk;
-    fprintf(output, "  VK=%15.7lE, XI=%15.7lE\n", cid->vk, xi);
+    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);
-    fprintf(output, "  XI=%15.7lE\n", 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) {
-    fprintf(output, "  STOP IN HJV\n");
+    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
   }
@@ -663,7 +731,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	  cid->c1, cid->c2, jer, lcalc, cid->arg
 	  );
       if (jer != 0) {
-	fprintf(output, "  STOP IN DME\n");
+	sprintf(virtual_line, "  STOP IN DME\n\0");
+	output->append_line(virtual_line);
 	return jer;
 	//break;
       }
@@ -722,9 +791,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   }
   // label 156: continue from here
   if (inpol == 0) {
-    fprintf(output, "   LIN\n");
+    sprintf(virtual_line, "   LIN\n\0");
+    output->append_line(virtual_line);
   } else { // label 158
-    fprintf(output, "  CIRC\n");
+    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);
@@ -734,7 +805,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   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) fprintf(output, "     SPHERES; LMX=LI\n");
+  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;
@@ -742,38 +816,52 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
       cid->c1->sqscs[i] *= cid->sqsfi;
       cid->c1->sqabs[i] *= cid->sqsfi;
       cid->c1->sqexs[i] *= cid->sqsfi;
-      fprintf(output, "     SPHERE %2d\n", i170);
+      sprintf(virtual_line, "     SPHERE %2d\n\0", i170);
+      output->append_line(virtual_line);
       if (cid->c1->nshl[i] != 1) {
-	fprintf(output, "  SIZE=%15.7lE\n", cid->c2->vsz[i]);
+	sprintf(virtual_line, "  SIZE=%15.7lE\n\0", cid->c2->vsz[i]);
+	output->append_line(virtual_line);
       } else { // label 162
-	fprintf(output, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i]));
+	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
-      fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
-      fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i], albeds);
-      fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
-      fprintf(output, " %14.7lE%15.7lE%15.7lE\n", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]);
-      fprintf(output, "  FSAS=%15.7lE%15.7lE\n", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i]));
+      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);
       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;
-      fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
+      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];
-      fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
-      fprintf(output, "  IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", cid->tqse[0][i], cid->tqss[0][i]);
-      fprintf(output, "  IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", cid->tqse[1][i], cid->tqss[1][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
-  fprintf(output, "  FSAT=%15.7lE%15.7lE\n", real(cid->c3->tfsas), imag(cid->c3->tfsas));
+  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;
-  fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
+  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));
   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);
@@ -782,7 +870,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 #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";
+  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
@@ -914,7 +1002,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	      }
 	    }
-	    fprintf(output, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
+	    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;
@@ -938,104 +1027,142 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	      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) {
-		fprintf(output, "   LIN %2d\n", ipol);
+		sprintf(virtual_line, "   LIN %2d\n\0", ipol);
+		output->append_line(virtual_line);
 	      } else { // label 206
-		fprintf(output, "  CIRC %2d\n", ipol);
+		sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
+		output->append_line(virtual_line);
 	      }
 	      // label 208
-	      fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
-	      fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
-	      fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
-	      fprintf(output, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
-	      fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
-	      fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
-	      fprintf(
-		      output, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
+	      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])
 		      );
-	      fprintf(
-		      output, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
+	      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
 		      );
-	      fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
+	      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;
-	      fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
-	      fprintf(output, "  Fk=%15.7lE\n", fz);
+	      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);
 	    } // 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]);
-	    fprintf(output, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
-	    fprintf(output, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
+	    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
-	  fprintf(output, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, jths, jphs);
-	  fprintf(output, "  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n", th, ph, ths, phs);
-	  fprintf(output, "  SCAND=%10.3lE\n", cid->scan);
-	  fprintf(output, "  CFMP=%15.7lE, SFMP=%15.7lE\n", cid->cfmp, cid->sfmp);
-	  fprintf(output, "  CFSP=%15.7lE, SFSP=%15.7lE\n", cid->cfsp, cid->sfsp);
+	  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) {
-	    fprintf(output, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n", cid->un[0], cid->un[1], cid->un[2]);
-	    fprintf(output, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n", cid->uns[0], cid->uns[1], cid->uns[2]);
+	    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]);
+	    output->append_line(virtual_line);
 	  } else { // label 214
-	    fprintf(output, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", cid->un[0], cid->un[1], cid->un[2]);
+	    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) {
-	    fprintf(output, "   LIN\n");
+	    sprintf(virtual_line, "   LIN\n\0");
+	    output->append_line(virtual_line);
 	  } else { // label 222
-	    fprintf(output, "  CIRC\n");
+	    sprintf(virtual_line, "  CIRC\n\0");
+	    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) fprintf(output, "     SPHERES; LMX=MIN0(LI,LE)\n");
+	  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) {
-	      fprintf(output, "     SPHERE %2d\n", i226);
-	      fprintf(
-		      output, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
+	      sprintf(virtual_line, "     SPHERE %2d\n\0", i226);
+	      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])
 		      );
-	      fprintf(
-		      output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
+	      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])
 		      );
+	      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);
-	      fprintf(output, "  MULS\n");
+	      sprintf(virtual_line, "  MULS\n\0");
+	      output->append_line(virtual_line);
 	      for (int i1 = 0; i1 < 4; i1++) {
-		fprintf(
-			output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		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
-	      fprintf(output, "  MULSLR\n");
+	      sprintf(virtual_line, "  MULSLR\n\0");
+	      output->append_line(virtual_line);
 	      for (int i1 = 0; i1 < 4; i1++) {
-		fprintf(
-			output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+		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
 	    }
 	  } // i226 loop
-	  fprintf(
-		  output, "  SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
+	  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])
 		  );
-	  fprintf(
-		  output, "  SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
+	  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])
 		  );
-	  fprintf(output, "     CLUSTER\n");
+	  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);
@@ -1137,44 +1264,56 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	    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) {
-	      fprintf(output, "   LIN %2d\n", ipol);
+	      sprintf(virtual_line, "   LIN %2d\n\0", ipol);
+	      output->append_line(virtual_line);
 	    } else { // label 273
-	      fprintf(output, "  CIRC %2d\n", ipol);
+	      sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
+	      output->append_line(virtual_line);
 	    }
 	    // label 275
-	    fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
-	    fprintf(
-		    output, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+	    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
 		    );
-	    fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
-	    fprintf(
-		    output, " %14.7lE%15.7lE%15.7lE\n",
+	    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
 		    );
-	    fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
-	    fprintf(
-		    output, " %14.7lE%15.7lE%15.7lE\n",
+	    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
 		    );
-	    fprintf(
-		    output, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
+	    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])
 		    );
-	    fprintf(
-		    output, "   SAC(%1d,%1d)=%15.7lE%15.7lE    SAC(%1d,%1d)=%15.7lE%15.7lE\n",
+	    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])
 		    );
-	    fprintf(
-		    output, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
+	    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
 		    );
-	    fprintf(
-		    output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+	    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);
 	    bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
 	    if (!goto190) {
 	      cid->gapv[0] = cid->gap[0][ilr290 - 1];
@@ -1184,9 +1323,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	      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);
-	      fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
-	      fprintf(output, "  Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
-	      fprintf(output, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", 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];
@@ -1195,35 +1337,45 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 	      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);
-	      fprintf(output, "   TQEl=%15.7lE,  TQEr=%15.7lE,  TQEk=%15.7lE\n", tep, ten, tek);
-	      fprintf(output, "   TQSl=%15.7lE,  TQSr=%15.7lE,  TQSk=%15.7lE\n", tsp, tsn, tsk);
-	      fprintf(
-		      output, "   TQEx=%15.7lE,  TQEy=%15.7lE,  TQEz=%15.7lE\n",
+	      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]
 		      );
-	      fprintf(
-		      output, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n",
+	      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]);
-	  fprintf(output, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
-	  fprintf(output, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
-	  fprintf(output, "  MULC\n");
+	  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++) {
-	    fprintf(
-		    output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+	    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);
 	  }
-	  fprintf(output, "  MULCLR\n");
+	  sprintf(virtual_line, "  MULCLR\n\0");
+	  output->append_line(virtual_line);
 	  for (int i = 0; i < 4; i++) {
-	    fprintf(
-		    output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+	    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);
@@ -1241,26 +1393,33 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	      }
 	    }
-	    fprintf(output, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
+	    sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm);
+	    output->append_line(virtual_line);
 	    if (inpol == 0) {
-	      fprintf(output, "   LIN\n");
+	      sprintf(virtual_line, "   LIN\n\0");
+	      output->append_line(virtual_line);
 	    } else { // label 316
-	      fprintf(output, "  CIRC\n");
+	      sprintf(virtual_line, "  CIRC\n\0");
+	      output->append_line(virtual_line);
 	    }
 	    // label 318
-	    fprintf(output, "  MULC\n");
+	    sprintf(virtual_line, "  MULC\n\0");
+	    output->append_line(virtual_line);
 	    for (int i = 0; i < 4; i++) {
-	      fprintf(
-		      output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+	      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);
 	    }
-	    fprintf(output, "  MULCLR\n");
+	    sprintf(virtual_line, "  MULCLR\n\0");
+	    output->append_line(virtual_line);
 	    for (int i = 0; i < 4; i++) {
-	      fprintf(
-		      output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
+	      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