Skip to content
Snippets Groups Projects
Select Git revision
  • 843c9960ebfae19d2f832e70ddc1ef442d7d5ea0
  • master default protected
  • ia2
  • adql2.1-ia2
  • private_rows
5 results

UnsupportedURIProtocolException.java

Blame
  • sphere.cpp 22.64 KiB
    /* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
    
       This program is free software: you can redistribute it and/or modify
       it under the terms of the GNU General Public License as published by
       the Free Software Foundation, either version 3 of the License, or
       (at your option) any later version.
       
       This program is distributed in the hope that it will be useful,
       but WITHOUT ANY WARRANTY; without even the implied warranty of
       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
       GNU General Public License for more details.
       
       A copy of the GNU General Public License is distributed along with
       this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
     */
    
    /*! \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;
        ParticleDescriptor *c1 = new ParticleDescriptorSphere(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
    	  }
    	  // label 125
    	  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, jxi488); // 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;
    }