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

#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

using namespace std;

/*! \brief C++ implementation of SPH
 *
 *  \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 sphere(string config_file, string data_file, string output_path) {
  complex<double> arg, s0, tfsas;
  double th, ph;
  printf("INFO: making legacy configuration...\n");
  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");
  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);
  }
  if (sconf->number_of_spheres == gconf->number_of_spheres) {
    int isq, ibf;
    double *argi, *args, *gaps;
    double cost, sint, cosp, sinp;
    double costs, sints, cosps, sinps;
    double scan;
    double *duk = new double[3];
    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 *upmp = new double[3];
    double *upsmp = new double[3];
    double *unmp = new double[3];
    double *unsmp = new double[3];
    double **cmul = new double*[4];
    double **cmullr = new double*[4];
    for (int i = 0; i < 4; i++) {
      cmul[i] = new double[4];
      cmullr[i] = new double[4];
    }
    complex<double> **tqspe, **tqsps;
    double **tqse, **tqss;
    tqse = new double*[2];
    tqss = new double*[2];
    tqspe = new std::complex<double>*[2];
    tqsps = new std::complex<double>*[2];
    for (int ti = 0; ti < 2; ti++) {
      tqse[ti] = new double[2]();
      tqss[ti] = new double[2]();
      tqspe[ti] = new std::complex<double>[2]();
      tqsps[ti] = new std::complex<double>[2]();
    }
    double frx = 0.0, fry = 0.0, frz = 0.0;
    double cfmp, cfsp, sfmp, sfsp;
    complex<double> *vint = new complex<double>[16];
    int jw;
    int nsph = gconf->number_of_spheres;
    C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
    for (int i = 0; i < nsph; i++) {
      c1->ros[i] = sconf->radii_of_spheres[i];
    }
    C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts);
    argi = new double[1];
    args = new double[1];
    gaps = new double[2];
    FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w");
    fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
    fprintf(
	    output,
	    " %5d%5d%5d%5d%5d%5d\n",
	    gconf->number_of_spheres,
	    gconf->l_max,
	    gconf->in_pol,
	    gconf->npnt,
	    gconf->npntts,
	    gconf->meridional_type
	    );
    fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
    fprintf(
	    output,
	    "  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
	    gconf->in_theta_start,
	    gconf->in_theta_step,
	    gconf->in_theta_end,
	    gconf->sc_theta_start,
	    gconf->sc_theta_step,
	    gconf->sc_theta_end
	    );
    fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
    fprintf(
	    output,
	    "  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
	    gconf->in_phi_start,
	    gconf->in_phi_step,
	    gconf->in_phi_end,
	    gconf->sc_phi_start,
	    gconf->sc_phi_step,
	    gconf->sc_phi_end
	    );
    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 \n");
    double sml = 1.0e-3;
    int nth = 0, nph = 0;
    if (gconf->in_theta_step != 0.0)
      nth = int((gconf->in_theta_end - gconf->in_theta_start) / gconf->in_theta_step + sml);
    nth += 1;
    if (gconf->in_phi_step != 0.0)
      nph = int((gconf->in_phi_end - gconf->in_phi_start) / gconf->in_phi_step + sml);
    nph += 1;
    int nths = 0, nphs = 0;
    double thsca;
    if (gconf->meridional_type > 1) { // ISAM > 1, fixed scattering angle
      nths = 1;
      thsca = gconf->sc_theta_start - gconf->in_theta_start;
    } else { //ISAM <= 1
      if (gconf->in_theta_step != 0.0)
	nths = int((gconf->sc_theta_end - gconf->sc_theta_start) / gconf->sc_theta_step + sml);
      nths += 1;
      if (gconf->meridional_type == 1) { // ISAM = 1
	nphs = 1;
      } else { // ISAM < 1
	if (gconf->sc_phi_step != 0.0)
	  nphs = int((gconf->sc_phi_end - gconf->sc_phi_start) / gconf->sc_phi_step + sml);
	nphs += 1;
      }
    }
    int nk = nth * nph;
    int nks = nths * nphs;
    int nkks = nk * nks;
    double th1 = gconf->in_theta_start;
    double ph1 = gconf->in_phi_start;
    double ths1 = gconf->sc_theta_start;
    double phs1 = gconf->sc_phi_start;
    const double half_pi = acos(0.0);
    const double pi = 2.0 * half_pi;
    double gcs = 0.0;
    for (int i116 = 0; i116 < nsph; i116++) {
      int i = i116 + 1;
      int iogi = c1->iog[i116];
      if (iogi >= i) {
	double gcss = pi * c1->ros[i116] * c1->ros[i116];
	c1->gcsv[i116] = gcss;
	int nsh = c1->nshl[i116];
	for (int j115 = 0; j115 < nsh; j115++) {
	  c1->rc[i116][j115] = sconf->rcf[i116][j115] * c1->ros[i116];
	}
      }
      gcs += c1->gcsv[iogi];
    }
    double ****zpv = new double***[gconf->l_max]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
    for (int zi = 0; zi < gconf->l_max; 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(gconf->l_max, zpv);
    double exri = sqrt(sconf->exdc);
    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::binary|ios::out);
    if (tppoan.is_open()) {
      int imode = 10;
      tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&(gconf->meridional_type)), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&(gconf->in_pol)), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&(sconf->number_of_scales)), 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));

      for (int nsi = 0; nsi < nsph; nsi++)
	tppoan.write(reinterpret_cast<char *>(&(sconf->iog_vec[nsi])), sizeof(int));
      if (gconf->in_pol == 0) fprintf(output, "   LIN\n");
      else fprintf(output, "  CIRC\n");
      fprintf(output, " \n");
      double wn = sconf->wp / 3.0e8;
      double sqsfi = 1.0;
      double vk, vkarg;
      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 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
	printf("INFO: running scale iteration...\n");
	int jxi = jxi488 + 1;
	fprintf(output, "========== JXI =%3d ====================\n", jxi);
	double xi = sconf->scale_vec[jxi488];
	if (sconf->idfc >= 0) {
	  vk = xi * wn;
	  vkarg = vk;
	  fprintf(output, "  VK=%15.7lE, XI=%15.7lE\n", xi, vk);
	} else { // IDFC < 0
	  vkarg = xi * vk;
	  sqsfi = 1.0 / (xi * xi);
	  fprintf(output, "  XI=%15.7lE\n", xi);
	}
	tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
	for (int i132 = 0; i132 < nsph; i132++) {
	  int i = i132 + 1;
	  int iogi = sconf->iog_vec[i132];
	  if (iogi != i) {
	    for (int l123 = 0; l123 < gconf->l_max; l123++) {
	      c1->rmi[l123][i132] = c1->rmi[l123][iogi - 1];
	      c1->rei[l123][i132] = c1->rei[l123][iogi - 1];
	    }
	    continue; // i132
	  }
	  int nsh = sconf->nshl_vec[i132];
	  int ici = (nsh + 1) / 2;
	  if (sconf->idfc == 0) {
	    for (int ic = 0; ic < ici; ic++)
	      c2->dc0[ic] = sconf->dc0_matrix[ic][i132][0]; // WARNING: IDFC=0 is not tested!
	  } else { // IDFC != 0
	    if (jxi == 1) {
	      for (int ic = 0; ic < ici; ic++) {
		c2->dc0[ic] = sconf->dc0_matrix[ic][i132][jxi488];
	      }
	    }
	  }
	  if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
	  int jer = 0;
	  int lcalc = 0;
	  dme(
	      gconf->l_max, i, gconf->npnt, gconf->npntts, vkarg,
	      sconf->exdc, exri, c1, c2, jer, lcalc, arg
	      );
	  if (jer != 0) {
	    fprintf(output, "  STOP IN DME\n");
	    fprintf(output, "  AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, arg.real(), arg.imag());
	    tppoan.close();
	    fclose(output);
	    return;
	  }
	} // i132
	if (sconf->idfc >= 0 and nsph == 1 and jxi == gconf->jwtm) {
	  // This is the condition that writes the transition matrix to output.
	  int is = 1111;
	  fstream ttms;
	  string ttms_name = output_path + "/c_TTMS";
	  ttms.open(ttms_name.c_str(), ios::binary | ios::out);
	  if (ttms.is_open()) {
	    ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
	    ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
	    ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
	    ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
	    for (int lmi = 0; lmi < gconf->l_max; lmi++) {
	      complex<double> element1 = -1.0 / c1->rmi[lmi][0];
	      complex<double> element2 = -1.0 / c1->rei[lmi][0];
	      double vreal, vimag;
	      vreal = element1.real();
	      vimag = element1.imag();
	      ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
	      ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
	      vreal = element2.real();
	      vimag = element2.imag();
	      ttms.write(reinterpret_cast<char *>(&vreal), sizeof(double));
	      ttms.write(reinterpret_cast<char *>(&vimag), sizeof(double));
	    }
	    ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
	    ttms.close();
	  } else { // Failed to open output file. Should never happen.
	    printf("ERROR: could not open TTMS file.\n");
	    tppoan.close();
	  }
	}
	double cs0 = 0.25 * vk * vk * vk / half_pi;
	//printf("DEBUG: cs0 = %lE\n", cs0);
	sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
	//printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
	double sqk = vk * vk * sconf->exdc;
	aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
	rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
	for (int i170 = 0; i170 < nsph; i170++) {
	  int i = i170 + 1;
	  if (c1->iog[i170] >= i) {
	    double albeds = c1->sscs[i170] / c1->sexs[i170];
	    c1->sqscs[i170] *= sqsfi;
	    c1->sqabs[i170] *= sqsfi;
	    c1->sqexs[i170] *= sqsfi;
	    fprintf(output, "     SPHERE %2d\n", i);
	    if (c1->nshl[i170] != 1) {
	      fprintf(output, "  SIZE=%15.7lE\n", c2->vsz[i170]);
	    } else {
	      fprintf(
		      output,
		      "  SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
		      c2->vsz[i170],
		      c2->vkt[i170].real(),
		      c2->vkt[i170].imag()
		      );
	    }
	    fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
	    fprintf(
		    output,
		    " %14.7lE%15.7lE%15.7lE%15.7lE\n",
		    c1->sscs[i170], c1->sabs[i170],
		    c1->sexs[i170], albeds
		    );
	    fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
	    fprintf(
		    output,
		    " %14.7lE%15.7lE%15.7lE\n",
		    c1->sqscs[i170], c1->sqabs[i170],
		    c1->sqexs[i170]
		    );
	    fprintf(output, "  FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
	    double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
	    s0 = c1->fsas[i170] * exri;
	    double qschu = csch * s0.imag();
	    double pschu = csch * s0.real();
	    double s0mag = cs0 * abs(s0);
	    fprintf(
		    output,
		    "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
		    qschu, pschu, s0mag
		    );
	    double rapr = c1->sexs[i170] - gaps[i170];
	    double cosav = gaps[i170] / c1->sscs[i170];
	    fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
	    fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170], tqss[0][i170]);
	    fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
	    tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
	    tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
	    double val = tqspe[0][i170].real();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqspe[0][i170].imag();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqsps[0][i170].real();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqsps[0][i170].imag();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
	    tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
	    val = tqspe[1][i170].real();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqspe[1][i170].imag();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqsps[1][i170].real();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = tqsps[1][i170].imag();
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	  } // End if iog[i170] >= i
	} // i170 loop
	if (nsph != 1) {
	  fprintf(output, "  FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
	  double csch = 2.0 * vk * sqsfi / gcs;
	  s0 = tfsas * exri;
	  double qschu = csch * s0.imag();
	  double pschu = csch * s0.real();
	  double s0mag = cs0 * abs(s0);
	  fprintf(
		  output,
		  "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
		  qschu, pschu, s0mag
		  );
	}
	th = th1;
	for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section
	  int jth = jth486 + 1;
	  ph = ph1;
	  for (int jph484 = 0; jph484 < nph; jph484++) {
	    int jph = jph484 + 1;
	    bool goto182 = (nk == 1) && (jxi > 1);
	    if (!goto182) {
	      upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
	    }
	    if (gconf->meridional_type >= 0) {
	      wmamp(0, cost, sint, cosp, sinp, gconf->in_pol, gconf->l_max, 0, nsph, argi, u, upmp, unmp, c1);
	      for (int i183 = 0; i183 < nsph; i183++) {
		double rapr = c1->sexs[i183] - gaps[i183];
		frx = rapr * u[0];
		fry = rapr * u[1];
		frz = rapr * u[2];
	      }
	      jw = 1;
	    }
	    double thsl = ths1;
	    double phsph = 0.0;
	    for (int jths482 = 0; jths482 < nths; jths482++) {
	      int jths = jths482 + 1;
	      double ths = thsl;
	      int icspnv = 0;
	      if (gconf->meridional_type > 1) ths = th + thsca;
	      if (gconf->meridional_type >= 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;
	      }
	      double phs = phs1;
	      for (int jphs480 = 0; jphs480 < nphs; jphs480++) {
		int jphs = jphs480 + 1;
		if (gconf->meridional_type >= 1) {
		  phs = ph + phsph;
		  if (phs >= 360.0) phs -= 360.0;
		}
		bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));
		if (!goto190) {
		  upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
		  if (gconf->meridional_type >= 0) {
		    wmamp(2, costs, sints, cosps, sinps, gconf->in_pol, gconf->l_max, 0, nsph, args, us, upsmp, unsmp, c1);
		  }
		}
		if (nkks != 0 || jxi == 1) {
		  upvsp(u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns, duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp);
		  if (gconf->meridional_type < 0) {
		    wmasp(
			  cost, sint, cosp, sinp, costs, sints, cosps, sinps,
			  u, up, un, us, ups, uns, isq, ibf, gconf->in_pol,
			  gconf->l_max, 0, nsph, argi, args, c1
			  );
		  }
		  for (int i193 = 0; i193 < 3; i193++) {
		    un[i193] = unmp[i193];
		    uns[i193] = unsmp[i193];
		  }
		}
		if (gconf->meridional_type < 0) jw = 1;
		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 (jw != 0) {
		  jw = 0;
		  tppoan.write(reinterpret_cast<char *>(&(u[0])), sizeof(double));
		  tppoan.write(reinterpret_cast<char *>(&(u[1])), sizeof(double));
		  tppoan.write(reinterpret_cast<char *>(&(u[2])), sizeof(double));
		}
		fprintf(
			output,
			"********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
			jth, jph, 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 (gconf->meridional_type >= 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 {
		  fprintf(output, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
		}
		sscr2(nsph, gconf->l_max, vk, exri, c1);
		for (int ns226 = 0; ns226 < nsph; ns226++) {
		  int ns = ns226 + 1;
		  fprintf(output, "     SPHERE %2d\n", ns);
		  fprintf(
			  output, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
			  c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
			  c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
			  );
		  fprintf(
			  output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
			  c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
			  c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
			  );
		  if (jths == 1 && jphs == 1)
		    fprintf(
			    output, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
			    frx, fry, frz
			    );
		  for (int i225 = 0; i225 < 16; i225++) vint[i225] = c1->vints[ns226][i225];
		  mmulc(vint, cmullr, cmul);
		  fprintf(output, "  MULS\n        ");
		  for (int imul = 0; imul < 4; imul++) {
		    for (int jmul = 0; jmul < 4; jmul++) {
		      fprintf(output, "%15.7lE", cmul[imul][jmul]);
		    }
		    if (imul < 3) fprintf(output, "\n        ");
		    else fprintf(output, "\n");
		  }
		  fprintf(output, "  MULSLR\n        ");
		  for (int imul = 0; imul < 4; imul++) {
		    for (int jmul = 0; jmul < 4; jmul++) {
		      fprintf(output, "%15.7lE", cmullr[imul][jmul]);
		    }
		    if (imul < 3) fprintf(output, "\n        ");
		    else fprintf(output, "\n");
		  }
		  for (int vi = 0; vi < 16; vi++) {
		    double value = vint[vi].real();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = vint[vi].imag();
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		  }
		  for (int imul = 0; imul < 4; imul++) {
		    for (int jmul = 0; jmul < 4; jmul++) {
		      tppoan.write(reinterpret_cast<char *>(&(cmul[imul][jmul])), sizeof(double));
		    }
		  }
		} // ns226 loop
		if (gconf->meridional_type < 1) phs += gconf->sc_phi_step;
	      } // jphs480 loop
	      if (gconf->meridional_type <= 1) thsl += gconf->sc_theta_step;
	    } // jths482 loop
	    ph += gconf->in_phi_step;
	  } // jph484 loop on elevation
	  th += gconf->in_theta_step;
	} // jth486 loop on azimuth
	printf("INFO: done scale.\n");
      } //jxi488 loop on scales
      tppoan.close();
    } else { // In case TPPOAN could not be opened. Should never happen.
      printf("ERROR: failed to open TPPOAN file.\n");
    }
    fclose(output);
    delete c1;
    delete c2;
    for (int zi = gconf->l_max - 1; zi > -1; zi--) {
      for (int zj = 0; zj < 3; zj++) {
	for (int zk = 0; zk < 2; zk++) {
	  delete[] zpv[zi][zj][zk];
	}
	delete[] zpv[zi][zj];
      }
      delete[] zpv[zi];
    }
    delete[] zpv;
    delete[] duk;
    delete[] u;
    delete[] us;
    delete[] un;
    delete[] uns;
    delete[] up;
    delete[] ups;
    delete[] upmp;
    delete[] upsmp;
    delete[] unmp;
    delete[] unsmp;
    delete[] vint;
    delete[] argi;
    delete[] args;
    delete[] gaps;
    for (int i = 3; i > -1; i--) {
      delete[] cmul[i];
      delete[] cmullr[i];
    }
    delete[] cmul;
    delete[] cmullr;
    for (int ti = 1; ti > -1; ti--) {
      delete[] tqse[ti];
      delete[] tqss[ti];
      delete[] tqspe[ti];
      delete[] tqsps[ti];
    }
    delete[] tqse;
    delete[] tqss;
    delete[] tqspe;
    delete[] tqsps;
    printf("Finished: output written to %s.\n", (output_path + "/c_OSPH").c_str());
  } else { // NSPH mismatch between geometry and scatterer configurations.
    throw UnrecognizedConfigurationException(
					     "Inconsistent geometry and scatterer configurations."
					     );
  }
  delete sconf;
  delete gconf;
}
