diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index daa0288cceafde941d5f3d8dd65f11465f2adbb5..e9f666b813bbb4cf952246b423f926f3af2cc8b0 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -53,7 +53,7 @@ using namespace std;
// I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, const string& output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, 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, 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, dcomplex arg, double wn, double vk, Logger *logger);
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, FILE *output, const string& output_path, fstream& tppoan, Logger *logger);
/*! \brief C++ implementation of CLU
*
@@ -99,88 +99,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
// Shortcuts to variables stored in configuration objects
ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
double wp = sconf->wp;
- // Global variables for CLU
- C1 *c1 = new C1(gconf, sconf);
- C3 *c3 = new C3();
- C4 *c4 = new C4(gconf);
- C1_AddOns *c1ao = new C1_AddOns(c4);
- // End of add-ons initialization
- C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
- dcomplex arg = 0.0 + 0.0 * I;
- dcomplex ccsam = 0.0 + 0.0 * I;
- C2 *c2 = new C2(gconf, sconf);
- const int ndi = c4->nsph * c4->nlim;
+ ClusterIterationData *cid = new ClusterIterationData(gconf, sconf);
+ const int ndi = cid->c4->nsph * cid->c4->nlim;
np_int ndit = 2 * ndi;
logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
- C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
- double *gaps = new double[gconf->number_of_spheres]();
- double *tqev = new double[3]();
- double *tqsv = new double[3]();
- double **tqse, **tqss, **tqce, **tqcs;
- dcomplex **tqspe, **tqsps, **tqcpe, **tqcps;
- tqse = new double*[2];
- tqspe = new dcomplex*[2];
- tqss = new double*[2];
- tqsps = new dcomplex*[2];
- tqce = new double*[2];
- tqcpe = new dcomplex*[2];
- tqcs = new double*[2];
- tqcps = new dcomplex*[2];
- for (int ti = 0; ti < 2; ti++) {
- tqse[ti] = new double[gconf->number_of_spheres]();
- tqspe[ti] = new dcomplex[gconf->number_of_spheres]();
- tqss[ti] = new double[gconf->number_of_spheres]();
- tqsps[ti] = new dcomplex[gconf->number_of_spheres]();
- tqce[ti] = new double[3]();
- tqcpe[ti] = new dcomplex[3]();
- tqcs[ti] = new double[3]();
- tqcps[ti] = new dcomplex[3]();
- }
- double *gapv = new double[3]();
- dcomplex **gapp, **gappm;
- double **gap, **gapm;
- gapp = new dcomplex*[3];
- gappm = new dcomplex*[3];
- gap = new double*[3];
- gapm = new double*[3];
- for (int gi = 0; gi < 3; gi++) {
- gapp[gi] = new dcomplex[2]();
- gappm[gi] = new dcomplex[2]();
- gap[gi] = new double[2]();
- gapm[gi] = new double[2]();
- }
- double *u = new double[3]();
- double *us = new double[3]();
- double *un = new double[3]();
- double *uns = new double[3]();
- double *up = new double[3]();
- double *ups = new double[3]();
- double *unmp = new double[3]();
- double *unsmp = new double[3]();
- double *upmp = new double[3]();
- double *upsmp = new double[3]();
- double *argi = new double[1]();
- double *args = new double[1]();
- double *duk = new double[3]();
- double **cextlr, **cext;
- double **cmullr, **cmul;
- cextlr = new double*[4];
- cext = new double*[4];
- cmullr = new double*[4];;
- cmul = new double*[4];
- for (int ci = 0; ci < 4; ci++) {
- cextlr[ci] = new double[4]();
- cext[ci] = new double[4]();
- cmullr[ci] = new double[4]();
- cmul[ci] = new double[4]();
- }
- double scan, cfmp, sfmp, cfsp, sfsp;
- // End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
- nsph, c4->li, c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
+ nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
gconf->npntts, gconf->iavm, gconf->iavm
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
@@ -210,21 +137,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
fprintf(output, " \n");
- str(sconf, c1, c1ao, c3, c4, c6);
- double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
- for (int zi = 0; zi < c4->lm; zi++) {
- zpv[zi] = new double**[3];
- for (int zj = 0; zj < 3; zj++) {
- zpv[zi][zj] = new double*[2];
- for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2]();
- }
- }
- }
- thdps(c4->lm, zpv);
+ str(sconf, cid->c1, cid->c1ao, cid->c3, cid->c4, cid->c6);
+ thdps(cid->c4->lm, cid->zpv);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
- double vk = 0.0;
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream *tppoanp = new fstream;
fstream &tppoan = *tppoanp;
@@ -252,18 +168,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
- double wn = wp / 3.0e8;
- double xip = sconf->xip;
- double sqsfi = 1.0;
if (sconf->idfc < 0) {
- vk = xip * wn;
- fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
+ cid->vk = cid->xip * cid->wn;
+ fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);
fprintf(output, " \n");
}
// do the first iteration on jxi488 separately, since it seems to be different from the others
int jxi488 = 1;
chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
- int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, arg, wn, vk, logger);
+ int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, output, output_path, tppoan, logger);
chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
elapsed = end_iter_1 - start_iter_1;
message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
@@ -283,226 +196,20 @@ void cluster(const string& config_file, const string& data_file, const string& o
if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
- C1 *c1_2 = NULL;
- C1_AddOns *c1ao_2 = NULL;
- C2 *c2_2 = NULL;
- C3 *c3_2 = NULL;
- C4 *c4_2 = NULL;
- C6 *c6_2 = NULL;
- C9 *c9_2 = NULL;
+ ClusterIterationData *cid_2 = NULL;
FILE *output_2 = NULL;
fstream *tppoanp_2 = NULL;
- double *gaps_2 = NULL;
- double **tqse_2 = NULL;
- double **tqce_2 = NULL;
- double **tqcs_2 = NULL;
- dcomplex **tqspe_2 = NULL;
- dcomplex **tqcpe_2 = NULL;
- dcomplex **tqcps_2 = NULL;
- double **tqss_2 = NULL;
- dcomplex **tqsps_2 = NULL;
- double ****zpv_2 = NULL;
- double **gapm_2 = NULL;
- dcomplex **gappm_2 = NULL;
- double **gap_2 = NULL;
- dcomplex **gapp_2 = NULL;
- double *argi_2 = NULL;
- double *args_2 = NULL;
- double *duk_2 = NULL;
- double **cextlr_2 = NULL;
- double **cext_2 = NULL;
- double **cmullr_2 = NULL;
- double **cmul_2 = NULL;
- double *gapv_2 = NULL;
- double *tqev_2 = NULL;
- double *tqsv_2 = NULL;
- double *u_2 = NULL;
- double *us_2 = NULL;
- double *un_2 = NULL;
- double *uns_2 = NULL;
- double *up_2 = NULL;
- double *ups_2 = NULL;
- double *unmp_2 = NULL;
- double *unsmp_2 = NULL;
- double *upmp_2 = NULL;
- double *upsmp_2 = NULL;
- double scan_2 = scan;
- double cfmp_2 = cfmp;
- double sfmp_2 = sfmp;
- double cfsp_2 = cfsp;
- double sfsp_2 = sfsp;
- double sqsfi_2 = sqsfi;
- dcomplex arg_2 = arg;
- double wn_2 = wn;
- double vk_2 = vk;
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
if (myompthread == 0) {
- c1_2 = c1;
- c1ao_2 = c1ao;
- c2_2 = c2;
- c3_2 = c3;
- c4_2 = c4;
- c6_2 = c6;
- c9_2 = c9;
+ cid_2 = cid;
output_2 = output;
tppoanp_2 = tppoanp;
- gaps_2 = gaps;
- tqse_2 = tqse;
- tqce_2 = tqce;
- tqcs_2 = tqcs;
- tqspe_2 = tqspe;
- tqcpe_2 = tqcpe;
- tqcps_2 = tqcps;
- tqss_2 = tqss;
- tqsps_2 = tqsps;
- zpv_2 = zpv;
- gapm_2 = gapm;
- gappm_2 = gappm;
- gap_2 = gap;
- gapp_2 = gapp;
- argi_2 = argi;
- args_2 = args;
- duk_2 = duk;
- cextlr_2 = cextlr;
- cext_2 = cext;
- cmullr_2 = cmullr;
- cmul_2 = cmul;
- gapv_2 = gapv;
- tqev_2 = tqev;
- tqsv_2 = tqsv;
- u_2 = u;
- us_2 = us;
- un_2 = un;
- uns_2 = uns;
- up_2 = up;
- ups_2 = ups;
- unmp_2 = unmp;
- unsmp_2 = unsmp;
- upmp_2 = upmp;
- upsmp_2 = upsmp;
- }
- else {
+ } else {
// this is not thread 0, so do create fresh copies of all local variables
- c1_2 = new C1(*c1);
- c1ao_2 = new C1_AddOns(*c1ao);
- c2_2 = new C2(*c2);
- c3_2 = new C3(*c3);
- c4_2 = new C4(*c4);
- c6_2 = new C6(*c6);
- c9_2 = new C9(*c9);
+ cid_2 = new ClusterIterationData(*cid);
output_2 = fopen((output_path + "/c_OCLU_" + to_string(myompthread)).c_str(), "w");
tppoanp_2 = new fstream;
tppoanp_2->open((output_path + "/c_TPPOAN_" + to_string(myompthread)).c_str(), ios::out | ios::binary);
- gaps_2 = new double[nsph]();
- tqse_2 = new double*[2];
- tqce_2 = new double*[2];
- tqcs_2 = new double*[2];
- tqspe_2 = new dcomplex*[2];
- tqcpe_2 = new dcomplex*[2];
- tqcps_2 = new dcomplex*[2];
- tqss_2 = new double*[2];
- tqsps_2 = new dcomplex*[2];
- for (int ti = 0; ti < 2; ti++) {
- tqse_2[ti] = new double[nsph]();
- tqce_2[ti] = new double[3]();
- tqcs_2[ti] = new double[3]();
- tqspe_2[ti] = new dcomplex[nsph]();
- tqcpe_2[ti] = new dcomplex[3]();
- tqcps_2[ti] = new dcomplex[3]();
- tqss_2[ti] = new double[nsph]();
- tqsps_2[ti] = new dcomplex[nsph]();
- for (int tj=0; tj<nsph; tj++) {
- tqse_2[ti][tj] = tqse[ti][tj];
- tqspe_2[ti][tj] = tqspe[ti][tj];
- tqss_2[ti][tj] = tqss[ti][tj];
- tqsps_2[ti][tj] = tqsps[ti][tj];
- }
- for (int tj=0; tj<3; tj++) {
- tqce_2[ti][tj] = tqce[ti][tj];
- tqcs_2[ti][tj] = tqcs[ti][tj];
- tqcpe_2[ti][tj] = tqcpe[ti][tj];
- tqcps_2[ti][tj] = tqcps[ti][tj];
- }
- }
- zpv_2 = new double***[c4->lm]; // [lm][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
- for (int zi = 0; zi < c4->lm; zi++) {
- zpv_2[zi] = new double**[3];
- for (int zj = 0; zj < 3; zj++) {
- zpv_2[zi][zj] = new double*[2];
- for (int zk = 0; zk < 2; zk++) {
- zpv_2[zi][zj][zk] = new double[2]();
- for (int zl = 0; zl < 2; zl++) zpv_2[zi][zj][zk][zl] = zpv[zi][zj][zk][zl];
- }
- }
- }
- gapm_2 = new double*[3];
- gappm_2 = new dcomplex*[3];
- gap_2 = new double*[3];
- gapp_2 = new dcomplex*[3];
- for (int gi = 0; gi < 3; gi++) {
- gap_2[gi] = new double[2]();
- gapm_2[gi] = new double[2]();
- gapp_2[gi] = new dcomplex[2]();
- gappm_2[gi] = new dcomplex[2]();
- for (int gj=0; gj<2; gj++) {
- gap_2[gi][gj] = gap[gi][gj];
- gapp_2[gi][gj] = gapp[gi][gj];
- gapm_2[gi][gj] = gapm[gi][gj];
- gappm_2[gi][gj] = gappm[gi][gj];
- }
- }
- argi_2 = new double[1]();
- argi_2[0] = argi[0];
- args_2 = new double[1]();
- args_2[0] = args[0];
- duk_2 = new double[3]();
- for (int di=0; di<3; di++) duk_2[di] = duk[di];
- cextlr_2 = new double*[4];
- cext_2 = new double*[4];
- cmullr_2 = new double*[4];
- cmul_2 = new double*[4];
- for (int ci = 0; ci < 4; ci++) {
- cextlr_2[ci] = new double[4]();
- cext_2[ci] = new double[4]();
- cmullr_2[ci] = new double[4]();
- cmul_2[ci] = new double[4]();
- for (int cj=0; cj<4; cj++) {
- cextlr_2[ci][cj] = cextlr[ci][cj];
- cext_2[ci][cj] = cext[ci][cj];
- cmullr_2[ci][cj] = cmullr[ci][cj];
- cmul_2[ci][cj] = cmul[ci][cj];
- }
- }
- gapv_2 = new double[3]();
- for (int gi=0; gi<3; gi++) gapv_2[gi] = gapv[gi];
- tqev_2 = new double[3]();
- tqsv_2 = new double[3]();
- for (int ti=0; ti<3; ti++) {
- tqev_2[ti] = tqev[ti];
- tqsv_2[ti] = tqsv[ti];
- }
- u_2 = new double[3]();
- us_2 = new double[3]();
- un_2 = new double[3]();
- uns_2 = new double[3]();
- up_2 = new double[3]();
- ups_2 = new double[3]();
- unmp_2 = new double[3]();
- unsmp_2 = new double[3]();
- upmp_2 = new double[3]();
- upsmp_2 = new double[3]();
- for (int ti=0; ti<3; ti++) {
- u_2[ti] = u[ti];
- us_2[ti] = us[ti];
- un_2[ti] = un[ti];
- uns_2[ti] = uns[ti];
- up_2[ti] = up[ti];
- ups_2[ti] = ups[ti];
- unmp_2[ti] = unmp[ti];
- unsmp_2[ti] = unsmp[ti];
- upmp_2[ti] = upmp[ti];
- upsmp_2[ti] = upsmp[ti];
- }
}
fstream &tppoan_2 = *tppoanp_2;
// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
@@ -511,87 +218,16 @@ void cluster(const string& config_file, const string& data_file, const string& o
// ok, now I can actually start the parallel calculations
#pragma omp for
for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
- int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, 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, 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, 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, arg_2, wn_2, vk_2, logger);
+ int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid_2, output_2, output_path, *tppoanp_2, logger);
}
#pragma omp barrier
// only threads different from 0 have to free local copies of variables and close local files
if (myompthread != 0) {
- delete c1_2;
- delete c1ao_2;
- delete c2_2;
- delete c3_2;
- delete c4_2;
- delete c6_2;
- delete c9_2;
+ delete cid_2;
fclose(output_2);
tppoanp_2->close();
delete tppoanp_2;
- delete[] gaps_2;
- for (int ti = 0; ti < 2; ti++) {
- delete[] tqse_2[ti];
- delete[] tqce_2[ti];
- delete[] tqcs_2[ti];
- delete[] tqspe_2[ti];
- delete[] tqcpe_2[ti];
- delete[] tqcps_2[ti];
- delete[] tqss_2[ti];
- delete[] tqsps_2[ti];
- }
- delete[] tqse_2;
- delete[] tqce_2;
- delete[] tqcs_2;
- delete[] tqspe_2;
- delete[] tqcpe_2;
- delete[] tqcps_2;
- delete[] tqss_2;
- delete[] tqsps_2;
- for (int zi = 0; zi < c4->lm; zi++) {
- for (int zj = 0; zj < 3; zj++) {
- for (int zk = 0; zk < 2; zk++) {
- delete[] zpv_2[zi][zj][zk];
- }
- delete[] zpv_2[zi][zj];
- }
- delete[] zpv_2[zi];
- }
- delete[] zpv_2;
- for (int gi = 0; gi < 3; gi++) {
- delete[] gap_2[gi];
- delete[] gapp_2[gi];
- delete[] gapm_2[gi];
- delete[] gappm_2[gi];
- }
- delete[] gap_2;
- delete[] gapp_2;
- delete[] gapm_2;
- delete[] gappm_2;
- delete[] argi_2;
- delete[] args_2;
- delete[] duk_2;
- for (int ci = 0; ci < 4; ci++) {
- delete[] cextlr_2[ci];
- delete[] cext_2[ci];
- delete[] cmullr_2[ci];
- delete[] cmul_2[ci];
- }
- delete[] cextlr_2;
- delete[] cext_2;
- delete[] cmullr_2;
- delete[] cmul_2;
- delete[] gapv_2;
- delete[] tqev_2;
- delete[] tqsv_2;
- delete[] u_2;
- delete[] us_2;
- delete[] un_2;
- delete[] uns_2;
- delete[] up_2;
- delete[] ups_2;
- delete[] unmp_2;
- delete[] unsmp_2;
- delete[] upmp_2;
- delete[] upsmp_2;
}
#pragma omp barrier
{
@@ -642,78 +278,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
fclose(output);
// Clean memory
- for (int zi = c4->lm - 1; zi > -1; zi--) {
- for (int zj = 2; zj > -1; zj--) {
- delete[] zpv[zi][zj][1];
- delete[] zpv[zi][zj][0];
- delete[] zpv[zi][zj];
- }
- delete[] zpv[zi];
- }
- delete[] zpv;
+ delete cid;
delete p_scattering_angles;
- delete c1;
- delete c1ao;
- delete c2;
- delete c3;
- delete c4;
- delete c6;
- delete c9;
- delete[] gaps;
- for (int ti = 1; ti > -1; ti--) {
- delete[] tqse[ti];
- delete[] tqss[ti];
- delete[] tqspe[ti];
- delete[] tqsps[ti];
- delete[] tqce[ti];
- delete[] tqcpe[ti];
- delete[] tqcs[ti];
- delete[] tqcps[ti];
- }
- delete[] tqse;
- delete[] tqss;
- delete[] tqspe;
- delete[] tqsps;
- delete[] tqce;
- delete[] tqcpe;
- delete[] tqcs;
- delete[] tqcps;
- delete[] tqev;
- delete[] tqsv;
- for (int gi = 2; gi > -1; gi--) {
- delete[] gapp[gi];
- delete[] gappm[gi];
- delete[] gap[gi];
- delete[] gapm[gi];
- }
- delete[] gapp;
- delete[] gappm;
- delete[] gap;
- delete[] gapm;
- delete[] gapv;
- delete[] u;
- delete[] us;
- delete[] un;
- delete[] uns;
- delete[] up;
- delete[] ups;
- delete[] unmp;
- delete[] unsmp;
- delete[] upmp;
- delete[] upsmp;
- delete[] argi;
- delete[] args;
- delete[] duk;
- for (int ci = 3; ci > -1; ci--) {
- delete[] cextlr[ci];
- delete[] cext[ci];
- delete[] cmullr[ci];
- delete[] cmul[ci];
- }
- delete[] cextlr;
- delete[] cext;
- delete[] cmullr;
- delete[] cmul;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
@@ -732,8 +298,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
delete logger;
}
-
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, const string& output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, 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, 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, dcomplex arg, double wn, double vk, Logger *logger)
+int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, FILE *output, const string& output_path, fstream& tppoan, Logger *logger)
{
int nxi = sconf->number_of_scales;
logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
@@ -753,7 +318,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int npntts = gconf->npntts;
int isam = gconf->iavm;
int jwtm = gconf->jwtm;
- np_int ndit = 2 * nsph * c4->nlim;
+ np_int ndit = 2 * nsph * cid->c4->nlim;
int isq, ibf;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
@@ -763,43 +328,43 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int idfc = (int)sconf->idfc;
double vkarg = 0.0;
if (idfc >= 0) {
- vk = xi * wn;
- vkarg = vk;
- fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
+ cid->vk = xi * cid->wn;
+ vkarg = cid->vk;
+ fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", cid->vk, xi);
} else {
- vkarg = xi * vk;
- sqsfi = 1.0 / (xi * xi);
+ vkarg = xi * cid->vk;
+ cid->sqsfi = 1.0 / (xi * xi);
fprintf(output, " XI=%15.7lE\n", xi);
}
- hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
+ hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4);
if (jer != 0) {
fprintf(output, " STOP IN HJV\n");
return jer;
// break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
}
for (int i132 = 1; i132 <= nsph; i132++) {
- int iogi = c1->iog[i132 - 1];
+ int iogi = cid->c1->iog[i132 - 1];
if (iogi != i132) {
for (int l123 = 1; l123 <= li; l123++) {
- c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
- c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
+ cid->c1->rmi[l123 - 1][i132 - 1] = cid->c1->rmi[l123 - 1][iogi - 1];
+ cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1];
} // l123 loop
} else {
- int nsh = c1->nshl[i132 - 1];
+ int nsh = cid->c1->nshl[i132 - 1];
int ici = (nsh + 1) / 2;
if (idfc == 0) {
for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
+ cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
} else {
if (jxi488 == 1) {
for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
+ cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
}
}
- if (nsh % 2 == 0) c2->dc0[ici] = exdc;
+ if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;
dme(
- c4->li, i132, npnt, npntts, vkarg, exdc, exri,
- c1, c2, jer, lcalc, arg
+ cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri,
+ cid->c1, cid->c2, jer, lcalc, cid->arg
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
@@ -812,29 +377,20 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
//break;
}
} // i132 loop
- dcomplex *am_vector = new dcomplex[ndit * ndit]();
- dcomplex **am = new dcomplex*[ndit];
- for (int ai = 0; ai < ndit; ai++) {
- am[ai] = (am_vector + ai * ndit);
- }
- cms(am, c1, c1ao, c4, c6);
- invert_matrix(am, ndit, jer, mxndm);
+ cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
+ invert_matrix(cid->am, ndit, jer, mxndm);
if (jer != 0) {
- delete[] am_vector;
- delete[] am;
return jer;
// break; // jxi488 loop: goes to memory clean
}
- ztm(am, c1, c1ao, c4, c6, c9);
- delete[] am_vector;
- delete[] am;
+ ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
if (idfc >= 0) {
if (jxi488 == jwtm) {
- int nlemt = 2 * c4->nlem;
+ int nlemt = 2 * cid->c4->nlem;
string ttms_name = output_path + "/c_TTMS.hd5";
- TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m, "HDF5");
+ TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5");
ttms_name = output_path + "/c_TTMS";
- TransitionMatrix::write_binary(ttms_name, nlemt, lm, vk, exri, c1ao->am0m);
+ TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m);
}
}
// label 156: continue from here
@@ -844,56 +400,56 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
fprintf(output, " CIRC\n");
}
// label 160
- double cs0 = 0.25 * vk * vk * vk / acos(0.0);
+ double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
dcomplex s0 = 0.0 + 0.0 * I;
- scr0(vk, exri, c1, c1ao, c3, c4);
- double sqk = vk * vk * exdc;
- aps(zpv, c4->li, nsph, c1, sqk, gaps);
- rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
- if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
+ scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4);
+ double sqk = cid->vk * cid->vk * exdc;
+ aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps);
+ rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);
+ if (cid->c4->li != cid->c4->le) fprintf(output, " SPHERES; LMX=LI\n");
for (int i170 = 1; i170 <= nsph; i170++) {
- if (c1->iog[i170 - 1] >= i170) {
+ if (cid->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;
+ double albeds = cid->c1->sscs[i] / cid->c1->sexs[i];
+ cid->c1->sqscs[i] *= cid->sqsfi;
+ cid->c1->sqabs[i] *= cid->sqsfi;
+ cid->c1->sqexs[i] *= cid->sqsfi;
fprintf(output, " SPHERE %2d\n", i170);
- if (c1->nshl[i] != 1) {
- fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
+ if (cid->c1->nshl[i] != 1) {
+ fprintf(output, " SIZE=%15.7lE\n", cid->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]));
+ fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->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, " %14.7lE%15.7lE%15.7lE%15.7lE\n", cid->c1->sscs[i], cid->c1->sabs[i], cid->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;
+ fprintf(output, " %14.7lE%15.7lE%15.7lE\n", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]);
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i]));
+ csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcsv[i];
+ s0 = cid->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];
+ double rapr = cid->c1->sexs[i] - cid->gaps[i];
+ double cosav = cid->gaps[i] / cid->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]);
+ fprintf(output, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", cid->tqse[0][i], cid->tqss[0][i]);
+ fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", cid->tqse[1][i], cid->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;
+ fprintf(output, " FSAT=%15.7lE%15.7lE\n", real(cid->c3->tfsas), imag(cid->c3->tfsas));
+ csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs;
+ s0 = cid->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);
+ tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double));
+ pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4);
+ apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm);
double th = sa->th;
for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
double ph = sa->ph;
@@ -901,22 +457,22 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
int jw = 0;
if (sa->nk != 1 || jxi488 <= 1) {
- upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
+ upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
if (isam >= 0) {
wmamp(
- 0, cost, sint, cosp, sinp, inpol, c4->le, 0,
- nsph, argi, u, upmp, unmp, c1
+ 0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0,
+ nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->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);
+ apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
+ raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->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);
+ apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
+ raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
jw = 1;
}
}
@@ -945,39 +501,39 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
// label 188
bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
if (!goto190) {
- upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
+ upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
if (isam >= 0)
wmamp(
- 2, costs, sints, cosps, sinps, inpol, c4->le,
- 0, nsph, args, us, upsmp, unsmp, c1
+ 2, costs, sints, cosps, sinps, inpol, cid->c4->le,
+ 0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1
);
}
// label 190
if (sa->nkks != 1 || jxi488 <= 1) {
upvsp(
- u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
- duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
+ cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,
+ cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->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
+ cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le,
+ 0, nsph, cid->argi, cid->args, cid->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];
+ cid->up[i193] = cid->upmp[i193];
+ cid->un[i193] = cid->unmp[i193];
+ cid->ups[i193] = cid->upsmp[i193];
+ cid->uns[i193] = cid->unsmp[i193];
}
}
}
// label 194
- if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
+ if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->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);
+ apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
+ raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
jw = 1;
}
// label 196
@@ -985,38 +541,38 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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));
+ tppoan.write(reinterpret_cast<char *>(&(cid->scan)), sizeof(double));
if (jaw != 0) {
jaw = 0;
- mextc(vk, exri, c1ao->fsacm, cextlr, cext);
+ mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->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];
+ double value = cid->cext[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
- double value = c1ao->scscm[i];
+ double value = cid->c1ao->scscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(c1ao->scscpm[i]);
+ value = real(cid->c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->scscpm[i]);
+ value = imag(cid->c1ao->scscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecscm[i];
+ value = cid->c1ao->ecscm[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(c1ao->ecscpm[i]);
+ value = real(cid->c1ao->ecscpm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->ecscpm[i]);
+ value = imag(cid->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];
+ double value = cid->gapm[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(gappm[i][j]);
+ value = real(cid->gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(gappm[i][j]);
+ value = imag(cid->gappm[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -1025,24 +581,24 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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 extsm = cid->c1ao->ecscm[ilr210 - 1];
+ double qextm = extsm * cid->sqsfi / cid->c3->gcs;
+ double extrm = extsm / cid->c3->ecs;
+ double scasm = cid->c1ao->scscm[ilr210 - 1];
double albdm = scasm / extsm;
- double qscam = scasm *sqsfi / c3->gcs;
- double scarm = scasm / c3->scs;
+ double qscam = scasm * cid->sqsfi / cid->c3->gcs;
+ double scarm = scasm / cid->c3->scs;
double abssm = extsm - scasm;
- double qabsm = abssm * sqsfi / c3->gcs;
- double absrm = abssm / c3->acs;
- double acsecs = c3->acs / c3->ecs;
+ double qabsm = abssm * cid->sqsfi / cid->c3->gcs;
+ double absrm = abssm / cid->c3->acs;
+ double acsecs = cid->c3->acs / cid->c3->ecs;
if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
- dcomplex s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
+ dcomplex s0m = cid->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);
+ double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas);
+ double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 206
@@ -1057,37 +613,37 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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])
+ ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
+ imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
+ real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->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 rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
+ double cosav = cid->gapm[2][ilr210 - 1] / cid->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]);
+ double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]);
+ double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->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);
+ fprintf(output, " SCAND=%10.3lE\n", cid->scan);
+ fprintf(output, " CFMP=%15.7lE, SFMP=%15.7lE\n", cid->cfmp, cid->sfmp);
+ fprintf(output, " CFSP=%15.7lE, SFSP=%15.7lE\n", cid->cfsp, cid->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]);
+ fprintf(output, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", cid->un[0], cid->un[1], cid->un[2]);
+ fprintf(output, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", cid->uns[0], cid->uns[1], cid->uns[2]);
} else { // label 214
- fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
+ fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", cid->un[0], cid->un[1], cid->un[2]);
}
// label 220
if (inpol == 0) {
@@ -1096,127 +652,127 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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");
+ scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4);
+ if (cid->c4->li != cid->c4->le) fprintf(output, " SPHERES; LMX=MIN0(LI,LE)\n");
for (int i226 = 1; i226 <= nsph; i226++) {
- if (c1->iog[i226 - 1] >= i226) {
+ if (cid->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])
+ real(cid->c1->sas[i226 - 1][0][0]), imag(cid->c1->sas[i226 - 1][0][0]),
+ real(cid->c1->sas[i226 - 1][1][0]), imag(cid->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])
+ real(cid->c1->sas[i226 - 1][0][1]), imag(cid->c1->sas[i226 - 1][0][1]),
+ real(cid->c1->sas[i226 - 1][1][1]), imag(cid->c1->sas[i226 - 1][1][1])
);
for (int j225 = 0; j225 < 16; j225++) {
- c1->vint[j225] = c1->vints[i226 - 1][j225];
+ cid->c1->vint[j225] = cid->c1->vints[i226 - 1][j225];
} // j225 loop
- mmulc(c1->vint, cmullr, cmul);
+ mmulc(cid->c1->vint, cid->cmullr, cid->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]
+ cid->cmul[i1][0], cid->cmul[i1][1], cid->cmul[i1][2], cid->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]
+ cid->cmullr[i1][0], cid->cmullr[i1][1], cid->cmullr[i1][2], cid->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])
+ real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]),
+ real(cid->c3->tsas[1][0]), imag(cid->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])
+ real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]),
+ real(cid->c3->tsas[1][1]), imag(cid->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);
+ pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4);
+ mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext);
+ mmulc(cid->c1->vint, cid->cmullr, cid->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];
+ double value = cid->cext[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 2; i++) {
- double value = c1ao->scsc[i];
+ double value = cid->c1ao->scsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(c1ao->scscp[i]);
+ value = real(cid->c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->scscp[i]);
+ value = imag(cid->c1ao->scscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = c1ao->ecsc[i];
+ value = cid->c1ao->ecsc[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(c1ao->ecscp[i]);
+ value = real(cid->c1ao->ecscp[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->ecscp[i]);
+ value = imag(cid->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];
+ double value = cid->gap[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(gapp[i][j]);
+ value = real(cid->gapp[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(gapp[i][j]);
+ value = imag(cid->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];
+ double value = cid->tqce[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(tqcpe[i][j]);
+ value = real(cid->tqcpe[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(tqcpe[i][j]);
+ value = imag(cid->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];
+ double value = cid->tqcs[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = real(tqcps[i][j]);
+ value = real(cid->tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(tqcps[i][j]);
+ value = imag(cid->tqcps[i][j]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
for (int i = 0; i < 3; i++) {
- double value = u[i];
+ double value = cid->u[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = up[i];
+ value = cid->up[i];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = un[i];
+ value = cid->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]);
+ double value = real(cid->c1->vint[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1->vint[i]);
+ value = imag(cid->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];
+ double value = cid->cmul[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -1224,24 +780,24 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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 extsec = cid->c1ao->ecsc[ilr290 - 1];
+ double qext = extsec * cid->sqsfi / cid->c3->gcs;
+ double extrat = extsec / cid->c3->ecs;
+ double scasec = cid->c1ao->scsc[ilr290 - 1];
double albedc = scasec / extsec;
- double qsca = scasec * sqsfi / c3->gcs;
- double scarat = scasec / c3->scs;
+ double qsca = scasec * cid->sqsfi / cid->c3->gcs;
+ double scarat = scasec / cid->c3->scs;
double abssec = extsec - scasec;
- double qabs = abssec * sqsfi / c3->gcs;
+ double qabs = abssec * cid->sqsfi / cid->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 ratio = cid->c3->acs / cid->c3->ecs;
+ if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs;
+ s0 = cid->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);
+ double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);
+ double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);
if (inpol == 0) {
fprintf(output, " LIN %2d\n", ipol);
} else { // label 273
@@ -1265,13 +821,13 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
);
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])
+ ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->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])
+ ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]),
+ jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->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",
@@ -1283,67 +839,67 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
);
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];
+ cid->gapv[0] = cid->gap[0][ilr290 - 1];
+ cid->gapv[1] = cid->gap[1][ilr290 - 1];
+ cid->gapv[2] = cid->gap[2][ilr290 - 1];
+ double extins = cid->c1ao->ecsc[ilr290 - 1];
+ double scatts = cid->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);
+ rftr(cid->u, cid->up, cid->un, cid->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];
+ cid->tqev[0] = cid->tqce[ilr290 - 1][0];
+ cid->tqev[1] = cid->tqce[ilr290 - 1][1];
+ cid->tqev[2] = cid->tqce[ilr290 - 1][2];
+ cid->tqsv[0] = cid->tqcs[ilr290 - 1][0];
+ cid->tqsv[1] = cid->tqcs[ilr290 - 1][1];
+ cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
double tep, ten, tek, tsp, tsn, tsk;
- tqr(u, up, un, tqev, tqsv, tep, ten, tek, tsp, tsn, tsk);
+ tqr(cid->u, cid->up, cid->un, cid->tqev, cid->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]
+ cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->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]
+ cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->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]);
+ double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]);
+ double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->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]
+ cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->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]
+ cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
);
}
if (iavm != 0) {
- mmulc(c1ao->vintm, cmullr, cmul);
+ mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul);
// Some implicit loops writing to binary.
for (int i = 0; i < 16; i++) {
double value;
- value = real(c1ao->vintm[i]);
+ value = real(cid->c1ao->vintm[i]);
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
- value = imag(c1ao->vintm[i]);
+ value = imag(cid->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];
+ double value = cid->cmul[i][j];
tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
}
}
@@ -1358,14 +914,14 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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]
+ cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->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]
+ cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
);
}
}
@@ -1378,20 +934,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
} // jph484 loop
th += sa->thstp;
} // jth486 loop
-
- // delete[] u;
- // delete[] us;
- // delete[] un;
- // delete[] uns;
- // delete[] up;
- // delete[] ups;
- // delete[] unmp;
- // delete[] unsmp;
- // delete[] upmp;
- // delete[] upsmp;
-
+
logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
return jer;
-
}
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index e74bfa662c96e5ba22ae8914ebf379ab9ecc76c8..1a4daaeb107360550b5ca69d62f91288ab224152 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -794,13 +794,11 @@ ClusterIterationData::~ClusterIterationData() {
delete c6;
delete c9;
delete[] gaps;
- for (int ti = nsph -1; ti > -1; ti--) {
+ for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
- }
- for (int ti = 1; ti > -1; ti--) {
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];