Something went wrong on our end
-
Giovanni La Mura authoredGiovanni La Mura authored
sphere.cpp 22.36 KiB
/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
/*! \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(string config_file, string data_file, 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 = (int)sconf->get_param("nsph");
int nsph = (int)gconf->get_param("nsph");
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->get_param("l_max");
int *iog_vec = sconf->get_iog_vec();
int *nshl_vec = sconf->get_nshl();
C1 *c1 = new C1(nsph, l_max, nshl_vec, iog_vec);
for (int i = 0; i < nsph; i++) {
c1->ros[i] = sconf->get_radius(i);
}
int npnt = (int)gconf->get_param("npnt");
int npntts = (int)gconf->get_param("npntts");
int in_pol = (int)gconf->get_param("in_pol");
int meridional_type = (int)gconf->get_param("meridional_type");
int jwtm = (int)gconf->get_param("jwtm");
double in_theta_start = gconf->get_param("in_theta_start");
double in_theta_step = gconf->get_param("in_theta_step");
double in_theta_end = gconf->get_param("in_theta_end");
double sc_theta_start = gconf->get_param("sc_theta_start");
double sc_theta_step = gconf->get_param("sc_theta_step");
double sc_theta_end = gconf->get_param("sc_theta_end");
double in_phi_start = gconf->get_param("in_phi_start");
double in_phi_step = gconf->get_param("in_phi_step");
double in_phi_end = gconf->get_param("in_phi_end");
double sc_phi_start = gconf->get_param("sc_phi_start");
double sc_phi_step = gconf->get_param("sc_phi_step");
double sc_phi_end = gconf->get_param("sc_phi_end"); C2 *c2 = new C2(nsph, 5, npnt, 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",
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];
}
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->get_param("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 *>(&(iog_vec[nsi])), sizeof(int));
if (in_pol == 0) fprintf(output, " LIN\n");
else fprintf(output, " CIRC\n");
fprintf(output, " \n");
double wp = sconf->get_param("wp");
double xip = sconf->get_param("xip");
double wn = wp / 3.0e8;
double sqsfi = 1.0;
double vk, vkarg;
int idfc = (int)sconf->get_param("idfc");
int nxi = (int)sconf->get_param("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 = iog_vec[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 = nshl_vec[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.
TransitionMatrix ttms(l_max, vk, exri, c1->rmi, c1->rei, sconf->get_radius(0));
string ttms_name = output_path + "/c_TTMS.hd5";
ttms.write_binary(ttms_name, "HDF5");
ttms_name = output_path + "/c_TTMS";
ttms.write_binary(ttms_name);
}
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;
}