/*! \file cluster.cpp
 */
#include <cstdio>
#include <complex>
#include <exception>
#include <fstream>
#include <string>

#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif

#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif

#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif

#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif

using namespace std;

/*! \brief C++ implementation of CLU
 *
 *  \param config_file: `string` Name of the configuration file.
 *  \param data_file: `string` Name of the input data file.
 *  \param output_path: `string` Directory to write the output files in.
 */
void cluster(string config_file, string data_file, string output_path) {
  printf("INFO: making legacy configuration...");
  ScattererConfiguration *sconf = NULL;
  try {
    sconf = ScattererConfiguration::from_dedfb(config_file);
  } catch(const OpenConfigurationFileException &ex) {
    printf("\nERROR: failed to open scatterer configuration file.\n");
    printf("FILE: %s\n", ex.what());
    exit(1);
  }
  sconf->write_formatted(output_path + "/c_OEDFB");
  sconf->write_binary(output_path + "/c_TEDF");
  sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");
  GeometryConfiguration *gconf = NULL;
  try {
    gconf = GeometryConfiguration::from_legacy(data_file);
  } catch (const OpenConfigurationFileException &ex) {
    printf("\nERROR: failed to open geometry configuration file.\n");
    printf("FILE: %s\n", ex.what());
    if (sconf) delete sconf;
    exit(1);
  }
  printf(" done.\n");
  if (sconf->number_of_spheres == gconf->number_of_spheres) {
    // Shortcuts to variables stored in configuration objects
    int nsph = gconf->number_of_spheres;
    int mxndm = gconf->mxndm;
    int inpol = gconf->in_pol;
    int npnt = gconf->npnt;
    int npntts = gconf->npntts;
    int iavm = gconf->iavm;
    int isam = gconf->meridional_type;
    int nxi = sconf->number_of_scales;
    double th = gconf->in_theta_start;
    double thstp = gconf->in_theta_step;
    double thlst = gconf->in_theta_end;
    double ths = gconf->sc_theta_start;
    double thsstp = gconf->sc_theta_step;
    double thslst = gconf->sc_theta_end;
    double ph = gconf->in_phi_start;
    double phstp = gconf->in_phi_step;
    double phlst = gconf->in_phi_end;
    double phs = gconf->sc_phi_start;
    double phsstp = gconf->sc_phi_step;
    double phslst = gconf->sc_phi_end;
    double wp = sconf->wp;
    // Global variables for CLU
    int lm = gconf->l_max;
    if (gconf->li > lm) lm = gconf->li;
    if (gconf->le > lm) lm = gconf->le;
    C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
    C3 *c3 = new C3();
    for (int c1i = 0; c1i < nsph; c1i++) {
      c1->rxx[c1i] = gconf->sph_x[c1i];
      c1->ryy[c1i] = gconf->sph_y[c1i];
      c1->rzz[c1i] = gconf->sph_z[c1i];
      c1->ros[c1i] = sconf->radii_of_spheres[c1i];
    }
    C4 *c4 = new C4;
    c4->li = gconf->li;
    c4->le = gconf->le;
    c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;
    c4->nv3j = (c4->lm * (c4->lm  + 1) * (2 * c4->lm + 7)) / 6;
    c4->nsph = nsph;
    // The following is needed to initialize C1_AddOns
    c4->litpo = c4->li + c4->li + 1;
    c4->litpos = c4->litpo * c4->litpo;
    c4->lmtpo = c4->li + c4->le + 1;
    c4->lmtpos = c4->lmtpo * c4->lmtpo;
    c4->nlim = c4->li * (c4->li + 2);
    c4->nlem = c4->le * (c4->le + 2);
    c4->lmpo = c4->lm + 1;
    C1_AddOns *c1ao = new C1_AddOns(c4);
    // End of add-ons initialization
    C6 *c6 = new C6(c4->lmtpo);
    FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
    int jer = 0, lcalc = 0;
    complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);
    int max_ici = 0;
    for (int insh = 0; insh < nsph; insh++) {
      int nsh = sconf->nshl_vec[insh];
      int ici = (nsh + 1) / 2;
      if (ici > max_ici) max_ici = ici;
    }
    C2 *c2 = new C2(nsph, max_ici, npnt, npntts);
    complex<double> **am = new complex<double>*[mxndm];
    for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm]();
    const int ndi = c4->nsph * c4->nlim;
    C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
    double *gaps = new double[nsph]();
    double *tqev = new double[3]();
    double *tqsv = new double[3]();
    double **tqse, **tqss, **tqce, **tqcs;
    complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
    tqse = new double*[2];
    tqspe = new complex<double>*[2];
    tqss = new double*[2];
    tqsps = new complex<double>*[2];
    tqce = new double*[2];
    tqcpe = new complex<double>*[2];
    tqcs = new double*[2];
    tqcps = new complex<double>*[2];
    for (int ti = 0; ti < 2; ti++) {
      tqse[ti] = new double[nsph]();
      tqspe[ti] = new complex<double>[nsph]();
      tqss[ti] = new double[nsph]();
      tqsps[ti] = new complex<double>[nsph]();
      tqce[ti] = new double[3]();
      tqcpe[ti] = new complex<double>[3]();
      tqcs[ti] = new double[3]();
      tqcps[ti] = new complex<double>[3]();
    }
    double *gapv = new double[3]();
    complex<double> **gapp, **gappm;
    double **gap, **gapm;
    gapp = new complex<double>*[3];
    gappm = new complex<double>*[3];
    gap = new double*[3];
    gapm = new double*[3];
    for (int gi = 0; gi < 3; gi++) {
      gapp[gi] = new complex<double>[2]();
      gappm[gi] = new complex<double>[2]();
      gap[gi] = new double[2]();
      gapm[gi] = new double[2]();
    }
    double *u = new double[3]();
    double *us = new double[3]();
    double *un = new double[3]();
    double *uns = new double[3]();
    double *up = new double[3]();
    double *ups = new double[3]();
    double *unmp = new double[3]();
    double *unsmp = new double[3]();
    double *upmp = new double[3]();
    double *upsmp = new double[3]();
    double *argi = new double[1]();
    double *args = new double[1]();
    double *duk = new double[3]();
    double **cextlr, **cext;
    double **cmullr, **cmul;
    cextlr = new double*[4];
    cext = new double*[4];
    cmullr = new double*[4];;
    cmul = new double*[4];
    for (int ci = 0; ci < 4; ci++) {
      cextlr[ci] = new double[4]();
      cext[ci] = new double[4]();
      cmullr[ci] = new double[4]();
      cmul[ci] = new double[4]();
    }
    int isq, ibf;
    double scan, cfmp, sfmp, cfsp, sfsp;
    // End of global variables for CLU
    fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
    fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
	    nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
	    );
    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->sph_x[ri], gconf->sph_y[ri], gconf->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",
	    th, thstp, thlst, ths, thsstp, 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",
	    ph, phstp, phlst, phs, phsstp, 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");
    double small = 1.0e-3;
    int nth = 0, nph = 0;
    if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
    nth++;
    if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
    nph++;
    int nths = 0, nphs = 0;
    double thsca = 0.0;
    if (isam > 1) {
      nths = 1;
      thsca = ths - th;
    } else { // ISAM <= 1
      if (thsstp == 0.0) nths = 0;
      else nths = int ((thslst - ths) / thsstp + small);
      nths++;
    }
    if (isam >= 1) {
      nphs = 1;
    } else {
      if (phsstp == 0.0) nphs = 0;
      else nphs = int((phslst - phs) / phsstp + small);
      nphs++;
    }
    int nk = nth * nph;
    int nks = nths * nphs;
    int nkks = nk * nks;
    double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
    str(sconf->rcf, c1, c1ao, c3, c4, c6);
    double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
    for (int zi = 0; zi < c4->lm; zi++) {
      zpv[zi] = new double**[3];
      for (int zj = 0; zj < 3; zj++) {
	zpv[zi][zj] = new double*[2];
	for (int zk = 0; zk < 2; zk++) {
	  zpv[zi][zj][zk] = new double[2]();
	}
      }
    }
    thdps(c4->lm, zpv);
    double exri = sqrt(sconf->exdc);
    double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
    fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
    fstream tppoan;
    string tppoan_name = output_path + "/c_TPPOAN";
    tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
    if (tppoan.is_open()) {
      tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&nxi), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
      double wn = wp / 3.0e8;
      double sqsfi = 1.0;
      if (sconf->idfc < 0) {
	vk = sconf->xip * wn;
	fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
	fprintf(output, " \n");
      }
      for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
	printf("INFO: running scale iteration %d...", jxi488);
	int jaw = 1;
	fprintf(output, "========== JXI =%3d ====================\n", jxi488);
	double xi = sconf->scale_vec[jxi488 - 1];
	double vkarg = 0.0;
	if (sconf->idfc >= 0) {
	  vk = xi * wn;
	  vkarg = vk;
	  fprintf(output, "  VK=%15.7lE, XI=%15.7lE\n", vk, xi);
	} else {
	  vkarg = xi * vk;
	  sqsfi = 1.0 / (xi * xi);
	  fprintf(output, "  XI=%15.7lE\n", xi);
	}
	hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
	if (jer != 0) {
	  fprintf(output, "  STOP IN HJV\n");
	  break; // jxi488 loop: goes to memory cleaning and  return
	}
	for (int i132 = 1; i132 <= nsph; i132++) {
	  int iogi = c1->iog[i132 - 1];
	  if (iogi != i132) {
	    for (int l123 = 1; l123 <= gconf->li; l123++) {
	      c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
	      c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
	    } // l123 loop
	  } else {
	    int nsh = sconf->nshl_vec[i132 - 1];
	    int ici = (nsh + 1) / 2;
	    if (sconf->idfc == 0) {
	      for (int ic = 0; ic < ici; ic++)
		c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
	    } else {
	      if (jxi488 == 1) {
		for (int ic = 0; ic < ici; ic++)
		  c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
	      }
	    }
	    if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
	    dme(
		c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,
		c1, c2, jer, lcalc, arg
		);
	    if (jer != 0) {
	      fprintf(output, "  STOP IN DME\n");
	      break;
	    }
	  }
	  if (jer != 0) break;
	} // i132 loop
	cms(am, c1, c1ao, c4, c6);
	int ndit = 2 * nsph * c4->nlim;
	lucin(am, mxndm, ndit, jer);
	if (jer != 0) break; // jxi488 loop: goes to memory clean
	ztm(am, c1, c1ao, c4, c6, c9);
	if (sconf->idfc >= 0) {
	  if (jxi488 == gconf->jwtm) {
	    int is = 1;
	    fstream ttms_file;
	    string ttms_name = output_path + "/c_TTMS";
	    ttms_file.open(ttms_name.c_str(), ios::out | ios::binary);
	    if (ttms_file.is_open()) {
	      ttms_file.write(reinterpret_cast<char *>(&is), sizeof(int));
	      ttms_file.write(reinterpret_cast<char *>(&lm), sizeof(int));
	      ttms_file.write(reinterpret_cast<char *>(&vk), sizeof(double));
	      ttms_file.write(reinterpret_cast<char *>(&exri), sizeof(double));
	      int nlemt = 2 * c4->nlem;
	      for (int ami = 0; ami < nlemt; ami++) {
		for (int amj = 0; amj < nlemt; amj++) {
		  complex<double> value = c1ao->am0m[ami][amj];
		  double rval = value.real();
		  double ival = value.imag();
		  ttms_file.write(reinterpret_cast<char *>(&rval), sizeof(double));
		  ttms_file.write(reinterpret_cast<char *>(&ival), sizeof(double));
		}
	      }
	      ttms_file.close();
	    } else { // Could not open TM file. Should never occur.
	      printf("\nERROR: failed to open TTMS file.\n");
	      break;
	    }
	  }
	}
	// label 156: continue from here
	if (inpol == 0) {
	  fprintf(output, "   LIN\n");
	} else { // label 158
	  fprintf(output, "  CIRC\n");
	}
	// label 160
	double cs0 = 0.25 * vk * vk * vk / acos(0.0);
	double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
	std::complex<double> s0(0.0, 0.0);
	scr0(vk, exri, c1, c1ao, c3, c4);
	//printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
	double sqk = vk * vk * sconf->exdc;
	aps(zpv, c4->li, nsph, c1, sqk, gaps);
	rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
	if (c4->li != c4->le) fprintf(output, "     SPHERES; LMX=LI\n");
	for (int i170 = 1; i170 <= nsph; i170++) {
	  if (c1->iog[i170 - 1] >= i170) {
	    int i = i170 - 1;
	    double albeds = c1->sscs[i] / c1->sexs[i];
	    c1->sqscs[i] *= sqsfi;
	    c1->sqabs[i] *= sqsfi;
	    c1->sqexs[i] *= sqsfi;
	    fprintf(output, "     SPHERE %2d\n", i170);
	    if (c1->nshl[i] != 1) {
	      fprintf(output, "  SIZE=%15.7lE\n", c2->vsz[i]);
	    } else { // label 162
	      fprintf(output, "  SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
	    }
	    // label 164
	    fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
	    fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
	    fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
	    fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
	    fprintf(output, "  FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
	    csch = 2.0 * vk * sqsfi / c1->gcsv[i];
	    s0 = c1->fsas[i] * exri;
	    qschu = s0.imag() * csch;
	    pschu = s0.real() * csch;
	    s0mag = abs(s0) * cs0;
	    fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
	    double rapr = c1->sexs[i] - gaps[i];
	    double cosav = gaps[i] / 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", tqse[0][i], tqss[0][i]);
	    fprintf(output, "  IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
	  }
	} // i170 loop
	fprintf(output, "  FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
	csch = 2.0 * vk * sqsfi / c3->gcs;
	s0 = c3->tfsas * exri;
	qschu = s0.imag() * csch;
	pschu = s0.real() * csch;
	s0mag = abs(s0) * cs0;
	fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
	tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
	pcrsm0(vk, exri, inpol, c1, c1ao, c4);
	apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
	th = th1;
	for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
	  ph = ph1;
	  double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
	  for (int jph484 = 1; jph484 <= nph; jph484++) {
	    int jw = 0;
	    if (nk != 1 || jxi488 <= 1) {
	      upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
	      if (isam >= 0) {
		wmamp(
		      0, cost, sint, cosp, sinp, inpol, c4->le, 0,
		      nsph, argi, u, upmp, unmp, c1
		      );
		// label 182
		apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
		raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
		jw = 1;
	      }
	    } else { // label 180, NK == 1 AND JXI488 == 1
	      if (isam >= 0) {
		// label 182
		apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
		raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
		jw = 1;
	      }
	    }
	    // label 184
	    double thsl = ths1;
	    double phsph = 0.0;
	    for (int jths = 1; jths <= nths; jths++) {
	      ths = thsl;
	      int icspnv = 0;
	      if (isam > 1) ths += 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
	      phs = phs1;
	      for (int jphs = 1; jphs <= nphs; jphs++) {
		double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
		if (isam >= 1) {
		  phs = ph + phsph;
		  if (phs > 360.0) phs -= 360.0;
		}
		// label 188
		bool goto190 = (nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
		if (!goto190) {
		  upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
		  if (isam >= 0)
		    wmamp(
			  2, costs, sints, cosps, sinps, inpol, c4->le,
			  0, nsph, args, us, upsmp, unsmp, c1
			  );
		}
		// label 190
		if (nkks != 1 || jxi488 <= 1) {
		  upvsp(
			u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
			duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
			);
		  if (isam < 0) {
		    wmasp(
			  cost, sint, cosp, sinp, costs, sints, cosps, sinps,
			  u, up, un, us, ups, uns, isq, ibf, inpol, c4->le,
			  0, nsph, argi, args, c1
			  );
		  } else { // label 192
		    for (int i193 = 0; i193 < 3; i193++) {
		      up[i193] = upmp[i193];
		      un[i193] = unmp[i193];
		      ups[i193] = upsmp[i193];
		      uns[i193] = unsmp[i193];
		    }
		  }
		}
		// label 194
		if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
		if (isam < 0) {
		  apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
		  raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
		  jw = 1;
		}
		// label 196
		tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&scan), sizeof(double));
		if (jaw != 0) {
		  jaw = 0;
		  mextc(vk, exri, c1ao->fsacm, cextlr, 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 = cext[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  for (int i = 0; i < 3; i++) {
		    double value = c1ao->scscm[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->scscpm[i].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->scscpm[i].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecscm[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecscpm[i].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecscpm[i].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		  for (int i = 0; i < 3; i++) {
		    for (int j = 0; j < 2; j++) {
		      double value = gapm[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = gappm[i][j].real();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = gappm[i][j].imag();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  fprintf(output, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
		  int jlr = 2;
		  for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
		    int ipol = (ilr210 % 2 == 0) ? 1 : -1;
		    if (ilr210 == 2) jlr = 1;
		    double extsm = c1ao->ecscm[ilr210 - 1];
		    double qextm = extsm * sqsfi / c3->gcs;
		    double extrm = extsm / c3->ecs;
		    double scasm = c1ao->scscm[ilr210 - 1];
		    double albdm = scasm / extsm;
		    double qscam = scasm *sqsfi / c3->gcs;
		    double scarm = scasm / c3->scs;
		    double abssm = extsm - scasm;
		    double qabsm = abssm * sqsfi / c3->gcs;
		    double absrm = abssm / c3->acs;
		    double acsecs = c3->acs / c3->ecs;
		    if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
		    complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
		    double qschum = s0m.imag() * csch;
		    double pschum = s0m.real() * csch;
		    double s0magm = abs(s0m) * cs0;
		    double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real();
		    double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag();
		    if (inpol == 0) {
		      fprintf(output, "   LIN %2d\n", ipol);
		    } else { // label 206
		      fprintf(output, "  CIRC %2d\n", ipol);
		    }
		    // 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",
			    ilr210, ilr210, c1ao->fsacm[ilr210 - 1][ilr210 - 1].real(),
			    c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,
			    c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()
			    );
		    fprintf(
			    output, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
			    ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
			    );
		    fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
		    double rapr = c1ao->ecscm[ilr210 - 1] - gapm[2][ilr210 - 1];
		    double cosav = gapm[2][ilr210 - 1] / 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);
		  } // ilr210 loop
		  double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();
		  double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();
		  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);
		}
		// 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", scan);
		fprintf(output, "  CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
		fprintf(output, "  CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
		if (isam >= 0) {
		  fprintf(output, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
		  fprintf(output, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
		} else { // label 214
		  fprintf(output, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
		}
		// label 220
		if (inpol == 0) {
		  fprintf(output, "   LIN\n");
		} else { // label 222
		  fprintf(output, "  CIRC\n");
		}
		// label 224
		scr2(vk, vkarg, exri, duk, c1, c1ao, c3, c4);
		if (c4->li != c4->le) fprintf(output, "     SPHERES; LMX=MIN0(LI,LE)\n");
		for (int i226 = 1; i226 <= nsph; i226++) {
		  if (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",
			    c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
			    c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
			    );
		    fprintf(
			    output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
			    c1->sas[i226 - 1][0][1].real(), c1->sas[i226 - 1][0][1].imag(),
			    c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()
			    );
		    for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
		      c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
		    } // j225 loop
		    mmulc(c1ao->vint, cmullr, cmul);
		    fprintf(output, "  MULS\n");
		    for (int i1 = 0; i1 < 4; i1++) {
		      fprintf(
			      output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			      cmul[i1][0], cmul[i1][1], cmul[i1][2], cmul[i1][3]
			      );
		    } // i1 loop
		    fprintf(output, "  MULSLR\n");
		    for (int i1 = 0; i1 < 4; i1++) {
		      fprintf(
			      output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			      cmullr[i1][0], cmullr[i1][1], cmullr[i1][2], cmullr[i1][3]
			      );
		    } // i1 loop
		  }
		} // i226 loop
		fprintf(
			output, "  SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
			c3->tsas[0][0].real(), c3->tsas[0][0].imag(),
			c3->tsas[1][0].real(), c3->tsas[1][0].imag()
			);
		fprintf(
			output, "  SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
			c3->tsas[0][1].real(), c3->tsas[0][1].imag(),
			c3->tsas[1][1].real(), c3->tsas[1][1].imag()
			);
		fprintf(output, "     CLUSTER\n");
		pcros(vk, exri, c1, c1ao, c4);
		mextc(vk, exri, c1ao->fsac, cextlr, cext);
		mmulc(c1ao->vint, cmullr, 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 = cext[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  for (int i = 0; i < 2; i++) {
		    double value = c1ao->scsc[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->scscp[i].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->scscp[i].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecsc[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecscp[i].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->ecscp[i].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		  for (int i = 0; i < 3; i++) {
		    for (int j = 0; j < 2; j++) {
		      double value = gap[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = gapp[i][j].real();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = gapp[i][j].imag();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  for (int i = 0; i < 2; i++) {
		    for (int j = 0; j < 3; j++) {
		      double value = tqce[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = tqcpe[i][j].real();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = tqcpe[i][j].imag();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  for (int i = 0; i < 2; i++) {
		    for (int j = 0; j < 3; j++) {
		      double value = tqcs[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = tqcps[i][j].real();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		      value = tqcps[i][j].imag();
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  for (int i = 0; i < 3; i++) {
		    double value = u[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = up[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = un[i];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		}
		// label 254
		for (int i = 0; i < 16; i++) {
		  double value = c1ao->vint[i].real();
		  tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  value = c1ao->vint[i].imag();
		  tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		for (int i = 0; i < 4; i++) {
		  for (int j = 0; j < 4; j++) {
		    double value = cmul[i][j];
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		}
		int jlr = 2;
		for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
		  int ipol = (ilr290 % 2 == 0) ? 1 : -1;
		  if (ilr290 == 2) jlr = 1;
		  double extsec = c1ao->ecsc[ilr290 - 1];
		  double qext = extsec * sqsfi / c3->gcs;
		  double extrat = extsec / c3->ecs;
		  double scasec = c1ao->scsc[ilr290 - 1];
		  double albedc = scasec / extsec;
		  double qsca = scasec * sqsfi / c3->gcs;
		  double scarat = scasec / c3->scs;
		  double abssec = extsec - scasec;
		  double qabs = abssec * sqsfi / c3->gcs;
		  double absrat = 1.0;
		  double ratio = c3->acs / c3->ecs;
		  if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
		  s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
		  double qschu = s0.imag() * csch;
		  double pschu = s0.real() * csch;
		  s0mag = abs(s0) * cs0;
		  double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
		  double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
		  if (inpol == 0) {
		    fprintf(output, "   LIN %2d\n", ipol);
		  } else { // label 273
		    fprintf(output, "  CIRC %2d\n", ipol);
		  }
		  // label 275
		  fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
		  fprintf(
			  output, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
			  scasec, abssec, extsec, albedc
			  );
		  fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
		  fprintf(
			  output, " %14.7lE%15.7lE%15.7lE\n",
			  qsca, qabs, qext
			  );
		  fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
		  fprintf(
			  output, " %14.7lE%15.7lE%15.7lE\n",
			  scarat, absrat, extrat
			  );
		  fprintf(
			  output, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
			  ilr290, ilr290, c1ao->fsac[ilr290 - 1][ilr290 - 1].real(), c1ao->fsac[ilr290 - 1][ilr290 - 1].imag(),
			  jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()
			  );
		  fprintf(
			  output, "   SAC(%1d,%1d)=%15.7lE%15.7lE    SAC(%1d,%1d)=%15.7lE%15.7lE\n",
			  ilr290, ilr290, c1ao->sac[ilr290 - 1][ilr290 - 1].real(), c1ao->sac[ilr290 - 1][ilr290 - 1].imag(),
			  jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()
			  );
		  fprintf(
			  output, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
			  ilr290, ilr290, refinr, ilr290, ilr290, extcor
			  );
		  fprintf(
			  output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
			  qschu, pschu, s0mag
			  );
		  bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
		  if (!goto190) {
		    gapv[0] = gap[0][ilr290 - 1];
		    gapv[1] = gap[1][ilr290 - 1];
		    gapv[2] = gap[2][ilr290 - 1];
		    double extins = c1ao->ecsc[ilr290 - 1];
		    double scatts = c1ao->scsc[ilr290 - 1];
		    double rapr, cosav, fp, fn, fk, fx, fy, fz;
		    rftr(u, up, un, 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);
		    tqev[0] = tqce[ilr290 - 1][0];
		    tqev[1] = tqce[ilr290 - 1][1];
		    tqev[2] = tqce[ilr290 - 1][2];
		    tqsv[0] = tqcs[ilr290 - 1][0];
		    tqsv[1] = tqcs[ilr290 - 1][1];
		    tqsv[2] = tqcs[ilr290 - 1][2];
		    double tep, ten, tek, tsp, tsn, tsk;
		    tqr(u, up, un, tqev, 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",
			    tqce[ilr290 - 1][0], tqce[ilr290 - 1][1], tqce[ilr290 - 1][2]
			    );
		    fprintf(
			    output, "   TQSx=%15.7lE,  TQSy=%15.7lE,  TQSz=%15.7lE\n",
			    tqcs[ilr290 - 1][0], tqcs[ilr290 - 1][1], tqcs[ilr290 - 1][2]
			    );
		  }
		} //ilr290 loop
		double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();
		double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();
		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");
		for (int i = 0; i < 4; i++) {
		  fprintf(
			  output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			  cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
			  );
		}
		fprintf(output, "  MULCLR\n");
		for (int i = 0; i < 4; i++) {
		  fprintf(
			  output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			  cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
			  );
		}
		if (iavm != 0) {
		  mmulc(c1ao->vintm, cmullr, cmul);
		  // Some implicit loops writing to binary.
		  for (int i = 0; i < 16; i++) {
		    double value;
		    value = c1ao->vintm[i].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = c1ao->vintm[i].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		  for (int i = 0; i < 4; i++) {
		    for (int j = 0; j < 4; j++) {
		      double value = cmul[i][j];
		      tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    }
		  }
		  fprintf(output, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
		  if (inpol == 0) {
		    fprintf(output, "   LIN\n");
		  } else { // label 316
		    fprintf(output, "  CIRC\n");
		  }
		  // label 318
		  fprintf(output, "  MULC\n");
		  for (int i = 0; i < 4; i++) {
		    fprintf(
			    output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			    cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
			    );
		  }
		  fprintf(output, "  MULCLR\n");
		  for (int i = 0; i < 4; i++) {
		    fprintf(
			    output, "        %15.7lE%15.7lE%15.7lE%15.7lE\n",
			    cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
			    );
		  }
		}
		// label 420, continues jphs loop
		if (isam < 1) phs += phsstp;
	      } // jphs loop, labeled 480
	      if (isam <= 1) thsl += thsstp;
	    } // jths loop, labeled 482
	    ph += phstp;
	  } // jph484 loop
	  th += thstp;
	} // jth486 loop
	printf(" done.\n");
      } // jxi488 loop
      tppoan.close();
    } else { // In case TPPOAN could not be opened. Should never happen.
      printf("\nERROR: failed to open TPPOAN file.\n");
    }
    fclose(output);
    // Clean memory
    delete c1;
    delete c1ao;
    delete c3;
    delete c4;
    delete c6;
    delete c9;
    for (int zi = c4->lm - 1; zi > -1; zi--) {
      for (int zj = 2; zj > -1; zj--) {
	delete[] zpv[zi][zj][1];
	delete[] zpv[zi][zj][0];
	delete[] zpv[zi][zj];
      }
      delete[] zpv[zi];
    }
    delete[] zpv;
    for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
    delete[] am;
    delete[] gaps;
    for (int ti = 1; ti > -1; ti--) {
      delete[] tqse[ti];
      delete[] tqss[ti];
      delete[] tqspe[ti];
      delete[] tqsps[ti];
      delete[] tqce[ti];
      delete[] tqcpe[ti];
      delete[] tqcs[ti];
      delete[] tqcps[ti];
    }
    delete[] tqse;
    delete[] tqss;
    delete[] tqspe;
    delete[] tqsps;
    delete[] tqce;
    delete[] tqcpe;
    delete[] tqcs;
    delete[] tqcps;
    delete[] tqev;
    delete[] tqsv;
    for (int gi = 2; gi > -1; gi--) {
      delete[] gapp[gi];
      delete[] gappm[gi];
      delete[] gap[gi];
      delete[] gapm[gi];
    }
    delete[] gapp;
    delete[] gappm;
    delete[] gap;
    delete[] gapm;
    delete[] gapv;
    delete[] u;
    delete[] us;
    delete[] un;
    delete[] uns;
    delete[] up;
    delete[] ups;
    delete[] unmp;
    delete[] unsmp;
    delete[] upmp;
    delete[] upsmp;
    delete[] argi;
    delete[] args;
    delete[] duk;
    for (int ci = 3; ci > -1; ci--) {
      delete[] cextlr[ci];
      delete[] cext[ci];
      delete[] cmullr[ci];
      delete[] cmul[ci];
    }
    delete[] cextlr;
    delete[] cext;
    delete[] cmullr;
    delete[] cmul;
  } else { // NSPH mismatch between geometry and scatterer configurations.
    throw UnrecognizedConfigurationException(
					     "Inconsistent geometry and scatterer configurations."
					     );
  }
  delete sconf;
  delete gconf;
  printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
}
