/* Copyright 2024 INAF - Osservatorio Astronomico di Cagliari

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
 */

/*! \file sphere.cpp
 *
 * \brief Implementation of the single sphere calculation.
 */
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>

#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif

#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif

#ifndef INCLUDE_LOGGING_H_
#include "../include/logging.h"
#endif

#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_TRANSITIONMATRIX_H_
#include "../include/TransitionMatrix.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(const string& config_file, const string& data_file, const string& output_path) {
  Logger *logger = new Logger(LOG_INFO);
  dcomplex arg, s0, tfsas;
  double th, ph;
  logger->log("INFO: making legacy configuration...\n");
  ScattererConfiguration *sconf = NULL;
  try {
    sconf = ScattererConfiguration::from_dedfb(config_file);
  } catch(const OpenConfigurationFileException &ex) {
    logger->err("\nERROR: failed to open scatterer configuration file.\n");
    string message = ex.what();
    logger->err("FILE: " + message + "\n");
    delete logger;
    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) {
    logger->err("\nERROR: failed to open geometry configuration file.\n");
    string message = ex.what();
    logger->err("FILE: " + message + "\n");
    if (sconf != NULL) delete sconf;
    delete logger;
    exit(1);
  }
  int s_nsph = sconf->number_of_spheres;
  int nsph = gconf->number_of_spheres;
  if (s_nsph == nsph) {
    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];
    }
    dcomplex **tqspe, **tqsps;
    double **tqse, **tqss;
    tqse = new double*[2];
    tqss = new double*[2];
    tqspe = new dcomplex*[2];
    tqsps = new dcomplex*[2];
    for (int ti = 0; ti < 2; ti++) {
      tqse[ti] = new double[2]();
      tqss[ti] = new double[2]();
      tqspe[ti] = new dcomplex[2]();
      tqsps[ti] = new dcomplex[2]();
    }
    double frx = 0.0, fry = 0.0, frz = 0.0;
    double cfmp, cfsp, sfmp, sfsp;
    int jw;
    int l_max = gconf->l_max;
    C1 *c1 = new C1(gconf, sconf);
    int npnt = gconf->npnt;
    int npntts = gconf->npntts;
    int in_pol = gconf->in_pol;
    int meridional_type = gconf->iavm;
    int jwtm = gconf->jwtm;
    double in_theta_start = gconf->in_theta_start;
    double in_theta_step = gconf->in_theta_step;
    double in_theta_end = gconf->in_theta_end;
    double sc_theta_start = gconf->sc_theta_start;
    double sc_theta_step = gconf->sc_theta_step;
    double sc_theta_end = gconf->sc_theta_end;
    double in_phi_start = gconf->in_phi_start;
    double in_phi_step = gconf->in_phi_step;
    double in_phi_end = gconf->in_phi_end;
    double sc_phi_start = gconf->sc_phi_start;
    double sc_phi_step = gconf->sc_phi_step;
    double sc_phi_end = gconf->sc_phi_end;
    C2 *c2 = new C2(gconf, sconf);
    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",
	    nsph,
	    l_max,
	    in_pol,
	    npnt,
	    npntts,
	    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",
	    in_theta_start,
	    in_theta_step,
	    in_theta_end,
	    sc_theta_start,
	    sc_theta_step,
	    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",
	    in_phi_start,
	    in_phi_step,
	    in_phi_end,
	    sc_phi_start,
	    sc_phi_step,
	    sc_phi_end
	    );
    fprintf(output, " READ(IR,*)JWTM\n");
    fprintf(output, " %5d\n", 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 (in_theta_step != 0.0)
      nth = int((in_theta_end - in_theta_start) / in_theta_step + sml);
    nth += 1;
    if (in_phi_step != 0.0)
      nph = int((in_phi_end - in_phi_start) / in_phi_step + sml);
    nph += 1;
    int nths = 0, nphs = 0;
    double thsca;
    if (meridional_type > 1) { // ISAM > 1, fixed scattering angle
      nths = 1;
      thsca = sc_theta_start - in_theta_start;
    } else { //ISAM <= 1
      if (in_theta_step != 0.0)
	nths = int((sc_theta_end - sc_theta_start) / sc_theta_step + sml);
      nths += 1;
      if (meridional_type == 1) { // ISAM = 1
	nphs = 1;
      } else { // ISAM < 1
	if (sc_phi_step != 0.0)
	  nphs = int((sc_phi_end - sc_phi_start) / sc_phi_step + sml);
	nphs += 1;
      }
    }
    int nk = nth * nph;
    int nks = nths * nphs;
    int nkks = nk * nks;
    double th1 = in_theta_start;
    double ph1 = in_phi_start;
    double ths1 = sc_theta_start;
    double phs1 = 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->get_rcf(i116, j115) * c1->ros[i116];
	}
      }
      gcs += c1->gcsv[iogi - 1];
    }
    double ****zpv = new double***[l_max]; //[l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
    for (int zi = 0; zi < 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(l_max, zpv);
    double exdc = sconf->exdc;
    double exri = sqrt(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 *>(&(meridional_type)), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&(in_pol)), sizeof(int));
      tppoan.write(reinterpret_cast<char *>(&s_nsph), 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));
      tppoan.write(reinterpret_cast<char *>(&nsph), sizeof(int));

      for (int nsi = 0; nsi < nsph; nsi++)
	tppoan.write(reinterpret_cast<char *>(&(c1->iog[nsi])), sizeof(int));
      if (in_pol == 0) fprintf(output, "   LIN\n");
      else fprintf(output, "  CIRC\n");
      fprintf(output, " \n");
      double wp = sconf->wp;
      double xip = sconf->xip;
      double wn = wp / 3.0e8;
      double sqsfi = 1.0;
      double vk, vkarg;
      int idfc = sconf->idfc;
      int nxi = sconf->number_of_scales;
      if (idfc < 0) {
	vk = xip * wn;
	fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
	fprintf(output, " \n");
      }
      for (int jxi488 = 0; jxi488 < nxi; jxi488++) {
	int jxi = jxi488 + 1;
	logger->log("INFO: running scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
	fprintf(output, "========== JXI =%3d ====================\n", jxi);
	double xi = sconf->get_scale(jxi488);
	if (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 = c1->iog[i132];
	  if (iogi != i) {
	    for (int l123 = 0; l123 < 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 = c1->nshl[i132];
	  int ici = (nsh + 1) / 2;
	  if (idfc == 0) {
	    for (int ic = 0; ic < ici; ic++)
	      c2->dc0[ic] = sconf->get_dielectric_constant(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->get_dielectric_constant(ic, i132, jxi488);
	      }
	    }
	  }
	  if (nsh % 2 == 0) c2->dc0[ici] = exdc;
	  int jer = 0;
	  int lcalc = 0;
	  dme(l_max, i, npnt, npntts, vkarg, 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, real(arg), imag(arg));
	    tppoan.close();
	    fclose(output);
	    delete sconf;
	    delete gconf;
	    delete c1;
	    delete c2;
	    for (int zi = 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[] 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;
	    delete logger;
	    return;
	  }
	} // i132
	if (idfc >= 0 and nsph == 1 and jxi == jwtm) {
	  // This is the condition that writes the transition matrix to output.
	  string ttms_name = output_path + "/c_TTMS.hd5";
	  TransitionMatrix::write_binary(
					 ttms_name, l_max, vk, exri, c1->rmi, c1->rei,
					 sconf->get_radius(0), "HDF5"
					 );
	  ttms_name = output_path + "/c_TTMS";
	  TransitionMatrix::write_binary(
					 ttms_name, l_max, vk, exri, c1->rmi, c1->rei,
					 sconf->get_radius(0)
					 );
	}
	double cs0 = 0.25 * vk * vk * vk / half_pi;
	sscr0(tfsas, nsph, l_max, vk, exri, c1);
	double sqk = vk * vk * exdc;
	aps(zpv, l_max, nsph, c1, sqk, gaps);
	rabas(in_pol, 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],
		      real(c2->vkt[i170]),
		      imag(c2->vkt[i170])
		      );
	    }
	    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", real(c1->fsas[i170]), imag(c1->fsas[i170]));
	    double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
	    s0 = c1->fsas[i170] * exri;
	    double qschu = csch * imag(s0);
	    double pschu = csch * real(s0);
	    double s0mag = cs0 * cabs(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 = real(tqspe[0][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = imag(tqspe[0][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = real(tqsps[0][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = imag(tqsps[0][i170]);
	    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 = real(tqspe[1][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = imag(tqspe[1][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = real(tqsps[1][i170]);
	    tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
	    val = imag(tqsps[1][i170]);
	    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", real(tfsas), imag(tfsas));
	  double csch = 2.0 * vk * sqsfi / gcs;
	  s0 = tfsas * exri;
	  double qschu = csch * imag(s0);
	  double pschu = csch * real(s0);
	  double s0mag = cs0 * cabs(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 (meridional_type >= 0) {
	      wmamp(0, cost, sint, cosp, sinp, in_pol, 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 (meridional_type > 1) ths = th + thsca;
	      if (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 (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 (meridional_type >= 0) {
		    wmamp(2, costs, sints, cosps, sinps, in_pol, 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 (meridional_type < 0) {
		    wmasp(
			  cost, sint, cosp, sinp, costs, sints, cosps, sinps,
			  u, up, un, us, ups, uns, isq, ibf, in_pol,
			  l_max, 0, nsph, argi, args, c1
			  );
		  }
		  for (int i193 = 0; i193 < 3; i193++) {
		    un[i193] = unmp[i193];
		    uns[i193] = unsmp[i193];
		  }
		}
		if (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 (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, 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",
			  real(c1->sas[ns226][0][0]), imag(c1->sas[ns226][0][0]),
			  real(c1->sas[ns226][1][0]), imag(c1->sas[ns226][1][0])
			  );
		  fprintf(
			  output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
			  real(c1->sas[ns226][0][1]), imag(c1->sas[ns226][0][1]),
			  real(c1->sas[ns226][1][1]), imag(c1->sas[ns226][1][1])
			  );
		  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++) c1->vint[i225] = c1->vints[ns226][i225];
		  mmulc(c1->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 = real(c1->vint[vi]);
		    tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		    value = imag(c1->vint[vi]);
		    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 (meridional_type < 1) phs += sc_phi_step;
	      } // jphs480 loop
	      if (meridional_type <= 1) thsl += sc_theta_step;
	    } // jths482 loop
	    ph += in_phi_step;
	  } // jph484 loop on elevation
	  th += in_theta_step;
	} // jth486 loop on azimuth
	logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
      } //jxi488 loop on scales
      tppoan.close();
    } else { // In case TPPOAN could not be opened. Should never happen.
      logger->err("ERROR: failed to open TPPOAN file.\n");
    }
    fclose(output);
    delete c1;
    delete c2;
    for (int zi = 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[] 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;
    logger->log("Finished: output written to " + output_path + "/c_OSPH.\n");
  } else { // NSPH mismatch between geometry and scatterer configurations.
    throw UnrecognizedConfigurationException(
					     "Inconsistent geometry and scatterer configurations."
					     );
  }
  delete sconf;
  delete gconf;
  delete logger;
}
