diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 8494bdd310826e4ac83c5e04b32f11f3e7b13636..4c6cbeaa955d481a75533419da69620d56bc5ded 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -86,36 +86,40 @@ void cluster(string config_file, string data_file, string output_path) {
exit(1);
}
logger->log(" done.\n", LOG_INFO);
- if (sconf->number_of_spheres == gconf->number_of_spheres) {
+ int s_nsph = sconf->get_param("nsph");
+ int nsph = gconf->get_param("nsph");
+ if (s_nsph == nsph) {
// Shortcuts to variables stored in configuration objects
- int nsph = gconf->number_of_spheres;
- np_int mxndm = (np_int)gconf->mxndm;
- int inpol = gconf->in_pol;
- int npnt = gconf->npnt;
- int npntts = gconf->npntts;
- int iavm = gconf->iavm;
- int isam = gconf->meridional_type;
- int nxi = sconf->number_of_scales;
- 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 wp = sconf->wp;
+ np_int mxndm = (np_int)gconf->get_param("mxndm");
+ int inpol = (int)gconf->get_param("in_pol");
+ int npnt = (int)gconf->get_param("npnt");
+ int npntts = (int)gconf->get_param("npntts");
+ int iavm = (int)gconf->get_param("iavm");
+ int isam = (int)gconf->get_param("meridional_type");
+ int nxi = (int)sconf->get_param("number_of_scales");
+ int idfc = (int)sconf->get_param("idfc");
+ double th = gconf->get_param("in_theta_start");
+ double thstp = gconf->get_param("in_theta_step");
+ double thlst = gconf->get_param("in_theta_end");
+ double ths = gconf->get_param("sc_theta_start");
+ double thsstp = gconf->get_param("sc_theta_step");
+ double thslst = gconf->get_param("sc_theta_end");
+ double ph = gconf->get_param("in_phi_start");
+ double phstp = gconf->get_param("in_phi_step");
+ double phlst = gconf->get_param("in_phi_end");
+ double phs = gconf->get_param("sc_phi_start");
+ double phsstp = gconf->get_param("sc_phi_step");
+ double phslst = gconf->get_param("sc_phi_end");
+ double wp = sconf->get_param("wp");
// Global variables for CLU
- int lm = gconf->l_max;
- if (gconf->li > lm) lm = gconf->li;
- if (gconf->le > lm) lm = gconf->le;
+ int lm = (int)gconf->get_param("l_max");
+ int li = (int)gconf->get_param("li");
+ int le = (int)gconf->get_param("le");
+ if (li > lm) lm = li;
+ if (le > lm) lm = le;
C1 *c1 = new C1(gconf, sconf);
C3 *c3 = new C3();
- C4 *c4 = new C4(gconf->li, gconf->le, nsph);
+ C4 *c4 = new C4(li, le, nsph);
C1_AddOns *c1ao = new C1_AddOns(c4);
// End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo);
@@ -124,13 +128,10 @@ void cluster(string config_file, string data_file, string output_path) {
int lcalc = 0;
dcomplex arg = 0.0 + 0.0 * I;
dcomplex ccsam = 0.0 + 0.0 * I;
- int configurations = 0;
- for (int ci = 1; ci <= nsph; ci++) {
- if (sconf->iog_vec[ci -1] >= ci) configurations++;
- }
+ int configurations = (int)sconf->get_param("configurations");
C2 *c2 = new C2(nsph, configurations, npnt, npntts);
np_int ndit = 2 * nsph * c4->nlim;
- printf("Size of matrices to invert: %ld x %ld\n", (int64_t) ndit, (int64_t) ndit);
+ logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
dcomplex *am_vector = new dcomplex[ndit * ndit]();
dcomplex **am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) {
@@ -200,6 +201,7 @@ void cluster(string config_file, string data_file, string output_path) {
cmul[ci] = new double[4]();
}
int isq, ibf;
+ int jwtm = (int)gconf->get_param("jwtm");
double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
@@ -208,7 +210,7 @@ void cluster(string config_file, string data_file, string output_path) {
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
- gconf->sph_x[ri], gconf->sph_y[ri], gconf->sph_z[ri]
+ gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
);
fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf(
@@ -221,7 +223,7 @@ void cluster(string config_file, string data_file, string output_path) {
ph, phstp, phlst, phs, phsstp, phslst
);
fprintf(output, " READ(IR,*)JWTM\n");
- fprintf(output, " %5d\n", gconf->jwtm);
+ fprintf(output, " %5d\n", jwtm);
fprintf(output, " READ(ITIN)NSPHT\n");
fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
@@ -256,7 +258,7 @@ void cluster(string config_file, string data_file, string output_path) {
int nks = nths * nphs;
int nkks = nk * nks;
double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
- str(sconf->rcf, c1, c1ao, c3, c4, c6);
+ 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];
@@ -268,7 +270,8 @@ void cluster(string config_file, string data_file, string output_path) {
}
}
thdps(c4->lm, zpv);
- double exri = sqrt(sconf->exdc);
+ double exdc = sconf->get_param("exdc");
+ double exri = sqrt(exdc);
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream *tppoanp = new fstream;
@@ -290,9 +293,10 @@ void cluster(string config_file, string data_file, string output_path) {
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
double wn = wp / 3.0e8;
+ double xip = sconf->get_param("xip");
double sqsfi = 1.0;
- if (sconf->idfc < 0) {
- vk = sconf->xip * wn;
+ if (idfc < 0) {
+ vk = xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
@@ -501,7 +505,7 @@ void cluster(string config_file, string data_file, string output_path) {
tqcps_2[ti][tj] = tqcps[ti][tj];
}
}
- zpv_2 = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
+ 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++) {
@@ -589,7 +593,7 @@ void cluster(string config_file, string data_file, string output_path) {
fstream &tppoan_2 = *tppoanp_2;
// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
#pragma omp barrier
- if (myompthread==0) printf("Syncing OpenMP threads and starting the loop on wavelengths\n");
+ if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengthsn\n");
// ok, now I can actually start the parallel calculations
#pragma omp for
for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
@@ -815,74 +819,18 @@ void cluster(string config_file, string data_file, string output_path) {
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, Logger *logger)
{
-
- // 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 jer = 0;
int jaw = 1;
+ int jwtm = (int)gconf->get_param("jwtm");
+ int li = (int)gconf->get_param("li");
+ int le = (int)gconf->get_param("le");
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
- double xi = sconf->scale_vec[jxi488 - 1];
+ double xi = sconf->get_scale(jxi488 - 1);
+ double exdc = sconf->get_param("exdc");
+ int idfc = (int)sconf->get_param("idfc");
double vkarg = 0.0;
- if (sconf->idfc >= 0) {
+ if (idfc >= 0) {
vk = xi * wn;
vkarg = vk;
fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
@@ -894,43 +842,31 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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++) {
+ for (int l123 = 1; l123 <= li; l123++) {
c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
} // l123 loop
} else {
- int nsh = sconf->nshl_vec[i132 - 1];
+ int nsh = c1->nshl[i132 - 1];
int ici = (nsh + 1) / 2;
- if (sconf->idfc == 0) {
+ if (idfc == 0) {
for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+ 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->dc0_matrix[ic][i132 - 1][0];
+ c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
}
}
- if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+ if (nsh % 2 == 0) c2->dc0[ici] = exdc;
dme(
- c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,
+ c4->li, i132, npnt, npntts, vkarg, exdc, exri,
c1, c2, jer, lcalc, arg
);
if (jer != 0) {
@@ -987,8 +923,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
// break; // jxi488 loop: goes to memory clean
}
ztm(am, c1, c1ao, c4, c6, c9);
- if (sconf->idfc >= 0) {
- if (jxi488 == gconf->jwtm) {
+ if (idfc >= 0) {
+ if (jxi488 == jwtm) {
int nlemt = 2 * c4->nlem;
string ttms_name = output_path + "/c_TTMS.hd5";
TransitionMatrix ttms(nlemt, lm, vk, exri, c1ao->am0m);
@@ -1008,7 +944,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
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;
+ 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");
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 81254833b4cffffd42bb1f4f586f639853091b22..4b80aca1056df77fdd5d57bbadfdf6560dec4cce 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -30,19 +30,6 @@
#ifndef INCLUDE_CONFIGURATION_H_
#define INCLUDE_CONFIGURATION_H_
-class ScattererConfiguration;
-class GeometryConfiguration;
-class Logger;
-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.
*
@@ -52,9 +39,6 @@ class C9;
* fields and their polarization properties.
*/
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, Logger *);
protected:
//! \brief Number of spherical components.
@@ -228,9 +212,6 @@ public:
* data to describe the scatterer properties.
*/
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, Logger *);
protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 4621148a45de3fcd4b0107d449dbff454fd66ed3..2360c0607913b3f2939c7ce3a4219943f02bdc5d 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -344,14 +344,14 @@ void scr2(
* to radial coordinates, then it calls `sphar()` to calculate the vector of spherical
* harmonics of the incident field.
*
- * \param rcf: `double **` Matrix of double.
+ * \param sconf: `ScattererConfiguration *` Pointer to scatterer configuration object.
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
* \param c3: `C3 *` Pointer to a C3 instance.
* \param c4: `C4 *` Pointer to a C4 structure.
* \param c6: `C6 *` Pointer to a C6 instance.
*/
-void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
+void str(ScattererConfiguration *sconf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
/*! \brief Compute radiation torques on particles on a k-vector oriented system.
*
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index e6a24d5f3b3dca325026e646c67a34c1dcaa5be4..937ac24d4c512133c1d6c224fd9f675a9698a9a0 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1818,7 +1818,7 @@ void scr2(
} // ipo1 loop
}
-void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
+void str(ScattererConfiguration *sconf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
dcomplex *ylm;
const double pi = acos(-1.0);
c3->gcs = 0.0;
@@ -1830,7 +1830,7 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
c1->gcsv[i18 - 1] = gcss;
int nsh = c1->nshl[i18 - 1];
for (int j16 = 1; j16 <= nsh; j16++) {
- c1->rc[i18 - 1][j16 - 1] = rcf[i18 - 1][j16 - 1] * c1->ros[i18 - 1];
+ c1->rc[i18 - 1][j16 - 1] = sconf->get_rcf(i18 - 1, j16 - 1) * c1->ros[i18 - 1];
} // j16 loop
}
c3->gcs += gcss;