From ef512fd28cff9941a1a984216e7ac63cea944b66 Mon Sep 17 00:00:00 2001
From: "Mulas, Giacomo" <gmulas@oa-cagliari.inaf.it>
Date: Fri, 15 Mar 2024 20:04:56 +0100
Subject: [PATCH] Move the content of the loop over scales in a separate
function, so it can be called with fresh arguments every time to check
parallelism and dependencies.
---
src/cluster/cluster.cpp | 1404 +++++++++++++++++++----------------
src/include/Configuration.h | 14 +
2 files changed, 778 insertions(+), 640 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index ce5de2c0..b53681b4 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -1,6 +1,6 @@
/* Distributed under the terms of GPLv3 or later. See COPYING for details. */
-/*! \file cluster.cpp
+/*! \file cluster.cp
*
* \brief Implementation of the calculation for a cluster of spheres.
*/
@@ -43,6 +43,8 @@
using namespace std;
+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);
+
/*! \brief C++ implementation of CLU
*
* \param config_file: `string` Name of the configuration file.
@@ -112,20 +114,20 @@ void cluster(string config_file, string data_file, string output_path) {
// End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
- int jer = 0, lcalc = 0;
- dcomplex arg = 0.0 + 0.0 * I;
+ int jer = 0;
+ // int lcalc = 0;
dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (sconf->iog_vec[ci -1] >= ci) configurations++;
}
C2 *c2 = new C2(nsph, configurations, npnt, npntts);
- np_int ndit = 2 * nsph * c4->nlim;
- 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);
- }
+ // np_int ndit = 2 * nsph * c4->nlim;
+ // 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);
+ // }
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph]();
@@ -164,16 +166,6 @@ void cluster(string config_file, string data_file, string output_path) {
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]();
@@ -189,8 +181,8 @@ void cluster(string config_file, string data_file, string output_path) {
cmullr[ci] = new double[4]();
cmul[ci] = new double[4]();
}
- int isq, ibf;
- double scan, cfmp, sfmp, cfsp, sfsp;
+ //int isq, ibf;
+ //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",
@@ -279,619 +271,14 @@ void cluster(string config_file, string data_file, string output_path) {
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
double wn = wp / 3.0e8;
- double sqsfi = 1.0;
if (sconf->idfc < 0) {
vk = sconf->xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
- printf("INFO: running scale iteration %d of %d...", jxi488, nxi);
- int jaw = 1;
- fprintf(output, "========== JXI =%3d ====================\n", jxi488);
- double xi = sconf->scale_vec[jxi488 - 1];
- double vkarg = 0.0;
- if (sconf->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");
- break; // jxi488 loop: goes to memory cleaning and return
- }
- for (int i132 = 1; i132 <= nsph; i132++) {
- int iogi = c1->iog[i132 - 1];
- if (iogi != i132) {
- for (int l123 = 1; l123 <= gconf->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 = sconf->nshl_vec[i132 - 1];
- int ici = (nsh + 1) / 2;
- if (sconf->idfc == 0) {
- for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
- } else {
- if (jxi488 == 1) {
- for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
- }
- }
- if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
- dme(
- c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,
- c1, c2, jer, lcalc, arg
- );
- if (jer != 0) {
- fprintf(output, " STOP IN DME\n");
- break;
- }
- }
- if (jer != 0) break;
- } // i132 loop
- cms(am, c1, c1ao, c4, c6);
- invert_matrix(am, ndit, jer, mxndm);
- if (jer != 0) break; // jxi488 loop: goes to memory clean
- ztm(am, c1, c1ao, c4, c6, c9);
- if (sconf->idfc >= 0) {
- if (jxi488 == gconf->jwtm) {
- int nlemt = 2 * c4->nlem;
- string ttms_name = output_path + "/c_TTMS.hd5";
- TransitionMatrix ttms(nlemt, lm, vk, exri, c1ao->am0m);
- ttms.write_binary(ttms_name, "HDF5");
- ttms_name = output_path + "/c_TTMS";
- ttms.write_binary(ttms_name);
- }
- }
- // 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 * sconf->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++) { // QUESTION: check that 16 is a fixed dimension
- c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
- } // j225 loop
- mmulc(c1ao->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(c1ao->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(c1ao->vint[i]);
- tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->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
- printf(" done.\n");
+ // 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
+ 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);
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
@@ -914,8 +301,8 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] zpv[zi];
}
delete[] zpv;
- delete[] am_vector;
- delete[] am;
+ //delete[] am_vector;
+ //delete[] am;
//delete[] tam;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
@@ -949,16 +336,6 @@ void cluster(string config_file, string data_file, string output_path) {
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;
@@ -981,3 +358,750 @@ void cluster(string config_file, string data_file, string output_path) {
delete gconf;
printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
}
+
+
+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 = sconf->number_of_scales;
+ int nsph = gconf->number_of_spheres;
+ np_int mxndm = (np_int)gconf->mxndm;
+ int inpol = gconf->in_pol;
+ int iavm = gconf->iavm;
+ int npnt = gconf->npnt;
+ int npntts = gconf->npntts;
+ int isam = gconf->meridional_type;
+ int lm = gconf->l_max;
+ double th = gconf->in_theta_start;
+ double thstp = gconf->in_theta_step;
+ double thlst = gconf->in_theta_end;
+ double ths = gconf->sc_theta_start;
+ double thsstp = gconf->sc_theta_step;
+ double thslst = gconf->sc_theta_end;
+ double ph = gconf->in_phi_start;
+ double phstp = gconf->in_phi_step;
+ double phlst = gconf->in_phi_end;
+ double phs = gconf->sc_phi_start;
+ double phsstp = gconf->sc_phi_step;
+ double phslst = gconf->sc_phi_end;
+ double th1 = th;
+ double ph1 = ph;
+ double ths1 = ths;
+ double phs1 = phs;
+ double thsca = 0.0;
+ if (isam > 1) thsca = ths - th;
+ 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 scan;
+ double cfmp;
+ double sfmp;
+ double cfsp;
+ double sfsp;
+ double sqsfi = 1.0;
+ double exri = sqrt(sconf->exdc);
+ int jer = 0;
+ int lcalc = 0;
+ dcomplex arg = 0.0 + 0.0 * I;
+ double wp = sconf->wp;
+ double wn = wp / 3.0e8;
+ double vk = 0.0;
+ if (sconf->idfc < 0) vk = sconf->xip * wn;
+ np_int ndit = 2 * nsph * c4->nlim;
+ 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);
+ }
+ int isq;
+ int ibf;
+
+
+
+ printf("INFO: running scale iteration %d of %d...", jxi488, nxi);
+ int jaw = 1;
+ fprintf(output, "========== JXI =%3d ====================\n", jxi488);
+ double xi = sconf->scale_vec[jxi488 - 1];
+ double vkarg = 0.0;
+ if (sconf->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");
+ 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; // 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 <= gconf->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 = sconf->nshl_vec[i132 - 1];
+ int ici = (nsh + 1) / 2;
+ if (sconf->idfc == 0) {
+ for (int ic = 0; ic < ici; ic++)
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+ } else {
+ if (jxi488 == 1) {
+ for (int ic = 0; ic < ici; ic++)
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
+ }
+ }
+ if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+ dme(
+ c4->li, i132, npnt, npntts, vkarg, sconf->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
+ 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);
+ if (sconf->idfc >= 0) {
+ if (jxi488 == gconf->jwtm) {
+ int nlemt = 2 * c4->nlem;
+ string ttms_name = output_path + "/c_TTMS.hd5";
+ TransitionMatrix ttms(nlemt, lm, vk, exri, c1ao->am0m);
+ ttms.write_binary(ttms_name, "HDF5");
+ ttms_name = output_path + "/c_TTMS";
+ ttms.write_binary(ttms_name);
+ }
+ }
+ // 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 * sconf->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++) { // QUESTION: check that 16 is a fixed dimension
+ c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
+ } // j225 loop
+ mmulc(c1ao->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(c1ao->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(c1ao->vint[i]);
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = imag(c1ao->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;
+ delete[] am_vector;
+ delete[] am;
+
+ printf(" done.\n");
+
+ return jer;
+
+}
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index e8e3b895..b0b8c63c 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -30,6 +30,18 @@
#ifndef INCLUDE_CONFIGURATION_H_
#define INCLUDE_CONFIGURATION_H_
+class ScattererConfiguration;
+class GeometryConfiguration;
+class C1;
+class C1_AddOns;
+class C2;
+class C3;
+class C4;
+class C6;
+class C9;
+
+#include <fstream>
+
/**
* \brief A class to represent the configuration of the scattering geometry.
*
@@ -41,6 +53,7 @@
class GeometryConfiguration {
//! Temporary work-around to allow cluster() and sphere() peeking in.
friend void cluster(std::string, std::string, std::string);
+ friend int cluster_jxi488_cycle(int, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *, C3 *, C4 *, C6 *, C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int, int, int, int, int, int, int, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *);
friend void sphere(std::string, std::string, std::string);
protected:
//! \brief Number of spherical components.
@@ -167,6 +180,7 @@ public:
class ScattererConfiguration {
//! Temporary work-around to allow cluster() and sphere() peeking in.
friend void cluster(std::string, std::string, std::string);
+ friend int cluster_jxi488_cycle(int, ScattererConfiguration *, GeometryConfiguration *, C1 *, C1_AddOns *, C2 *, C3 *, C4 *, C6 *, C9 *, std::FILE *, std::string, double *, double **, dcomplex **, double **, dcomplex **, double ****, double **, dcomplex **, int, int, int, int, int, int, int, double *, double *, double **, dcomplex **, double **, dcomplex **, double **, dcomplex **, double *, std::fstream &, double **, double **, double **, double **, double *, double *, double *);
friend void sphere(std::string, std::string, std::string);
protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
--
GitLab