diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 971d83f6139a615b210a1bdb60077fde51867485..aae0e7436184dce99559440f729eb1135ca742b9 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.
*/
@@ -8,6 +8,9 @@
#include <exception>
#include <fstream>
#include <string>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
@@ -47,6 +50,57 @@
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...
+//
+// struct clu_jxi488_cycle_args {
+// 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 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, dcomplex **am, int isq, int ibf);
+
/*! \brief C++ implementation of CLU
*
* \param config_file: `string` Name of the configuration file.
@@ -112,27 +166,17 @@ void cluster(string config_file, string data_file, string output_path) {
c1->rxx[c1i] = gconf->sph_x[c1i];
c1->ryy[c1i] = gconf->sph_y[c1i];
c1->rzz[c1i] = gconf->sph_z[c1i];
+ }
+ for (int c1i = 0; c1i < c1->configurations; c1i++) {
c1->ros[c1i] = sconf->radii_of_spheres[c1i];
}
- C4 *c4 = new C4;
- c4->li = gconf->li;
- c4->le = gconf->le;
- c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;
- c4->nv3j = (c4->lm * (c4->lm + 1) * (2 * c4->lm + 7)) / 6;
- c4->nsph = nsph;
- // The following is needed to initialize C1_AddOns
- c4->litpo = c4->li + c4->li + 1;
- c4->litpos = c4->litpo * c4->litpo;
- c4->lmtpo = c4->li + c4->le + 1;
- c4->lmtpos = c4->lmtpo * c4->lmtpo;
- c4->nlim = c4->li * (c4->li + 2);
- c4->nlem = c4->le * (c4->le + 2);
- c4->lmpo = c4->lm + 1;
+ C4 *c4 = new C4(gconf->li, gconf->le, nsph);
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, lcalc = 0;
+ int jer = 0;
+ int lcalc = 0;
dcomplex arg = 0.0 + 0.0 * I;
dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = 0;
@@ -305,629 +349,339 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
- for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
- string message = ("INFO: running scale iteration " + to_string(jxi488)
- + " of " + to_string(nxi) + "...");
- logger.log(message, LOG_INFO);
- 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
+ // 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, am, isq, ibf);
+
+ // 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
+ {
+ // 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
+ ScattererConfiguration *sconf_2 = new ScattererConfiguration(*sconf);
+ GeometryConfiguration *gconf_2 = new GeometryConfiguration(*gconf);
+ C1 *c1_2 = new C1(*c1);
+ C1_AddOns *c1ao_2 = new C1_AddOns(*c1ao);
+ C2 *c2_2 = new C2(*c2);
+ C3 *c3_2 = new C3(*c3);
+ C4 *c4_2 = new C4(*c4);
+ C6 *c6_2 = new C6(*c6);
+ C9 *c9_2 = new C9(*c9);
+ // 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();
+ FILE *output_2 = fopen((output_path + "/c_OCLU_" + to_string(myompthread)).c_str(), "w");
+ fstream tppoan_2;
+ tppoan_2.open((output_path + "/c_TPPOAN_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
+#else
+ // If OpenMP is not enabled, just keep doing the output to the global files
+ FILE *output_2 = output;
+ fstream &tppoan_2 = tppoan;
+#endif
+ double *gaps_2 = new double[nsph]();
+ double **tqse_2 = new double*[2];
+ double **tqce_2 = new double*[2];
+ double **tqcs_2 = new double*[2];
+ dcomplex **tqspe_2 = new dcomplex*[2];
+ dcomplex **tqcpe_2 = new dcomplex*[2];
+ dcomplex **tqcps_2 = new dcomplex*[2];
+ double **tqss_2 = new double*[2];
+ dcomplex **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];
+ }
}
- 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;
+ double ****zpv_2 = 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_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];
}
}
- if (jer != 0) break;
- } // i132 loop
- logger.push("DEBUG: " + TOSTRING(cms(am, c1, c1ao, c4, c6);));
- 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);
+ }
+ double **gapm_2 = new double*[3];
+ dcomplex **gappm_2 = new dcomplex*[3];
+ double **gap_2 = new double*[3];
+ dcomplex **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];
+ }
+ }
+ double *argi_2 = new double[1]();
+ argi_2[0] = argi[0];
+ double *args_2 = new double[1]();
+ args_2[0] = args[0];
+ double *duk_2 = new double[3]();
+ for (int di=0; di<3; di++) duk_2[di] = duk[di];
+ double **cextlr_2 = new double*[4];
+ double **cext_2 = new double*[4];
+ double **cmullr_2 = new double*[4];
+ double **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];
}
}
- // label 156: continue from here
- if (inpol == 0) {
- fprintf(output, " LIN\n");
- } else { // label 158
- fprintf(output, " CIRC\n");
+ double *gapv_2 = new double[3]();
+ for (int gi=0; gi<3; gi++) gapv_2[gi] = gapv[gi];
+ double *tqev_2 = new double[3]();
+ double *tqsv_2 = new double[3]();
+ for (int ti=0; ti<3; ti++) {
+ tqev_2[ti] = tqev[ti];
+ tqsv_2[ti] = tqsv[ti];
+ }
+ 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 = new double[3]();
+ double *us_2 = new double[3]();
+ double *un_2 = new double[3]();
+ double *uns_2 = new double[3]();
+ double *up_2 = new double[3]();
+ double *ups_2 = new double[3]();
+ double *unmp_2 = new double[3]();
+ double *unsmp_2 = new double[3]();
+ double *upmp_2 = new double[3]();
+ double *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];
}
- // 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]));
+ 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;
+ dcomplex **am_2 = new dcomplex*[ndit];
+ dcomplex *am_vector_2 = new dcomplex[ndit * ndit]();
+ for (int ai = 0; ai < ndit; ai++) {
+ am_2[ai] = (am_vector_2 + ai * 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;
+
+
+#pragma omp barrier
+#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, am_2, isq_2, ibf_2);
+ }
+
+#pragma omp barrier
+ 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;
+#ifdef _OPENMP
+ fclose(output_2);
+ tppoan_2.close();
+#pragma omp barrier
+ {
+ string message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
+ logger.log(message);
+ }
+#endif
+ delete[] gaps_2;
+ for (int ti = 0; ti <2 -1; 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];
}
- // 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]);
+ delete[] zpv_2[zi][zj];
}
- } // 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
- 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
- logger.log(" done.\n", LOG_INFO);
- } // jxi488 loop
+ 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;
+ delete[] am_vector_2;
+ delete[] am_2;
+ } // closes pragma omp parallel
+#ifdef _OPENMP
+#pragma omp barrier
+ {
+ for (int ri = 0; 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");
+ }
+ for (int ri = 0; ri < ompnumthreads; ri++) {
+ string partial_file_name = output_path + "/c_TPPOAN_" + to_string(ri);
+ string 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
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
logger.err("\nERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
- delete c1;
- delete c1ao;
- delete c3;
- delete c4;
- delete c6;
- delete c9;
for (int zi = c4->lm - 1; zi > -1; zi--) {
for (int zj = 2; zj > -1; zj--) {
delete[] zpv[zi][zj][1];
@@ -937,6 +691,12 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] zpv[zi];
}
delete[] zpv;
+ delete c1;
+ delete c1ao;
+ delete c3;
+ delete c4;
+ delete c6;
+ delete c9;
delete[] am_vector;
delete[] am;
//delete[] tam;
@@ -1005,3 +765,749 @@ void cluster(string config_file, string data_file, string output_path) {
logger.flush(LOG_INFO);
logger.log("Finished: output written to " + output_path + "/c_OCLU\n");
}
+
+
+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, dcomplex **am, int isq, int ibf)
+{
+
+ // 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;
+
+
+ logger.log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
+ 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;
+
+ logger.log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
+
+ return jer;
+
+}
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 3bd517d08e33b377e81b3730e743c5ceb2ca0038..381f9ec0ed651f664d3fd1fbcb77e9d339551403 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -40,8 +40,6 @@ protected:
int nsph;
//! \brief Maximum order of field expansion.
int lm;
- //! \brief NLMMT = 2 * LM * (LM + 2)
- int nlmmt;
//! \brief Contiguous RMI vector.
dcomplex *vec_rmi;
//! \brief Contiguous REI vector.
@@ -52,6 +50,10 @@ protected:
dcomplex *vec_vints;
public:
+ //! \brief NLMMT = 2 * LM * (LM + 2).
+ int nlmmt;
+ //! \brief Number of configurations
+ int configurations;
//! \brief QUESTION: definition?
dcomplex **rmi;
//! \brief QUESTION: definition?
@@ -104,6 +106,12 @@ public:
*/
C1(int ns, int l_max, int *nshl, int *iog);
+ /*! \brief C1 instance constructor copying all contents from a preexisting template
+ *
+ * \param rhs: `C1` preexisting template.
+ */
+ C1(const C1& rhs);
+
//! \brief C1 instance destroyer.
~C1();
};
@@ -112,10 +120,14 @@ public:
*
*/
class C2 {
+protected:
//! \brief Number of spheres.
int nsph;
//! \brief Number of required orders.
int nhspo;
+ //! \brief QUESTION: what is nl?
+ int nl;
+
public:
//! \brief QUESTION: definition?
dcomplex *ris;
@@ -137,6 +149,12 @@ public:
*/
C2(int ns, int nl, int npnt, int npntts);
+ /*! \brief C2 instance constructor copying its contents from preexisting instance.
+ *
+ * \param rhs: `C2` object to copy contents from
+ */
+ C2(const C2& rhs);
+
//! \brief C2 instance destroyer.
~C2();
};
@@ -162,6 +180,10 @@ public:
*/
C3();
+ /*! \brief C3 instance constructor copying its contents from a preexisting object.
+ */
+ C3(const C3& rhs);
+
/*! \brief C3 instance destroyer.
*/
~C3();
@@ -169,7 +191,8 @@ public:
/*! \brief Representation of the FORTRAN C4 blocks.
*/
-struct C4 {
+class C4 {
+public:
//! \brief QUESTION: definition?
int litpo;
//! \brief QUESTION: definition?
@@ -194,6 +217,17 @@ struct C4 {
int nsph;
//! \brief QUESTION: definition?
int nv3j;
+
+ /*! \brief C3 instance constructor.
+ */
+ C4(int li, int le, int nsph);
+ /*! \brief C3 instance constructor copying its contents from a preexisting object.
+ */
+ C4(const C4& rhs);
+
+ /*! \brief C3 instance destroyer.
+ */
+ ~C4();
};
/*! \brief Vectors and matrices that are specific to cluster C1 blocks.
@@ -207,6 +241,18 @@ protected:
int nlemt;
//! \brief Maximum expansion order plus one. QUESTION: correct?
int lmpo;
+ //! \brief QUESTION: definition?
+ int litpo;
+ //! \brief QUESTION: definition?
+ int lmtpo;
+ //! \brief QUESTION: definition?
+ int litpos;
+ //! \brief QUESTION: definition?
+ int lmtpos;
+ //! \brief QUESTION: definition?
+ int nv3j;
+ //! \brief QUESTION: definition?
+ int lm;
public:
//! \brief QUESTION: definition?
@@ -224,6 +270,10 @@ public:
//! \brief QUESTION: definition?
dcomplex *vintm;
//! \brief QUESTION: definition?
+ dcomplex *vint;
+ //! \brief QUESTION: definition?
+ dcomplex *vintm;
+ //! \brief QUESTION: definition?
dcomplex *vintt;
//! \brief QUESTION: definition?
dcomplex **fsac;
@@ -260,6 +310,12 @@ public:
*/
C1_AddOns(C4 *c4);
+ /*! \brief C1_AddOns instance constructor copying contents from a preexisting object.
+ *
+ * \param rhs: `C1_AddOns` preexisting object to copy from.
+ */
+ C1_AddOns(const C1_AddOns& rhs);
+
//! \brief C1_AddOns instance destroyer.
~C1_AddOns();
};
@@ -268,6 +324,8 @@ public:
*/
class C6 {
public:
+ //! \brief LMTPO = 2 * LM + 1.
+ int lmtpo;
//! \brief QUESTION: definition?
double *rac3j;
@@ -277,6 +335,12 @@ public:
*/
C6(int lmtpo);
+ /*! \brief C6 instance constructor copying contents from preexisting object.
+ *
+ * \param lmtpo: `int` QUESTION: definition?
+ */
+ C6(const C6& rhs);
+
/*! \brief C6 instance destroyer.
*/
~C6();
@@ -290,6 +354,11 @@ protected:
int gis_size_0;
//! \brief Number of rows in the SAM matrix
int sam_size_0;
+ //! \brief QUESTION: definition?
+ int nlem;
+ //! \brief QUESTION: definition?
+ int nlemt;
+
public:
//! \brief QUESTION: definition?
dcomplex **gis;
@@ -307,6 +376,12 @@ public:
*/
C9(int ndi, int nlem, int ndit, int nlemt);
+ /*! \brief C9 instance constructor copying contents from preexisting object.
+ *
+ * \param rhs: `C9` preexisting object to copy from
+ */
+ C9(const C9& rhs);
+
/*! \brief C9 instance destroyer.
*/
~C9();
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index e8e3b8951b4b01538cbf4f541b13b100d84e3e94..591fc0bad9d6b428e579896b447764a556115624 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 *, int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex **, int, int);
friend void sphere(std::string, std::string, std::string);
protected:
//! \brief Number of spherical components.
@@ -139,6 +152,11 @@ public:
double sc_ph_start, double sc_ph_step, double sc_ph_end,
int jwtm
);
+ /*! \brief Build a scattering geometry configuration structure copying it from an existing one.
+ *
+ * \param rhs: `GeometryConfiguration` preexisting object to copy from.
+ */
+ GeometryConfiguration(const GeometryConfiguration& rhs);
/*! \brief Destroy a GeometryConfiguration instance.
*/
@@ -167,6 +185,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 *, int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex **, int, int);
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].
@@ -185,9 +204,11 @@ protected:
std::string reference_variable_name;
//! \brief Number of spherical components.
int number_of_spheres;
+ //! \brief Number of configurations.
+ int configurations;
//! \brief Number of scales to use in calculation.
int number_of_scales;
- //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 contants).
+ //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 constants).
int idfc;
//! \brief External medium dielectric constant. QUESTION: correct?
double exdc;
@@ -268,6 +289,7 @@ public:
*/
ScattererConfiguration(
int nsph,
+ int configs,
double *scale_vector,
int nxi,
std::string variable_name,
@@ -282,6 +304,29 @@ public:
double wp,
double xip
);
+ /*! \brief Build a scatterer configuration structure copying its contents from a preexisting one.
+ *
+ * Prepare a default configuration structure by allocating the necessary
+ * memory structures.
+ *
+ * \param nsph: `int` The number of spheres in the simulation.
+ * \param scale_vector: `double*` The radiation-particle scale vector.
+ * \param nxi: `int` The number of radiation-particle scalings.
+ * \param variable_name: `string` The name of the radiation-particle scaling type.
+ * \param iog_vector: `int*` Array of sphere identification numbers. QUESTION: correct?
+ * \param ros_vector: `double*` Sphere radius array.
+ * \param nshl_vector: `int*` Array of layer numbers.
+ * \param rcf_vector: `double**` Array of fractional break radii. QUESTION: correct?
+ * \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,
+ * \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor
+ * for dimensions).
+ * \param dc_matrix: `complex double ***` Matrix of reference dielectric constants.
+ * \param has_external: `bool` Flag to set whether to add an external spherical layer.
+ * \param exdc: `double` External medium dielectric constant.
+ * \param wp: `double` wp
+ * \param xip: `double` xip
+ */
+ ScattererConfiguration(const ScattererConfiguration& rhs);
/*! \brief Destroy a scatterer configuration instance.
*/
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index e93af0c274c345764b45ad29f2befc13386a89b2..b76e3cc1baab2db561d86689a920b07fd21e0103 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -36,15 +36,15 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
vec_w = new dcomplex[4 * nlmmt]();
w = new dcomplex*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
- int configurations = 0;
+ configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
vint = new dcomplex[16]();
vec_vints = new dcomplex[nsph * 16]();
vints = new dcomplex*[nsph];
- rc = new double*[configurations];
- nshl = new int[configurations]();
+ rc = new double*[nsph];
+ nshl = new int[nsph]();
iog = new int[nsph]();
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
@@ -77,6 +77,84 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
}
}
+C1::C1(const C1& rhs) {
+ nlmmt = rhs.nlmmt;
+ nsph = rhs.nsph;
+ lm = rhs.lm;
+
+ vec_rmi = new dcomplex[nsph * lm]();
+ vec_rei = new dcomplex[nsph * lm]();
+ rmi = new dcomplex*[lm];
+ rei = new dcomplex*[lm];
+ for (int ri = 0; ri < lm; ri++) {
+ rmi[ri] = &(vec_rmi[nsph * ri]);
+ rei[ri] = &(vec_rei[nsph * ri]);
+ /*! Copy the contents from the template */
+ for (int rj=0; rj<nsph; rj++) {
+ rmi[ri][rj] = rhs.rmi[ri][rj];
+ rei[ri][rj] = rhs.rei[ri][rj];
+ }
+ }
+ vec_w = new dcomplex[4 * nlmmt]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) {
+ w[wi] = &(vec_w[4 * wi]);
+ for (int wj=0; wj<4; wj++) w[wi][wj] = rhs.w[wi][wj];
+ }
+ configurations = rhs.configurations;
+ vec_vints = new dcomplex[nsph * 16]();
+ vints = new dcomplex*[nsph];
+ rc = new double*[nsph];
+ nshl = new int[nsph]();
+ for (int ni = 0; ni < nsph; ni++) nshl[ni] = rhs.nshl[ni];
+ iog = new int[nsph]();
+ for (int ni = 0; ni < nsph; ni++) iog[ni] = rhs.iog[ni];
+ int conf_index = 0;
+ for (int vi = 0; vi < nsph; vi++) {
+ vints[vi] = &(vec_vints[16 * vi]);
+ for (int vj = 0; vj < 16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
+ }
+ for (int ri=0; ri<configurations; ri++) {
+ rc[ri] = new double[nshl[ri]]();
+ for (int rj=0; rj<nshl[ri]; rj++) rc[ri][rj] = rhs.rc[ri][rj];
+ }
+ fsas = new dcomplex[nsph]();
+ sscs = new double[nsph]();
+ sexs = new double[nsph]();
+ sabs = new double[nsph]();
+ sqscs = new double[nsph]();
+ sqexs = new double[nsph]();
+ sqabs = new double[nsph]();
+ gcsv = new double[nsph]();
+ rxx = new double[nsph]();
+ ryy = new double[nsph]();
+ rzz = new double[nsph]();
+ ros = new double[nsph]();
+
+ sas = new dcomplex**[nsph];
+ for (int si = 0; si < nsph; si++) {
+ sas[si] = new dcomplex*[2];
+ for (int sj=0; sj<2; sj++) {
+ sas[si][sj] = new dcomplex[2]();
+ for (int sk=0; sk<2; sk++) sas[si][sj][sk] = rhs.sas[si][sj][sk];
+ }
+ fsas[si] = rhs.fsas[si];
+ sscs[si] = rhs.sscs[si];
+ sexs[si] = rhs.sexs[si];
+ sabs[si] = rhs.sabs[si];
+ sqscs[si] = rhs.sqscs[si];
+ sqexs[si] = rhs.sqexs[si];
+ sqabs[si] = rhs.sqabs[si];
+ gcsv[si] = rhs.gcsv[si];
+ rxx[si] = rhs.rxx[si];
+ ryy[si] = rhs.ryy[si];
+ rzz[si] = rhs.rzz[si];
+ }
+ for (int si = 0; si < configurations; si++) {
+ ros[si] = rhs.ros[si];
+ }
+}
+
C1::~C1() {
delete[] vec_rmi;
delete[] vec_rei;
@@ -85,11 +163,8 @@ C1::~C1() {
delete[] vec_w;
delete[] w;
int conf_index = 0;
- for (int ci = 1; ci <= nsph; ci++) {
- if (iog[ci] >= ci) {
- delete[] rc[conf_index];
- conf_index++;
- }
+ for (int ci = 0; ci < configurations; ci++) {
+ delete[] rc[ci];
}
delete[] rc;
delete[] vint;
@@ -121,15 +196,22 @@ C1_AddOns::C1_AddOns(C4 *c4) {
nsph = c4->nsph;
lmpo = c4->lmpo;
nlemt = 2 * c4->nlem;
- vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();
- vj0 = new dcomplex[nsph * c4->lmtpo]();
- vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
- vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();
- vyj0 = new dcomplex[nsph * c4->lmtpos]();
+ litpo = c4->litpo;
+ lmtpo = c4->lmtpo;
+ litpos = c4->litpos;
+ lmtpos = c4->lmtpos;
+ nv3j = c4->nv3j;
+ lm = c4->lm;
+ vh = new dcomplex[(nsph * nsph - 1) * litpo]();
+ vj0 = new dcomplex[nsph * lmtpo]();
+ vj = new dcomplex[1]();
+ vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
+ vyj0 = new dcomplex[nsph * lmtpos]();
am0m = new dcomplex*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
am0m[ai] = new dcomplex[nlemt]();
}
+ vint = new dcomplex[16]();
vintm = new dcomplex[16]();
vintt = new dcomplex[16]();
// This calculates the size of v3j0
@@ -149,11 +231,93 @@ C1_AddOns::C1_AddOns(C4 *c4) {
ecscp = new dcomplex[2]();
scscpm = new dcomplex[2]();
ecscpm = new dcomplex[2]();
+ v3j0 = new double[nv3j]();
+ ind3j = new int*[lmpo];
+ for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[lm]();
+ sscs = new double[nsph]();
+ ecscm = new double[2]();
+ scscm = new double[2]();
+ scsc = new double[2]();
+ ecsc = new double[2]();
+}
+
+C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
+ nsph = rhs.nsph;
+ lmpo = rhs.lmpo;
+ nlemt = rhs.nlemt;
+ litpo = rhs.litpo;
+ lmtpo = rhs.lmtpo;
+ litpos = rhs.litpos;
+ lmtpos = rhs.lmtpos;
+ nv3j = rhs.nv3j;
+ lm = rhs.lm;
+ int vhsize=(nsph * nsph - 1) * litpo;
+ vh = new dcomplex[vhsize]();
+ for (int hi=0; hi<vhsize; hi++) vh[hi] = rhs.vh[hi];
+ int vj0size=nsph * lmtpo;
+ vj0 = new dcomplex[vj0size]();
+ for (int hi=0; hi<vj0size; hi++) vj0[hi] = rhs.vj0[hi];
+ vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
+ vj[0] = rhs.vj[0];
+ int vyhjsize=(nsph * nsph - 1) * litpos;
+ vyhj = new dcomplex[vyhjsize]();
+ for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi];
+ int vyj0size = nsph * lmtpos;
+ vyj0 = new dcomplex[vyj0size]();
+ for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
+ am0m = new dcomplex*[nlemt];
+ for (int ai = 0; ai < nlemt; ai++) {
+ am0m[ai] = new dcomplex[nlemt]();
+ for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];
+ }
+ vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
+ vintm = new dcomplex[16]();
+ vintt = new dcomplex[16]();
+ for (int hi = 0; hi < 16; hi++) {
+ vint[hi] = rhs.vint[hi];
+ vintm[hi] = rhs.vintm[hi];
+ vintt[hi] = rhs.vintt[hi];
+ }
+ fsac = new dcomplex*[2];
+ sac = new dcomplex*[2];
+ fsacm = new dcomplex*[2];
+ for (int fi = 0; fi < 2; fi++) {
+ fsac[fi] = new dcomplex[2]();
+ sac[fi] = new dcomplex[2]();
+ fsacm[fi] = new dcomplex[2]();
+ for (int fj = 0; fj < 2; fj++) {
+ fsac[fi][fj] = rhs.fsac[fi][fj];
+ sac[fi][fj] = rhs.sac[fi][fj];
+ fsacm[fi][fj] = rhs.fsacm[fi][fj];
+ }
+ }
+ scscp = new dcomplex[2]();
+ ecscp = new dcomplex[2]();
+ scscpm = new dcomplex[2]();
+ ecscpm = new dcomplex[2]();
+ v3j0 = new double[nv3j]();
+ for (int hi = 0; hi < nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
+ ind3j = new int*[lmpo];
+ for (int ii = 0; ii < lmpo; ii++) {
+ ind3j[ii] = new int[lm]();
+ for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];
+ }
sscs = new double[nsph]();
+ for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];
ecscm = new double[2]();
scscm = new double[2]();
scsc = new double[2]();
ecsc = new double[2]();
+ for (int hi=0; hi<2; hi++) {
+ scscp[hi] = rhs.scscp[hi];
+ ecscp[hi] = rhs.ecscp[hi];
+ scscpm[hi] = rhs.scscpm[hi];
+ ecscpm[hi] = rhs.ecscpm[hi];
+ ecscm[hi] = rhs.ecscm[hi];
+ scscm[hi] = rhs.scscm[hi];
+ scsc[hi] = rhs.scsc[hi];
+ ecsc[hi] = rhs.ecsc[hi];
+ }
}
C1_AddOns::~C1_AddOns() {
@@ -168,7 +332,8 @@ C1_AddOns::~C1_AddOns() {
for (int ai = nlemt - 1; ai > -1; ai--) {
delete[] am0m[ai];
}
- delete am0m;
+ delete[] am0m;
+ delete[] vint;
delete[] vintm;
delete[] vintt;
for (int fi = 1; fi > -1; fi--) {
@@ -189,10 +354,11 @@ C1_AddOns::~C1_AddOns() {
delete[] ecsc;
}
-C2::C2(int ns, int nl, int npnt, int npntts) {
+C2::C2(int ns, int _nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
+ nl = _nl;
ris = new dcomplex[nhspo]();
dlri = new dcomplex[nhspo]();
vkt = new dcomplex[nsph]();
@@ -200,6 +366,26 @@ C2::C2(int ns, int nl, int npnt, int npntts) {
vsz = new double[nsph]();
}
+C2::C2(const C2& rhs) {
+ nsph = rhs.nsph;
+ nhspo = rhs.nhspo;
+ nl = rhs.nl;
+ ris = new dcomplex[nhspo]();
+ dlri = new dcomplex[nhspo]();
+ for (int ind=0; ind<nhspo; ind++) {
+ ris[ind] = rhs.ris[ind];
+ dlri[ind] = rhs.dlri[ind];
+ }
+ vkt = new dcomplex[nsph]();
+ vsz = new double[nsph]();
+ for (int ind=0; ind<nsph; ind++) {
+ vkt[ind] = rhs.vkt[ind];
+ vsz[ind] = rhs.vsz[ind];
+ }
+ dc0 = new dcomplex[nl]();
+ for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];
+}
+
C2::~C2() {
delete[] ris;
delete[] dlri;
@@ -219,21 +405,80 @@ C3::C3() {
acs = 0.0;
}
+C3::C3(const C3& rhs) {
+ tsas = new dcomplex*[2];
+ for (int ti=0; ti<2; ti++) {
+ tsas[ti] = new dcomplex[2];
+ for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];
+ }
+ tfsas = rhs.tfsas;
+ gcs = rhs.gcs;
+ scs = rhs.scs;
+ ecs = rhs.ecs;
+ acs = rhs.acs;
+}
+
C3::~C3() {
delete[] tsas[1];
delete[] tsas[0];
delete[] tsas;
}
-C6::C6(int lmtpo) {
+C4::C4(int _li, int _le, int _nsph) {
+ li = _li;
+ le = _le;
+ lm = (li > le) ? li : le;
+ nv3j = (lm * (lm + 1) * (2 * lm + 7)) / 6;
+ nsph = _nsph;
+ // The following is needed to initialize C1_AddOns
+ litpo = li + li + 1;
+ litpos = litpo * litpo;
+ lmtpo = li + le + 1;
+ lmtpos = lmtpo * lmtpo;
+ nlim = li * (li + 2);
+ nlem = le * (le + 2);
+ lmpo = lm + 1;
+}
+
+C4::C4(const C4& rhs) {
+ li = rhs.li;
+ le = rhs.le;
+ lm = rhs.lm;
+ nv3j = rhs.nv3j;
+ nsph = rhs.nsph;
+ // The following is needed to initialize C1_AddOns
+ litpo = rhs.litpo;
+ litpos = rhs.litpos;
+ lmtpo = rhs.lmtpo;
+ lmtpos = rhs.lmtpos;
+ nlim = rhs.nlim;
+ nlem = rhs.nlem;
+ lmpo = rhs.lmpo;
+}
+
+C4::~C4() {
+}
+
+C6::C6(int _lmtpo) {
+ lmtpo = _lmtpo;
rac3j = new double[lmtpo]();
}
+C6::C6(const C6& rhs) {
+ lmtpo = rhs.lmtpo;
+ rac3j = new double[lmtpo]();
+ for (int ri = 0; ri < lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
+}
+
C6::~C6() {
delete[] rac3j;
}
-C9::C9(int ndi, int nlem, int ndit, int nlemt) {
+C9::C9(int ndi, int _nlem, int ndit, int _nlemt) {
+ gis_size_0 = ndi;
+ sam_size_0 = ndit;
+ nlem = _nlem;
+ nlemt = _nlemt;
gis = new dcomplex*[ndi];
gls = new dcomplex*[ndi];
for (int gi = 0; gi < ndi; gi++) {
@@ -242,8 +487,28 @@ C9::C9(int ndi, int nlem, int ndit, int nlemt) {
}
sam = new dcomplex*[ndit];
for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();
- gis_size_0 = ndi;
- sam_size_0 = ndit;
+}
+
+C9::C9(const C9& rhs) {
+ gis_size_0 = rhs.gis_size_0;
+ sam_size_0 = rhs.sam_size_0;
+ nlem = rhs.nlem;
+ nlemt = rhs.nlemt;
+ gis = new dcomplex*[gis_size_0];
+ gls = new dcomplex*[gis_size_0];
+ for (int gi = 0; gi < gis_size_0; gi++) {
+ gis[gi] = new dcomplex[nlem]();
+ gls[gi] = new dcomplex[nlem]();
+ for (int gj=0; gj<nlem; gj++) {
+ gis[gi][gj] = rhs.gis[gi][gj];
+ gls[gi][gj] = rhs.gls[gi][gj];
+ }
+ }
+ sam = new dcomplex*[sam_size_0];
+ for (int si = 0; si < sam_size_0; si++) {
+ sam[si] = new dcomplex[nlemt]();
+ for (int sj=0; sj<nlemt; sj++) sam[si][sj] = rhs.sam[si][sj];
+ }
}
C9::~C9() {
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 1bf78b7db3e571858bee6ad893f0c575dbb33644..fa91d648c506b59a474105e039ec95e7114dfbbd 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -76,6 +76,40 @@ GeometryConfiguration::GeometryConfiguration(
sph_y = y;
sph_z = z;
}
+GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)
+{
+ number_of_spheres = rhs.number_of_spheres;
+ l_max = rhs.l_max;
+ in_pol = rhs.in_pol;
+ npnt = rhs.npnt;
+ npntts = rhs.npntts;
+ meridional_type = rhs.meridional_type;
+ li = rhs.li;
+ le = rhs.le;
+ mxndm = rhs.mxndm;
+ iavm = rhs.iavm;
+ in_theta_start = rhs.in_theta_start;
+ in_theta_step = rhs.in_theta_step;
+ in_theta_end = rhs.in_theta_end;
+ in_phi_start = rhs.in_phi_start;
+ in_phi_step = rhs.in_phi_step ;
+ in_phi_end = rhs.in_phi_end;
+ sc_theta_start = rhs.sc_theta_start;
+ sc_theta_step = rhs.sc_theta_step;
+ sc_theta_end = rhs.sc_theta_end;
+ sc_phi_start = rhs.sc_phi_start;
+ sc_phi_step = rhs.sc_phi_step;
+ sc_phi_end = rhs.sc_phi_end;
+ jwtm = rhs.jwtm;
+ sph_x = new double[number_of_spheres]();
+ sph_y = new double[number_of_spheres]();
+ sph_z = new double[number_of_spheres]();
+ for (int ni=0; ni<number_of_spheres; ni++) {
+ sph_x[ni] = rhs.sph_x[ni];
+ sph_y[ni] = rhs.sph_y[ni];
+ sph_z[ni] = rhs.sph_z[ni];
+ }
+}
GeometryConfiguration::~GeometryConfiguration() {
delete[] sph_x;
@@ -199,6 +233,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
ScattererConfiguration::ScattererConfiguration(
int nsph,
+ int configs,
double *scale_vector,
int nxi,
string variable_name,
@@ -214,6 +249,7 @@ ScattererConfiguration::ScattererConfiguration(
double x
) {
number_of_spheres = nsph;
+ configurations = configs;
number_of_scales = nxi;
reference_variable_name = variable_name;
iog_vec = iog_vector;
@@ -241,6 +277,40 @@ ScattererConfiguration::ScattererConfiguration(
}
}
}
+ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs)
+{
+ number_of_spheres = rhs.number_of_spheres;
+ configurations = rhs.configurations;
+ number_of_scales = rhs.number_of_scales;
+ reference_variable_name = rhs.reference_variable_name;
+ idfc = rhs.idfc;
+ use_external_sphere = rhs.use_external_sphere;
+ exdc = rhs.exdc;
+ wp = rhs.wp;
+ xip = rhs.xip;
+ iog_vec = new int[number_of_spheres]();
+ radii_of_spheres = new double[number_of_spheres]();
+ nshl_vec = new int[number_of_spheres]();
+ rcf = new double*[number_of_spheres];
+ scale_vec = new double[number_of_scales]();
+ dc0_matrix = new dcomplex**[configurations];
+ for (int si = 0; si < number_of_scales; si++) scale_vec[si] = rhs.scale_vec[si];
+ for (int si=0; si<number_of_spheres; si++) {
+ iog_vec[si] = rhs.iog_vec[si];
+ }
+ int dim3 = (idfc == 0) ? number_of_scales : 1;
+ for (int si=0; si<configurations; si++) {
+ radii_of_spheres[si] = rhs.radii_of_spheres[si];
+ nshl_vec[si] = rhs.nshl_vec[si];
+ rcf[si] = new double[nshl_vec[si]]();
+ dc0_matrix[si] = new dcomplex*[number_of_spheres];
+ for (int sj=0; sj<nshl_vec[si]; sj++) rcf[si][sj] = rhs.rcf[si][sj];
+ for (int sj=0; sj<number_of_spheres; sj++) {
+ dc0_matrix[si][sj] = new dcomplex[dim3]();
+ for (int sk=0; sk<dim3; sk++) dc0_matrix[si][sj][sk] = rhs.dc0_matrix[si][sj][sk];
+ }
+ }
+}
ScattererConfiguration::~ScattererConfiguration() {
int configurations = 0;
@@ -453,9 +523,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
}
}
int index = 0;
- double *ros_vector = new double[configurations]();
- int *nshl_vector = new int[configurations]();
- double **rcf_vector = new double*[configurations];
+ double *ros_vector = new double[nsph]();
+ int *nshl_vector = new int[nsph]();
+ double **rcf_vector = new double*[nsph];
for (int i113 = 1; i113 <= nsph; i113++) {
if (iog_vector[i113 - 1] < i113) continue;
str_target = file_lines[++last_read_line];
@@ -519,6 +589,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
ScattererConfiguration *config = new ScattererConfiguration(
nsph,
+ configurations,
variable_vector,
nxi,
variable_name,
@@ -568,9 +639,9 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
if (iog[ci] < ci + 1) continue;
configuration_count++;
}
- nshl_vector = new int[configuration_count]();
- ros_vector = new double[configuration_count]();
- rcf_vector = new double*[configuration_count];
+ nshl_vector = new int[nsph]();
+ ros_vector = new double[nsph]();
+ rcf_vector = new double*[nsph];
for (int i115 = 1; i115 <= nsph; i115++) {
if (iog[i115 - 1] < i115) continue;
str_name = "NSHL_" + to_string(iog[i115 - 1]);
@@ -614,6 +685,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
status = hdf_file->close();
conf = new ScattererConfiguration(
nsph,
+ configuration_count,
xi_vec,
nxi,
"XIV",
@@ -664,9 +736,9 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
for (int i = 1; i <= _nsph; i++) {
if (_iog_vec[i - 1] >= i) configurations++;
}
- _nshl_vec = new int[configurations]();
- _ros_vec = new double[configurations]();
- _rcf_vec = new double*[configurations];
+ _nshl_vec = new int[_nsph]();
+ _ros_vec = new double[_nsph]();
+ _rcf_vec = new double*[_nsph];
for (int i115 = 1; i115 <= _nsph; i115++) {
if (_iog_vec[i115 - 1] < i115) continue;
input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));
@@ -702,6 +774,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
ScattererConfiguration *conf = new ScattererConfiguration(
_nsph,
+ configurations,
_xi_vec,
_nxi,
"XIV",
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index ea5bd6fa53b80720b1531c04fff4c9a1f7f85100..2473d63602e9d589c412f8ffba01425ff0e6c6cb 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -834,8 +834,8 @@ void hjv(
const int rfj_size = (lit > lmt) ? lit : lmt;
const int rfn_size = c4->litpo;
double *rfj, *rfn;
- rfj = new double[rfj_size]();
- rfn = new double[rfn_size]();
+ rfj = new double[rfj_size+1]();
+ rfn = new double[rfn_size+1]();
jer = 0;
int ivhb = 0;
for (int nf40 = 1; nf40 <= nsphmo; nf40++) { // GPU portable?
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 72f8cb1024912d78193dedbb3d71467e8e35de53..f3a7d54d569085df555ad3443440ebb3cc1b840c 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -668,7 +668,7 @@ void rbf(int n, double x, int &nm, double sj[]) {
double abs_sb = (sb < 0.0) ? -sb : sb;
if (abs_sa > abs_sb) cs = sa / f;
else cs = sb / f0;
- for (int k = 0; k <= nm; k++) {
+ for (int k = 0; k <= nm ; k++) {
sj[k] = cs * sj[k];
}
}
diff --git a/src/make.inc b/src/make.inc
index 1d7b99e0bbeb6052302f009f7d3b2f3f57547383..e854d59e55315a2a57c46c8e3c4e58d5d8678e54 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -74,6 +74,9 @@ endif
# CXXFLAGS defines the default compilation options for the C++ compiler
ifndef CXXFLAGS
override CXXFLAGS=-O3 -ggdb -pg -coverage -I$(HDF5_INCLUDE)
+ifdef USE_OPENMP
+override CXXFLAGS+= -fopenmp
+endif
ifdef USE_LAPACK
override CXXFLAGS+= -DUSE_LAPACK -DLAPACK_ILP64
ifdef USE_MKL
diff --git a/src/scripts/pycompare.py b/src/scripts/pycompare.py
index 2add23e6d0964ba5cfa671ca2d16088f6b549bca..0ebeff4eaf0f19458a2a8478d9ba7e67ccf40532 100755
--- a/src/scripts/pycompare.py
+++ b/src/scripts/pycompare.py
@@ -141,6 +141,9 @@ def compare_files(config):
print("ERROR: {0:s} and {1:s} have different numbers of lines!".format(
config['fortran_file_name'], config['c_file_name']
))
+ if (config['log_html']):
+ print("Different file sizes. No log produced.")
+ config['log_html'] = False
return mismatch_count
## \brief Perform the comparison of two file lines.