Something went wrong on our end
-
Giovanni La Mura authoredGiovanni La Mura authored
cluster.cpp 53.73 KiB
/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
/*! \file cluster.cp
*
* \brief Implementation of the calculation for a cluster of spheres.
*/
#include <chrono>
#include <cstdio>
#include <exception>
#include <fstream>
#include <string>
#ifdef _OPENMP
#include <omp.h>
#endif
#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_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif
#ifndef INCLUDE_TRANSITIONMATRIX_H_
#include "../include/TransitionMatrix.h"
#endif
#ifndef INCLUDE_ALGEBRAIC_H_
#include "../include/algebraic.h"
#endif
using namespace std;
// I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int nth, int nths, int nph, int nphs, int nk, int nks, int nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double th, double thstp, double thlst, double ths, double thsstp, double thslst, double ph, double phstp, double phlst, double phs, double phsstp, double phslst, double th1, double ph1, double ths1, double phs1, double thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, int isq, int ibf, Logger *logger);
/*! \brief C++ implementation of CLU
*
* \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 cluster(const string& config_file, const string& data_file, const string& output_path) {
chrono::time_point<chrono::high_resolution_clock> t_start = chrono::high_resolution_clock::now();
string timing_name = output_path + "/c_timing.log";
FILE *timing_file = fopen(timing_name.c_str(), "w");
Logger *time_logger = new Logger(LOG_DEBG, timing_file);
Logger *logger = new Logger(LOG_INFO);
logger->log("INFO: making legacy configuration...", LOG_INFO); 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 = "FILE: " + string(ex.what()) + "\n";
logger->err(message);
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 = "FILE: " + string(ex.what()) + "\n";
logger->err(message);
if (sconf) delete sconf;
exit(1);
}
logger->log(" done.\n", LOG_INFO);
int s_nsph = sconf->get_param("nsph");
int nsph = gconf->get_param("nsph");
if (s_nsph == nsph) {
// Shortcuts to variables stored in configuration objects
np_int mxndm = (np_int)gconf->get_param("mxndm");
int inpol = (int)gconf->get_param("in_pol");
int npnt = (int)gconf->get_param("npnt");
int npntts = (int)gconf->get_param("npntts");
int iavm = (int)gconf->get_param("iavm");
int isam = (int)gconf->get_param("meridional_type");
int nxi = (int)sconf->get_param("number_of_scales");
int idfc = (int)sconf->get_param("idfc");
double th = gconf->get_param("in_theta_start");
double thstp = gconf->get_param("in_theta_step");
double thlst = gconf->get_param("in_theta_end");
double ths = gconf->get_param("sc_theta_start");
double thsstp = gconf->get_param("sc_theta_step");
double thslst = gconf->get_param("sc_theta_end");
double ph = gconf->get_param("in_phi_start");
double phstp = gconf->get_param("in_phi_step");
double phlst = gconf->get_param("in_phi_end");
double phs = gconf->get_param("sc_phi_start");
double phsstp = gconf->get_param("sc_phi_step");
double phslst = gconf->get_param("sc_phi_end");
double wp = sconf->get_param("wp");
// Global variables for CLU
int lm = (int)gconf->get_param("l_max");
int li = (int)gconf->get_param("li");
int le = (int)gconf->get_param("le");
if (li > lm) lm = li;
if (le > lm) lm = le;
C1 *c1 = new C1(gconf, sconf);
C3 *c3 = new C3();
C4 *c4 = new C4(gconf);
C1_AddOns *c1ao = new C1_AddOns(c4);
// End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0;
int lcalc = 0;
dcomplex arg = 0.0 + 0.0 * I;
dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = (int)sconf->get_param("configurations");
C2 *c2 = new C2(gconf, sconf);
np_int ndit = 2 * nsph * c4->nlim;
logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n"); const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph]();
double *tqev = new double[3]();
double *tqsv = new double[3]();
double **tqse, **tqss, **tqce, **tqcs;
dcomplex **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2];
tqspe = new dcomplex*[2];
tqss = new double*[2];
tqsps = new dcomplex*[2];
tqce = new double*[2];
tqcpe = new dcomplex*[2];
tqcs = new double*[2];
tqcps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph]();
tqspe[ti] = new dcomplex[nsph]();
tqss[ti] = new double[nsph]();
tqsps[ti] = new dcomplex[nsph]();
tqce[ti] = new double[3]();
tqcpe[ti] = new dcomplex[3]();
tqcs[ti] = new double[3]();
tqcps[ti] = new dcomplex[3]();
}
double *gapv = new double[3]();
dcomplex **gapp, **gappm;
double **gap, **gapm;
gapp = new dcomplex*[3];
gappm = new dcomplex*[3];
gap = new double*[3];
gapm = new double*[3];
for (int gi = 0; gi < 3; gi++) {
gapp[gi] = new dcomplex[2]();
gappm[gi] = new dcomplex[2]();
gap[gi] = new double[2]();
gapm[gi] = new double[2]();
}
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 *unmp = new double[3]();
double *unsmp = new double[3]();
double *upmp = new double[3]();
double *upsmp = new double[3]();
double *argi = new double[1]();
double *args = new double[1]();
double *duk = new double[3]();
double **cextlr, **cext;
double **cmullr, **cmul;
cextlr = new double*[4];
cext = new double*[4];
cmullr = new double*[4];;
cmul = new double*[4];
for (int ci = 0; ci < 4; ci++) {
cextlr[ci] = new double[4]();
cext[ci] = new double[4]();
cmullr[ci] = new double[4]();
cmul[ci] = new double[4]();
}
int isq, ibf;
int jwtm = (int)gconf->get_param("jwtm");
double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam );
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
);
fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
th, thstp, thlst, ths, thsstp, thslst
);
fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
ph, phstp, phlst, phs, phsstp, phslst
);
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");
fprintf(output, " \n");
double small = 1.0e-3;
int nth = 0, nph = 0;
if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
nth++;
if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
nph++;
int nths = 0, nphs = 0;
double thsca = 0.0;
if (isam > 1) {
nths = 1;
thsca = ths - th;
} else { // ISAM <= 1
if (thsstp == 0.0) nths = 0;
else nths = int ((thslst - ths) / thsstp + small);
nths++;
}
if (isam >= 1) {
nphs = 1;
} else {
if (phsstp == 0.0) nphs = 0;
else nphs = int((phslst - phs) / phsstp + small);
nphs++;
}
int nk = nth * nph;
int nks = nths * nphs;
int nkks = nk * nks;
double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
str(sconf, c1, c1ao, c3, c4, c6);
double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
for (int zi = 0; zi < c4->lm; 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(c4->lm, zpv);
double exdc = sconf->get_param("exdc");
double exri = sqrt(exdc);
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream *tppoanp = new fstream;
fstream &tppoan = *tppoanp;
string tppoan_name = output_path + "/c_TPPOAN"; tppoan.open(tppoan_name.c_str(), ios::out | ios::binary);
if (tppoan.is_open()) {
#ifdef USE_LAPACK
logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
#else
logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
#endif
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nxi), 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));
double wn = wp / 3.0e8;
double xip = sconf->get_param("xip");
double sqsfi = 1.0;
if (idfc < 0) {
vk = xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
// do the first iteration on jxi488 separately, since it seems to be different from the others
int jxi488 = 1;
jer = cluster_jxi488_cycle(jxi488, sconf, gconf, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, nxi, nsph, mxndm, inpol, iavm, npnt, npntts, isam, lm, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, isq, ibf, logger);
// Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
int ompnumthreads = 1;
#pragma omp parallel
{
// Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
int myompthread = 0;
#ifdef _OPENMP
// If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files
myompthread = omp_get_thread_num();
if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
ScattererConfiguration *sconf_2 = NULL;
GeometryConfiguration *gconf_2 = NULL;
C1 *c1_2 = NULL;
C1_AddOns *c1ao_2 = NULL;
C2 *c2_2 = NULL;
C3 *c3_2 = NULL;
C4 *c4_2 = NULL;
C6 *c6_2 = NULL;
C9 *c9_2 = NULL;
FILE *output_2 = NULL;
fstream *tppoanp_2 = NULL;
double *gaps_2 = NULL;
double **tqse_2 = NULL;
double **tqce_2 = NULL;
double **tqcs_2 = NULL;
dcomplex **tqspe_2 = NULL;
dcomplex **tqcpe_2 = NULL;
dcomplex **tqcps_2 = NULL;
double **tqss_2 = NULL;
dcomplex **tqsps_2 = NULL;
double ****zpv_2 = NULL;
double **gapm_2 = NULL;
dcomplex **gappm_2 = NULL;
double **gap_2 = NULL;
dcomplex **gapp_2 = NULL;
double *argi_2 = NULL;
double *args_2 = NULL;
double *duk_2 = NULL;
double **cextlr_2 = NULL;
double **cext_2 = NULL; double **cmullr_2 = NULL;
double **cmul_2 = NULL;
double *gapv_2 = NULL;
double *tqev_2 = NULL;
double *tqsv_2 = NULL;
int nxi_2 = nxi;
int nsph_2 = nsph;
np_int mxndm_2 = mxndm;
int inpol_2 = inpol;
int iavm_2 = iavm;
int npnt_2 = npnt;
int npntts_2 = npntts;
int isam_2 = isam;
int lm_2 = lm;
double th_2 = th;
double thstp_2 = thstp;
double thlst_2 = thlst;
double ths_2 = ths;
double thsstp_2 = thsstp;
double thslst_2 = thslst;
double ph_2 = ph;
double phstp_2 = phstp;
double phlst_2 = phlst;
double phs_2 = phs;
double phsstp_2 = phsstp;
double phslst_2 = phslst;
double th1_2 = th1;
double ph1_2 = ph1;
double ths1_2 = ths1;
double phs1_2 = phs1;
double thsca_2 = thsca;
double *u_2 = NULL;
double *us_2 = NULL;
double *un_2 = NULL;
double *uns_2 = NULL;
double *up_2 = NULL;
double *ups_2 = NULL;
double *unmp_2 = NULL;
double *unsmp_2 = NULL;
double *upmp_2 = NULL;
double *upsmp_2 = NULL;
double scan_2 = scan;
double cfmp_2 = cfmp;
double sfmp_2 = sfmp;
double cfsp_2 = cfsp;
double sfsp_2 = sfsp;
double sqsfi_2 = sqsfi;
double exri_2 = exri;
int lcalc_2 = lcalc;
dcomplex arg_2 = arg;
double wn_2 = wn;
double vk_2 = vk;
np_int ndit_2 = ndit;
int isq_2 = isq;
int ibf_2 = ibf;
int nth_2 = nth;
int nths_2 = nths;
int nph_2 = nph;
int nphs_2 = nph;
int nk_2 = nk;
int nks_2 = nks;
int nkks_2 = nkks;
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
if (myompthread == 0) {
sconf_2 = sconf;
gconf_2 = gconf;
c1_2 = c1;
c1ao_2 = c1ao;
c2_2 = c2;
c3_2 = c3; c4_2 = c4;
c6_2 = c6;
c9_2 = c9;
output_2 = output;
tppoanp_2 = tppoanp;
gaps_2 = gaps;
tqse_2 = tqse;
tqce_2 = tqce;
tqcs_2 = tqcs;
tqspe_2 = tqspe;
tqcpe_2 = tqcpe;
tqcps_2 = tqcps;
tqss_2 = tqss;
tqsps_2 = tqsps;
zpv_2 = zpv;
gapm_2 = gapm;
gappm_2 = gappm;
gap_2 = gap;
gapp_2 = gapp;
argi_2 = argi;
args_2 = args;
duk_2 = duk;
cextlr_2 = cextlr;
cext_2 = cext;
cmullr_2 = cmullr;
cmul_2 = cmul;
gapv_2 = gapv;
tqev_2 = tqev;
tqsv_2 = tqsv;
u_2 = u;
us_2 = us;
un_2 = un;
uns_2 = uns;
up_2 = up;
ups_2 = ups;
unmp_2 = unmp;
unsmp_2 = unsmp;
upmp_2 = upmp;
upsmp_2 = upsmp;
}
else {
// this is not thread 0, so do create fresh copies of all local variables
sconf_2 = new ScattererConfiguration(*sconf);
gconf_2 = new GeometryConfiguration(*gconf);
c1_2 = new C1(*c1);
c1ao_2 = new C1_AddOns(*c1ao);
c2_2 = new C2(*c2);
c3_2 = new C3(*c3);
c4_2 = new C4(*c4);
c6_2 = new C6(*c6);
c9_2 = new C9(*c9);
output_2 = fopen((output_path + "/c_OCLU_" + to_string(myompthread)).c_str(), "w");
tppoanp_2 = new fstream;
tppoanp_2->open((output_path + "/c_TPPOAN_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
gaps_2 = new double[nsph]();
tqse_2 = new double*[2];
tqce_2 = new double*[2];
tqcs_2 = new double*[2];
tqspe_2 = new dcomplex*[2];
tqcpe_2 = new dcomplex*[2];
tqcps_2 = new dcomplex*[2];
tqss_2 = new double*[2];
tqsps_2 = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse_2[ti] = new double[nsph]();
tqce_2[ti] = new double[3]();
tqcs_2[ti] = new double[3]();
tqspe_2[ti] = new dcomplex[nsph]();
tqcpe_2[ti] = new dcomplex[3]();
tqcps_2[ti] = new dcomplex[3](); tqss_2[ti] = new double[nsph]();
tqsps_2[ti] = new dcomplex[nsph]();
for (int tj=0; tj<nsph; tj++) {
tqse_2[ti][tj] = tqse[ti][tj];
tqspe_2[ti][tj] = tqspe[ti][tj];
tqss_2[ti][tj] = tqss[ti][tj];
tqsps_2[ti][tj] = tqsps[ti][tj];
}
for (int tj=0; tj<3; tj++) {
tqce_2[ti][tj] = tqce[ti][tj];
tqcs_2[ti][tj] = tqcs[ti][tj];
tqcpe_2[ti][tj] = tqcpe[ti][tj];
tqcps_2[ti][tj] = tqcps[ti][tj];
}
}
zpv_2 = new double***[c4->lm]; // [lm][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
for (int zi = 0; zi < c4->lm; zi++) {
zpv_2[zi] = new double**[3];
for (int zj = 0; zj < 3; zj++) {
zpv_2[zi][zj] = new double*[2];
for (int zk = 0; zk < 2; zk++) {
zpv_2[zi][zj][zk] = new double[2]();
for (int zl = 0; zl < 2; zl++) zpv_2[zi][zj][zk][zl] = zpv[zi][zj][zk][zl];
}
}
}
gapm_2 = new double*[3];
gappm_2 = new dcomplex*[3];
gap_2 = new double*[3];
gapp_2 = new dcomplex*[3];
for (int gi = 0; gi < 3; gi++) {
gap_2[gi] = new double[2]();
gapm_2[gi] = new double[2]();
gapp_2[gi] = new dcomplex[2]();
gappm_2[gi] = new dcomplex[2]();
for (int gj=0; gj<2; gj++) {
gap_2[gi][gj] = gap[gi][gj];
gapp_2[gi][gj] = gapp[gi][gj];
gapm_2[gi][gj] = gapm[gi][gj];
gappm_2[gi][gj] = gappm[gi][gj];
}
}
argi_2 = new double[1]();
argi_2[0] = argi[0];
args_2 = new double[1]();
args_2[0] = args[0];
duk_2 = new double[3]();
for (int di=0; di<3; di++) duk_2[di] = duk[di];
cextlr_2 = new double*[4];
cext_2 = new double*[4];
cmullr_2 = new double*[4];
cmul_2 = new double*[4];
for (int ci = 0; ci < 4; ci++) {
cextlr_2[ci] = new double[4]();
cext_2[ci] = new double[4]();
cmullr_2[ci] = new double[4]();
cmul_2[ci] = new double[4]();
for (int cj=0; cj<4; cj++) {
cextlr_2[ci][cj] = cextlr[ci][cj];
cext_2[ci][cj] = cext[ci][cj];
cmullr_2[ci][cj] = cmullr[ci][cj];
cmul_2[ci][cj] = cmul[ci][cj];
}
}
gapv_2 = new double[3]();
for (int gi=0; gi<3; gi++) gapv_2[gi] = gapv[gi];
tqev_2 = new double[3]();
tqsv_2 = new double[3]();
for (int ti=0; ti<3; ti++) {
tqev_2[ti] = tqev[ti]; tqsv_2[ti] = tqsv[ti];
}
u_2 = new double[3]();
us_2 = new double[3]();
un_2 = new double[3]();
uns_2 = new double[3]();
up_2 = new double[3]();
ups_2 = new double[3]();
unmp_2 = new double[3]();
unsmp_2 = new double[3]();
upmp_2 = new double[3]();
upsmp_2 = new double[3]();
for (int ti=0; ti<3; ti++) {
u_2[ti] = u[ti];
us_2[ti] = us[ti];
un_2[ti] = un[ti];
uns_2[ti] = uns[ti];
up_2[ti] = up[ti];
ups_2[ti] = ups[ti];
unmp_2[ti] = unmp[ti];
unsmp_2[ti] = unsmp[ti];
upmp_2[ti] = upmp[ti];
upsmp_2[ti] = upsmp[ti];
}
}
fstream &tppoan_2 = *tppoanp_2;
// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
#pragma omp barrier
if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
// ok, now I can actually start the parallel calculations
#pragma omp for
for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output_2, output_path, gaps_2, tqse_2, tqspe_2, tqss_2, tqsps_2, zpv_2, gapm_2, gappm_2, nth_2, nths_2, nph_2, nphs_2, nk_2, nks_2, nkks_2, argi_2, args_2, gap_2, gapp_2, tqce_2, tqcpe_2, tqcs_2, tqcps_2, duk_2, tppoan_2, cextlr_2, cext_2, cmullr_2, cmul_2, gapv_2, tqev_2, tqsv_2, nxi_2, nsph_2, mxndm_2, inpol_2, iavm_2, npnt_2, npntts_2, isam_2, lm_2, th_2, thstp_2, thlst_2, ths_2, thsstp_2, thslst_2, ph_2, phstp_2, phlst_2, phs_2, phsstp_2, phslst_2, th1_2, ph1_2, ths1_2, phs1_2, thsca_2, u_2, us_2, un_2, uns_2, up_2, ups_2, unmp_2, unsmp_2, upmp_2, upsmp_2, scan_2, cfmp_2, sfmp_2, cfsp_2, sfsp_2, sqsfi_2, exri_2, lcalc_2, arg_2, wn_2, vk_2, ndit_2, isq_2, ibf_2, logger);
}
#pragma omp barrier
// only threads different from 0 have to free local copies of variables and close local files
if (myompthread != 0) {
delete sconf_2;
delete gconf_2;
delete c1_2;
delete c1ao_2;
delete c2_2;
delete c3_2;
delete c4_2;
delete c6_2;
delete c9_2;
fclose(output_2);
tppoanp_2->close();
delete tppoanp_2;
delete[] gaps_2;
for (int ti = 0; ti < 2; ti++) {
delete[] tqse_2[ti];
delete[] tqce_2[ti];
delete[] tqcs_2[ti];
delete[] tqspe_2[ti];
delete[] tqcpe_2[ti];
delete[] tqcps_2[ti];
delete[] tqss_2[ti];
delete[] tqsps_2[ti];
}
delete[] tqse_2;
delete[] tqce_2;
delete[] tqcs_2;
delete[] tqspe_2;
delete[] tqcpe_2;
delete[] tqcps_2;
delete[] tqss_2;
delete[] tqsps_2;
for (int zi = 0; zi < c4->lm; zi++) { for (int zj = 0; zj < 3; zj++) {
for (int zk = 0; zk < 2; zk++) {
delete[] zpv_2[zi][zj][zk];
}
delete[] zpv_2[zi][zj];
}
delete[] zpv_2[zi];
}
delete[] zpv_2;
for (int gi = 0; gi < 3; gi++) {
delete[] gap_2[gi];
delete[] gapp_2[gi];
delete[] gapm_2[gi];
delete[] gappm_2[gi];
}
delete[] gap_2;
delete[] gapp_2;
delete[] gapm_2;
delete[] gappm_2;
delete[] argi_2;
delete[] args_2;
delete[] duk_2;
for (int ci = 0; ci < 4; ci++) {
delete[] cextlr_2[ci];
delete[] cext_2[ci];
delete[] cmullr_2[ci];
delete[] cmul_2[ci];
}
delete[] cextlr_2;
delete[] cext_2;
delete[] cmullr_2;
delete[] cmul_2;
delete[] gapv_2;
delete[] tqev_2;
delete[] tqsv_2;
delete[] u_2;
delete[] us_2;
delete[] un_2;
delete[] uns_2;
delete[] up_2;
delete[] ups_2;
delete[] unmp_2;
delete[] unsmp_2;
delete[] upmp_2;
delete[] upsmp_2;
}
#pragma omp barrier
{
string message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
logger->log(message);
}
} // closes pragma omp parallel
#ifdef _OPENMP
#pragma omp barrier
{
// thread 0 already wrote on global files, skip it and take care of appending the others
for (int ri = 1; ri < ompnumthreads; ri++) {
// Giovanni, please add here in this loop code to reopen the temporary files, reread them and append them respectively to the global output and tppoan, before closing them
string partial_file_name = output_path + "/c_OCLU_" + to_string(ri);
string message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
logger->log(message);
FILE *partial_output = fopen(partial_file_name.c_str(), "r");
char c = fgetc(partial_output);
while (c != EOF) {
fputc(c, output);
c = fgetc(partial_output);
}
fclose(partial_output);
remove(partial_file_name.c_str());
logger->log("done.\n"); partial_file_name = output_path + "/c_TPPOAN_" + to_string(ri);
message = "Copying binary output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
logger->log(message);
fstream partial_tppoan;
partial_tppoan.open(partial_file_name.c_str(), ios::in | ios::binary);
partial_tppoan.seekg(0, ios::end);
long buffer_size = partial_tppoan.tellg();
char *binary_buffer = new char[buffer_size];
partial_tppoan.seekg(0, ios::beg);
partial_tppoan.read(binary_buffer, buffer_size);
tppoan.write(binary_buffer, buffer_size);
partial_tppoan.close();
delete[] binary_buffer;
remove(partial_file_name.c_str());
logger->log("done.\n");
}
}
#endif
tppoanp->close();
delete tppoanp;
} else { // In case TPPOAN could not be opened. Should never happen.
logger->err("\nERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
for (int zi = c4->lm - 1; zi > -1; zi--) {
for (int zj = 2; zj > -1; zj--) {
delete[] zpv[zi][zj][1];
delete[] zpv[zi][zj][0];
delete[] zpv[zi][zj];
}
delete[] zpv[zi];
}
delete[] zpv;
delete c1;
delete c1ao;
delete c2;
delete c3;
delete c4;
delete c6;
delete c9;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
}
delete[] tqse;
delete[] tqss;
delete[] tqspe;
delete[] tqsps;
delete[] tqce;
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
delete[] tqev;
delete[] tqsv;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
delete[] gap[gi];
delete[] gapm[gi];
}
delete[] gapp;
delete[] gappm; delete[] gap;
delete[] gapm;
delete[] gapv;
delete[] u;
delete[] us;
delete[] un;
delete[] uns;
delete[] up;
delete[] ups;
delete[] unmp;
delete[] unsmp;
delete[] upmp;
delete[] upsmp;
delete[] argi;
delete[] args;
delete[] duk;
for (int ci = 3; ci > -1; ci--) {
delete[] cextlr[ci];
delete[] cext[ci];
delete[] cmullr[ci];
delete[] cmul[ci];
}
delete[] cextlr;
delete[] cext;
delete[] cmullr;
delete[] cmul;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
);
}
delete sconf;
delete gconf;
chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();
const chrono::duration<double> elapsed = t_end - t_start;
string message = "Calculation lasted " + to_string(elapsed.count()) + ".\n";
logger->log(message);
logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
time_logger->log(message);
fclose(timing_file);
delete time_logger;
delete logger;
}
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, int nth, int nths, int nph, int nphs, int nk, int nks, int nkks, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double th, double thstp, double thlst, double ths, double thsstp, double thslst, double ph, double phstp, double phlst, double phs, double phsstp, double phslst, double th1, double ph1, double ths1, double phs1, double thsca, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, int isq, int ibf, Logger *logger)
{
logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
int jer = 0;
int jaw = 1;
int jwtm = (int)gconf->get_param("jwtm");
int li = (int)gconf->get_param("li");
int le = (int)gconf->get_param("le");
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->get_param("exdc");
int idfc = (int)sconf->get_param("idfc");
double vkarg = 0.0;
if (idfc >= 0) {
vk = xi * wn;
vkarg = vk;
fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
} else {
vkarg = xi * vk;
sqsfi = 1.0 / (xi * xi);
fprintf(output, " XI=%15.7lE\n", xi);
}
hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
if (jer != 0) {
fprintf(output, " STOP IN HJV\n"); return jer;
// break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
}
for (int i132 = 1; i132 <= nsph; i132++) {
int iogi = c1->iog[i132 - 1];
if (iogi != i132) {
for (int l123 = 1; l123 <= li; l123++) {
c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
} // l123 loop
} else {
int nsh = c1->nshl[i132 - 1];
int ici = (nsh + 1) / 2;
if (idfc == 0) {
for (int ic = 0; ic < ici; ic++)
c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
} else {
if (jxi488 == 1) {
for (int ic = 0; ic < ici; ic++)
c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
}
}
if (nsh % 2 == 0) c2->dc0[ici] = exdc;
dme(
c4->li, i132, npnt, npntts, vkarg, exdc, exri,
c1, c2, jer, lcalc, arg
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
// delete[] u;
// delete[] us;
// delete[] un;
// delete[] uns;
// delete[] up;
// delete[] ups;
// delete[] unmp;
// delete[] unsmp;
// delete[] upmp;
// delete[] upsmp;
// delete[] am_vector;
// delete[] am;
return jer;
//break;
}
}
if (jer != 0) {
// delete[] u;
// delete[] us;
// delete[] un;
// delete[] uns;
// delete[] up;
// delete[] ups;
// delete[] unmp;
// delete[] unsmp;
// delete[] upmp;
// delete[] upsmp;
// delete[] am_vector;
// delete[] am;
return jer;
//break;
}
} // i132 loop
dcomplex *am_vector = new dcomplex[ndit * ndit]();
dcomplex **am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) {
am[ai] = (am_vector + ai * ndit);
}
cms(am, c1, c1ao, c4, c6);
invert_matrix(am, ndit, jer, mxndm);
if (jer != 0) { // delete[] u;
// delete[] us;
// delete[] un;
// delete[] uns;
// delete[] up;
// delete[] ups;
// delete[] unmp;
// delete[] unsmp;
// delete[] upmp;
// delete[] upsmp;
delete[] am_vector;
delete[] am;
return jer;
// break; // jxi488 loop: goes to memory clean
}
ztm(am, c1, c1ao, c4, c6, c9);
delete[] am_vector;
delete[] am;
if (idfc >= 0) {
if (jxi488 == jwtm) {
int nlemt = 2 * c4->nlem;
string ttms_name = output_path + "/c_TTMS.hd5";
TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m, "HDF5");
ttms_name = output_path + "/c_TTMS";
TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m);
}
}
// label 156: continue from here
if (inpol == 0) {
fprintf(output, " LIN\n");
} else { // label 158
fprintf(output, " CIRC\n");
}
// label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0);
double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
dcomplex s0 = 0.0 + 0.0 * I;
scr0(vk, exri, c1, c1ao, c3, c4);
double sqk = vk * vk * exdc;
aps(zpv, c4->li, nsph, c1, sqk, gaps);
rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
for (int i170 = 1; i170 <= nsph; i170++) {
if (c1->iog[i170 - 1] >= i170) {
int i = i170 - 1;
double albeds = c1->sscs[i] / c1->sexs[i];
c1->sqscs[i] *= sqsfi;
c1->sqabs[i] *= sqsfi;
c1->sqexs[i] *= sqsfi;
fprintf(output, " SPHERE %2d\n", i170);
if (c1->nshl[i] != 1) {
fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
} else { // label 162
fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], real(c2->vkt[i]), imag(c2->vkt[i]));
}
// label 164
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i]), imag(c1->fsas[i]));
csch = 2.0 * vk * sqsfi / c1->gcsv[i];
s0 = c1->fsas[i] * exri;
qschu = imag(s0) * csch;
pschu = real(s0) * csch;
s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
double rapr = c1->sexs[i] - gaps[i];
double cosav = gaps[i] / c1->sscs[i];
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr); fprintf(output, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[0][i], tqss[0][i]);
fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
}
} // i170 loop
fprintf(output, " FSAT=%15.7lE%15.7lE\n", real(c3->tfsas), imag(c3->tfsas));
csch = 2.0 * vk * sqsfi / c3->gcs;
s0 = c3->tfsas * exri;
qschu = imag(s0) * csch;
pschu = real(s0) * csch;
s0mag = cabs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
th = th1;
for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
ph = ph1;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
for (int jph484 = 1; jph484 <= nph; jph484++) {
int jw = 0;
if (nk != 1 || jxi488 <= 1) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
if (isam >= 0) {
wmamp(
0, cost, sint, cosp, sinp, inpol, c4->le, 0,
nsph, argi, u, upmp, unmp, c1
);
// label 182
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
} else { // label 180, NK == 1 AND JXI488 == 1
if (isam >= 0) {
// label 182
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
}
// label 184
double thsl = ths1;
double phsph = 0.0;
for (int jths = 1; jths <= nths; jths++) {
ths = thsl;
int icspnv = 0;
if (isam > 1) ths += thsca;
if (isam >= 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;
}
// label 186
phs = phs1;
for (int jphs = 1; jphs <= nphs; jphs++) {
double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
if (isam >= 1) {
phs = ph + phsph;
if (phs > 360.0) phs -= 360.0;
}
// label 188
bool goto190 = (nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
if (isam >= 0)
wmamp(
2, costs, sints, cosps, sinps, inpol, c4->le,
0, nsph, args, us, upsmp, unsmp, c1 );
}
// label 190
if (nkks != 1 || jxi488 <= 1) {
upvsp(
u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
);
if (isam < 0) {
wmasp(
cost, sint, cosp, sinp, costs, sints, cosps, sinps,
u, up, un, us, ups, uns, isq, ibf, inpol, c4->le,
0, nsph, argi, args, c1
);
} else { // label 192
for (int i193 = 0; i193 < 3; i193++) {
up[i193] = upmp[i193];
un[i193] = unmp[i193];
ups[i193] = upsmp[i193];
uns[i193] = unsmp[i193];
}
}
}
// label 194
if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
if (isam < 0) {
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
// label 196
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 (jaw != 0) {
jaw = 0;
mextc(vk, exri, c1ao->fsacm, cextlr, cext);
// We now have some implicit loops writing to binary
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
double value = cext[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
double value = c1ao->scscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gapm[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
} fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
int jlr = 2;
for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
int ipol = (ilr210 % 2 == 0) ? 1 : -1;
if (ilr210 == 2) jlr = 1;
double extsm = c1ao->ecscm[ilr210 - 1];
double qextm = extsm * sqsfi / c3->gcs;
double extrm = extsm / c3->ecs;
double scasm = c1ao->scscm[ilr210 - 1];
double albdm = scasm / extsm;
double qscam = scasm *sqsfi / c3->gcs;
double scarm = scasm / c3->scs;
double abssm = extsm - scasm;
double qabsm = abssm * sqsfi / c3->gcs;
double absrm = abssm / c3->acs;
double acsecs = c3->acs / c3->ecs;
if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
dcomplex s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
double qschum = imag(s0m) * csch;
double pschum = real(s0m) * csch;
double s0magm = cabs(s0m) * cs0;
double rfinrm = real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(c3->tfsas);
double extcrm = imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 206
fprintf(output, " CIRC %2d\n", ipol);
}
// label 208
fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
ilr210, ilr210, real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
real(c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(c1ao->fsacm[jlr - 1][ilr210 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
);
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
double rapr = c1ao->ecscm[ilr210 - 1] - gapm[2][ilr210 - 1];
double cosav = gapm[2][ilr210 - 1] / c1ao->scscm[ilr210 - 1];
double fz = rapr;
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " Fk=%15.7lE\n", fz);
} // ilr210 loop
double rmbrif = (real(c1ao->fsacm[0][0]) - real(c1ao->fsacm[1][1])) / real(c1ao->fsacm[0][0]);
double rmdchr = (imag(c1ao->fsacm[0][0]) - imag(c1ao->fsacm[1][1])) / imag(c1ao->fsacm[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
}
// label 212
fprintf(output, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, 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 (isam >= 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 { // label 214
fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
}
// label 220 if (inpol == 0) {
fprintf(output, " LIN\n");
} else { // label 222
fprintf(output, " CIRC\n");
}
// label 224
scr2(vk, vkarg, exri, duk, c1, c1ao, c3, c4);
if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=MIN0(LI,LE)\n");
for (int i226 = 1; i226 <= nsph; i226++) {
if (c1->iog[i226 - 1] >= i226) {
fprintf(output, " SPHERE %2d\n", i226);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
real(c1->sas[i226 - 1][0][0]), imag(c1->sas[i226 - 1][0][0]),
real(c1->sas[i226 - 1][1][0]), imag(c1->sas[i226 - 1][1][0])
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
real(c1->sas[i226 - 1][0][1]), imag(c1->sas[i226 - 1][0][1]),
real(c1->sas[i226 - 1][1][1]), imag(c1->sas[i226 - 1][1][1])
);
for (int j225 = 0; j225 < 16; j225++) {
c1->vint[j225] = c1->vints[i226 - 1][j225];
} // j225 loop
mmulc(c1->vint, cmullr, cmul);
fprintf(output, " MULS\n");
for (int i1 = 0; i1 < 4; i1++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmul[i1][0], cmul[i1][1], cmul[i1][2], cmul[i1][3]
);
} // i1 loop
fprintf(output, " MULSLR\n");
for (int i1 = 0; i1 < 4; i1++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmullr[i1][0], cmullr[i1][1], cmullr[i1][2], cmullr[i1][3]
);
} // i1 loop
}
} // i226 loop
fprintf(
output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
real(c3->tsas[0][0]), imag(c3->tsas[0][0]),
real(c3->tsas[1][0]), imag(c3->tsas[1][0])
);
fprintf(
output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
real(c3->tsas[0][1]), imag(c3->tsas[0][1]),
real(c3->tsas[1][1]), imag(c3->tsas[1][1])
);
fprintf(output, " CLUSTER\n");
pcros(vk, exri, c1, c1ao, c4);
mextc(vk, exri, c1ao->fsac, cextlr, cext);
mmulc(c1->vint, cmullr, cmul);
if (jw != 0) {
jw = 0;
// Some implicit loops writing to binary.
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
double value = cext[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
double value = c1ao->scsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1ao->scscp[i]); tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = c1ao->ecsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 2; j++) {
double value = gap[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 3; j++) {
double value = tqce[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 3; j++) {
double value = tqcs[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = real(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 3; i++) {
double value = u[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = up[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = un[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
// label 254
for (int i = 0; i < 16; i++) {
double value = real(c1->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
double value = cmul[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
int jlr = 2;
for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
int ipol = (ilr290 % 2 == 0) ? 1 : -1;
if (ilr290 == 2) jlr = 1;
double extsec = c1ao->ecsc[ilr290 - 1];
double qext = extsec * sqsfi / c3->gcs;
double extrat = extsec / c3->ecs;
double scasec = c1ao->scsc[ilr290 - 1];
double albedc = scasec / extsec;
double qsca = scasec * sqsfi / c3->gcs; double scarat = scasec / c3->scs;
double abssec = extsec - scasec;
double qabs = abssec * sqsfi / c3->gcs;
double absrat = 1.0;
double ratio = c3->acs / c3->ecs;
if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
double qschu = imag(s0) * csch;
double pschu = real(s0) * csch;
s0mag = cabs(s0) * cs0;
double refinr = real(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(c3->tfsas);
double extcor = imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 273
fprintf(output, " CIRC %2d\n", ipol);
}
// label 275
fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
fprintf(
output, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
scasec, abssec, extsec, albedc
);
fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
fprintf(
output, " %14.7lE%15.7lE%15.7lE\n",
qsca, qabs, qext
);
fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
fprintf(
output, " %14.7lE%15.7lE%15.7lE\n",
scarat, absrat, extrat
);
fprintf(
output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
ilr290, ilr290, real(c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]),
jlr, ilr290, real(c1ao->fsac[jlr - 1][ilr290 - 1]), imag(c1ao->fsac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",
ilr290, ilr290, real(c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(c1ao->sac[ilr290 - 1][ilr290 - 1]),
jlr, ilr290, real(c1ao->sac[jlr - 1][ilr290 - 1]), imag(c1ao->sac[jlr - 1][ilr290 - 1])
);
fprintf(
output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
ilr290, ilr290, refinr, ilr290, ilr290, extcor
);
fprintf(
output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
qschu, pschu, s0mag
);
bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
if (!goto190) {
gapv[0] = gap[0][ilr290 - 1];
gapv[1] = gap[1][ilr290 - 1];
gapv[2] = gap[2][ilr290 - 1];
double extins = c1ao->ecsc[ilr290 - 1];
double scatts = c1ao->scsc[ilr290 - 1];
double rapr, cosav, fp, fn, fk, fx, fy, fz;
rftr(u, up, un, gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
fprintf(output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz);
tqev[0] = tqce[ilr290 - 1][0];
tqev[1] = tqce[ilr290 - 1][1];
tqev[2] = tqce[ilr290 - 1][2];
tqsv[0] = tqcs[ilr290 - 1][0];
tqsv[1] = tqcs[ilr290 - 1][1];
tqsv[2] = tqcs[ilr290 - 1][2];
double tep, ten, tek, tsp, tsn, tsk; tqr(u, up, un, tqev, tqsv, tep, ten, tek, tsp, tsn, tsk);
fprintf(output, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n", tep, ten, tek);
fprintf(output, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n", tsp, tsn, tsk);
fprintf(
output, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
tqce[ilr290 - 1][0], tqce[ilr290 - 1][1], tqce[ilr290 - 1][2]
);
fprintf(
output, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
tqcs[ilr290 - 1][0], tqcs[ilr290 - 1][1], tqcs[ilr290 - 1][2]
);
}
} //ilr290 loop
double rbirif = (real(c1ao->fsac[0][0]) - real(c1ao->fsac[1][1])) / real(c1ao->fsac[0][0]);
double rdichr = (imag(c1ao->fsac[0][0]) - imag(c1ao->fsac[1][1])) / imag(c1ao->fsac[0][0]);
fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
fprintf(output, " MULC\n");
for (int i = 0; i < 4; i++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
);
}
fprintf(output, " MULCLR\n");
for (int i = 0; i < 4; i++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
);
}
if (iavm != 0) {
mmulc(c1ao->vintm, cmullr, cmul);
// Some implicit loops writing to binary.
for (int i = 0; i < 16; i++) {
double value;
value = real(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
value = imag(c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
double value = cmul[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
if (inpol == 0) {
fprintf(output, " LIN\n");
} else { // label 316
fprintf(output, " CIRC\n");
}
// label 318
fprintf(output, " MULC\n");
for (int i = 0; i < 4; i++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
);
}
fprintf(output, " MULCLR\n");
for (int i = 0; i < 4; i++) {
fprintf(
output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
);
}
}
// label 420, continues jphs loop if (isam < 1) phs += phsstp;
} // jphs loop, labeled 480
if (isam <= 1) thsl += thsstp;
} // jths loop, labeled 482
ph += phstp;
} // jph484 loop
th += thstp;
} // jth486 loop
// delete[] u;
// delete[] us;
// delete[] un;
// delete[] uns;
// delete[] up;
// delete[] ups;
// delete[] unmp;
// delete[] unsmp;
// delete[] upmp;
// delete[] upsmp;
logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
return jer;
}