/* 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; }