From 4c68446ea1dc67ccaa73027e577ade9e3acd0be6 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 6 Jun 2024 11:22:26 +0200
Subject: [PATCH] Protect MPI-dependent closing brace with MPI pre-compiler
 directive

---
 src/cluster/cluster.cpp | 1428 ++++++++++++++++++++-------------------
 1 file changed, 715 insertions(+), 713 deletions(-)

diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 67dda9a6..a6d640b8 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -482,7 +482,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
 #endif
 #pragma omp barrier
 
-	} // close strided loop running on MPI processes
+	} // close strided loop running on MPI processes, ixi488 loop
 	// delete the shared arrays I used to make available to thread 0 the virtual files of other threads
 #pragma omp barrier
 	if (myompthread == 0) {
@@ -633,769 +633,818 @@ void cluster(const string& config_file, const string& data_file, const string& o
     fclose(timing_file);
     delete time_logger;
     delete logger;
+#ifdef MPI_VERSION
   }
+#endif
 }
 
 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 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);
+  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);
-    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");
+    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 {
-      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
+      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);
       } else {
-	int nsh = cid->c1->nshl[i132 - 1];
-	int ici = (nsh + 1) / 2;
-	if (idfc == 0) {
+	if (jxi488 == 1) {
 	  for (int ic = 0; ic < ici; ic++)
-	    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;
+	    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;
       }
-    } // i132 loop
+    }
+    if (jer != 0) {
+      return jer;
+      //break;
+    }
+  } // 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);
+  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 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");
+  }
+  // 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);
-    }
-    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);
+      if (cid->c1->nshl[i] != 1) {
+	sprintf(virtual_line, "  SIZE=%15.7lE\n\0", cid->c2->vsz[i]);
 	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]);
+      } 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);
       }
-    } // 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);
+      // 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
+  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
 		  );
-	    // 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;
+	    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];
+	      }
+	    }
 	  }
-	} else { // label 180, NK == 1 AND JXI488 == 1
-	  if (isam >= 0) {
-	    // label 182
+	  // 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;
 	  }
-	}
-	// 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];
-		}
+	  // 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 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 < 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 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 = cid->c1ao->ecscm[i];
+	    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->c1ao->ecscpm[i]);
+		value = real(cid->gappm[i][j]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->c1ao->ecscpm[i]);
+		value = imag(cid->gappm[i][j]);
 		// 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->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));
-		}
-	      }
-	      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);
-		}
-		// 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);
+	    }
+	    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);
-		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);
+	      } else { // label 206
+		sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
 		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);
+	      }
+	      // label 208
+	      sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0");
 	      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);
+	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", scasm, abssm, extsm, albdm);
 	      output->append_line(virtual_line);
-	    }
-	    // label 212
-	    sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n\0", jth486, jph484, jths, jphs);
+	      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);
-	    sprintf(virtual_line, "  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n\0", th, ph, ths, phs);
+	    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);
-	    sprintf(virtual_line, "  SCAND=%10.3lE\n\0", cid->scan);
+	  }
+	  // 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, "  CFMP=%15.7lE, SFMP=%15.7lE\n\0", cid->cfmp, cid->sfmp);
+	    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, "  CFSP=%15.7lE, SFSP=%15.7lE\n\0", cid->cfsp, cid->sfsp);
+	  } 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);
-	    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]);
-	      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]);
+	  }
+	  // label 220
+	  if (inpol == 0) {
+	    sprintf(virtual_line, "   LIN\n\0");
+	    output->append_line(virtual_line);
+	  } else { // label 222
+	    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) {
+	    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);
-	    }
-	    // label 220
-	    if (inpol == 0) {
-	      sprintf(virtual_line, "   LIN\n\0");
+	      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])
+		      );
 	      output->append_line(virtual_line);
-	    } else { // label 222
-	      sprintf(virtual_line, "  CIRC\n\0");
+	      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);
-	    }
-	    // 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");
+	      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 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);
+	      for (int i1 = 0; i1 < 4; i1++) {
 		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])
+			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, "  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])
+			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);
-		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
-	      }
-	    } // 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));
-		}
+	      } // 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 i = 0; i < 2; i++) {
-		double value = cid->c1ao->scsc[i];
+	    }
+	    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];
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = real(cid->c1ao->scscp[i]);
+		value = real(cid->gapp[i][j]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->c1ao->scscp[i]);
+		value = imag(cid->gapp[i][j]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = cid->c1ao->ecsc[i];
+	      }
+	    }
+	    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->c1ao->ecscp[i]);
+		value = real(cid->tqcpe[i][j]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = imag(cid->c1ao->ecscp[i]);
+		value = imag(cid->tqcpe[i][j]);
 		// 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];
-		  // 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];
+	    }
+	    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 = cid->up[i];
+		value = real(cid->tqcps[i][j]);
 		// tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 		vtppoanp->append_line(VirtualBinaryLine(value));
-		value = cid->un[i];
+		value = imag(cid->tqcps[i][j]);
 		// 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]);
+	    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];
 	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
 	      vtppoanp->append_line(VirtualBinaryLine(value));
-	      value = imag(cid->c1->vint[i]);
+	      value = cid->un[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));
-	      }
+	  }
+	  // 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));
 	    }
-	    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);
-	      sprintf(
-		      virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0",
-		      scasec, abssec, extsec, albedc
-		      );
+	  }
+	  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);
-	      sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0");
+	    } else { // label 273
+	      sprintf(virtual_line, "  CIRC %2d\n\0", ipol);
 	      output->append_line(virtual_line);
-	      sprintf(
-		      virtual_line, " %14.7lE%15.7lE%15.7lE\n\0",
-		      qsca, qabs, qext
-		      );
+	    }
+	    // 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);
 	      output->append_line(virtual_line);
-	      sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0");
+	      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, " %14.7lE%15.7lE%15.7lE\n\0",
-		      scarat, absrat, extrat
-		      );
+	      sprintf(virtual_line, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n\0", fx, fy, fz);
 	      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])
-		      );
+	      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, "   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])
-		      );
+	      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, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0",
-		      ilr290, ilr290, refinr, ilr290, ilr290, extcor
+		      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, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0",
-		      qschu, pschu, s0mag
+		      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);
-	      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);
-		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);
+	    }
+	  } //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, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rdichr);
+	  }
+	  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");
+	      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++) {
@@ -1414,74 +1463,27 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
 		      );
 	      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");
-		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
+	  }
+	  // 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;
+}
-- 
GitLab