src/include/Commons.h

@@ -386,6 +386,85 @@ public:
   ~C9();
 };
 
+/*! \brief A data structure representing the information used for a single scale
+ * of the CLUSTER case.
+ */
+class ClusterIterationData {
+public:
+  //! \brief Pointer to a C1 structure.
+  C1 *c1;
+  //! \brief Pointer to a C1_AddOns structure.
+  C1_AddOns *c1ao;
+  //! \brief Pointer to a C2 structure.
+  C2 *c2;
+  //! \brief Pointer to a C3 structure.
+  C3 *c3;
+  //! brief Pointer to a C4 structure.
+  C4 *c4;
+  //! \brief Pointer to a C6 structure.
+  C6 *c6;
+  //! \brief Pointer to a C9 structure.
+  C9 *c9;
+  //! \brief Pointer to a formatted output file.
+  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;
+  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 qsfi;
+  double sqsfi;
+  dcomplex *am_vector;
+  dcomplex **am;
+  dcomplex arg;
+  //! \brief Wave vector.
+  double vk;
+  //! \brief Wave number.
+  double wn;
+  double xip;
+  //! \brief Pointer to a logger instance.
+  Logger *logger;
+
+  ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
+
+  ClusterIterationData(const ClusterIterationData& rhs);
+
+  ~ClusterIterationData();
+};
+
 /*! \brief A data structure representing the angles to be evaluated in the problem.
  *
  */

src/libnptm/Commons.cpp

@@ -8,6 +8,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif
@@ -529,6 +533,323 @@ C9::~C9() {
   delete[] sam;
 }
 
+ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
+  c1 = new C1(gconf, sconf);
+  c2 = new C2(gconf, sconf);
+  c3 = new C3();
+  c4 = new C4(gconf);
+  c1ao = new C1_AddOns(c4);
+  c6 = new C6(c4->lmtpo);
+  const int ndi = c4->nsph * c4->nlim;
+  const np_int ndit = 2 * ndi;
+  c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
+  gaps = new double[c4->nsph]();
+  tqev = new double[3]();
+  tqsv = new double[3]();
+  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[c4->nsph]();
+    tqspe[ti] = new dcomplex[c4->nsph]();
+    tqss[ti] = new double[c4->nsph]();
+    tqsps[ti] = new dcomplex[c4->nsph]();
+    tqce[ti] = new double[3]();
+    tqcpe[ti] = new dcomplex[3]();
+    tqcs[ti] = new double[3]();
+    tqcps[ti] = new dcomplex[3]();
+  }
+  gapv = new double[3]();
+  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]();
+  }
+  u = new double[3]();
+  us = new double[3]();
+  un = new double[3]();
+  uns = new double[3]();
+  up = new double[3]();
+  ups = new double[3]();
+  unmp = new double[3]();
+  unsmp = new double[3]();
+  upmp = new double[3]();
+  upsmp = new double[3]();
+  argi = new double[1]();
+  args = new double[1]();
+  duk = new double[3]();
+  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]();
+  }
+  zpv = new double***[c4->lm];
+  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]();
+      }
+    }
+  }
+  am_vector = new dcomplex[ndit * ndit]();
+  am = new dcomplex*[ndit];
+  for (int ai = 0; ai < ndit; ai++) {
+    am[ai] = (am_vector + ai * ndit);
+  }
+  
+  arg = 0.0 + 0.0 * I;
+  // These are suspect initializations
+  scan = 0.0;
+  cfmp = 0.0;
+  sfmp = 0.0;
+  cfsp = 0.0;
+  sfsp = 0.0;
+  qsfi = 0.0;
+  // End of suspect initializations
+  wn = sconf->wp / 3.0e8;
+  xip = sconf->xip;
+  sqsfi = 1.0;
+  vk = 0.0;
+}
+
+ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
+  c1 = new C1(*(rhs.c1));
+  c2 = new C2(*(rhs.c2));
+  c3 = new C3(*(rhs.c3));
+  c4 = new C4(*(rhs.c4));
+  c1ao = new C1_AddOns(*(rhs.c1ao));
+  c6 = new C6(*(rhs.c6));
+  const int ndi = c4->nsph * c4->nlim;
+  const np_int ndit = 2 * ndi;
+  c9 = new C9(*(rhs.c9));
+  gaps = new double[c4->nsph]();
+  for (int gi = 0; gi < c4->nsph; gi++) gaps[gi] = rhs.gaps[gi];
+  tqev = new double[3]();
+  tqsv = new double[3]();
+  for (int ti = 0; ti < 3; ti++) {
+    tqev[ti] = rhs.tqev[ti];
+    tqsv[ti] = rhs.tqsv[ti];
+  }
+  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[c4->nsph]();
+    tqspe[ti] = new dcomplex[c4->nsph]();
+    tqss[ti] = new double[c4->nsph]();
+    tqsps[ti] = new dcomplex[c4->nsph]();
+    for (int tj = 0; tj < c4->nsph; tj++) {
+      tqse[ti][tj] = rhs.tqse[ti][tj];
+      tqspe[ti][tj] = rhs.tqspe[ti][tj];
+      tqss[ti][tj] = rhs.tqss[ti][tj];
+      tqsps[ti][tj] = rhs.tqsps[ti][tj];
+    }
+    tqce[ti] = new double[3]();
+    tqcpe[ti] = new dcomplex[3]();
+    tqcs[ti] = new double[3]();
+    tqcps[ti] = new dcomplex[3]();
+    for (int tj = 0; tj < 3; tj++) {
+      tqce[ti][tj] = rhs.tqce[ti][tj];
+      tqcpe[ti][tj] = rhs.tqcpe[ti][tj];
+      tqcs[ti][tj] = rhs.tqcs[ti][tj];
+      tqcps[ti][tj] = rhs.tqcps[ti][tj];
+    }
+  }
+  gapv = new double[3]();
+  gapp = new dcomplex*[3];
+  gappm = new dcomplex*[3];
+  gap = new double*[3];
+  gapm = new double*[3];
+  for (int gi = 0; gi < 3; gi++) {
+    gapv[gi] = rhs.gapv[gi];
+    gapp[gi] = new dcomplex[2]();
+    gappm[gi] = new dcomplex[2]();
+    gap[gi] = new double[2]();
+    gapm[gi] = new double[2]();
+    for (int gj = 0; gj < 2; gj++) {
+      gapp[gi][gj] = rhs.gapp[gi][gj];
+      gappm[gi][gj] = rhs.gappm[gi][gj];
+      gap[gi][gj] = rhs.gap[gi][gj];
+      gapm[gi][gj] = rhs.gapm[gi][gj];
+    }
+  }
+  u = new double[3]();
+  us = new double[3]();
+  un = new double[3]();
+  uns = new double[3]();
+  up = new double[3]();
+  ups = new double[3]();
+  unmp = new double[3]();
+  unsmp = new double[3]();
+  upmp = new double[3]();
+  upsmp = new double[3]();
+  duk = new double[3]();
+  for (int ui = 0; ui < 3; ui++) {
+    u[ui] = rhs.u[ui];
+    us[ui] = rhs.us[ui];
+    un[ui] = rhs.un[ui];
+    uns[ui] = rhs.uns[ui];
+    up[ui] = rhs.up[ui];
+    ups[ui] = rhs.ups[ui];
+    unmp[ui] = rhs.unmp[ui];
+    unsmp[ui] = rhs.unsmp[ui];
+    upmp[ui] = rhs.upmp[ui];
+    upsmp[ui] = rhs.upsmp[ui];
+    duk[ui] = rhs.duk[ui];
+  }
+  argi = new double[1]();
+  args = new double[1]();
+  argi[0] = rhs.argi[0];
+  args[0] = rhs.args[0];
+  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]();
+    for (int cj = 0; cj < 4; cj++) {
+      cextlr[ci][cj] = rhs.cextlr[ci][cj];
+      cext[ci][cj] = rhs.cext[ci][cj];
+      cmullr[ci][cj] = rhs.cmullr[ci][cj];
+      cmul[ci][cj] = rhs.cmul[ci][cj];
+    }
+  }
+  zpv = new double***[c4->lm];
+  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]();
+	zpv[zi][zj][zk][0] = rhs.zpv[zi][zj][zk][0];
+	zpv[zi][zj][zk][1] = rhs.zpv[zi][zj][zk][1];
+      }
+    }
+  }
+  am_vector = new dcomplex[ndit * ndit]();
+  for (np_int ai = 0; ai < ndit * ndit; ai++) am_vector[ai] = rhs.am_vector[ai];
+  am = new dcomplex*[ndit];
+  for (np_int ai = 0; ai < ndit; ai++) {
+    am[ai] = (am_vector + ai * ndit);
+  }
+  
+  arg = rhs.arg;
+  // These are suspect initializations
+  scan = rhs.scan;
+  cfmp = rhs.cfmp;
+  sfmp = rhs.sfmp;
+  cfsp = rhs.cfsp;
+  sfsp = rhs.sfsp;
+  qsfi = rhs.qsfi;
+  // End of suspect initializations
+  wn = rhs.wn;
+  xip = rhs.xip;
+  sqsfi = rhs.sqsfi;
+  vk = rhs.vk;
+}
+
+ClusterIterationData::~ClusterIterationData() {
+  const int nsph = c4->nsph;
+  delete[] am_vector;
+  delete[] am;
+  for (int zi = c4->lm - 1; zi > -1; zi--) {
+    for (int zj = 2; zj > -1; zj--) {
+      delete[] zpv[zi][zj][1];
+      delete[] zpv[zi][zj][0];
+      delete[] zpv[zi][zj];
+    }
+    delete[] zpv[zi];
+  }
+  delete[] zpv;
+  delete c1;
+  delete c1ao;
+  delete c2;
+  delete c3;
+  delete c4;
+  delete c6;
+  delete c9;
+  delete[] gaps;
+  for (int ti = 1; ti > -1; ti--) {
+    delete[] tqse[ti];
+    delete[] tqss[ti];
+    delete[] tqspe[ti];
+    delete[] tqsps[ti];
+    delete[] tqce[ti];
+    delete[] tqcpe[ti];
+    delete[] tqcs[ti];
+    delete[] tqcps[ti];
+  }
+  delete[] tqse;
+  delete[] tqss;
+  delete[] tqspe;
+  delete[] tqsps;
+  delete[] tqce;
+  delete[] tqcpe;
+  delete[] tqcs;
+  delete[] tqcps;
+  delete[] tqev;
+  delete[] tqsv;
+  for (int gi = 2; gi > -1; gi--) {
+    delete[] gapp[gi];
+    delete[] gappm[gi];
+    delete[] gap[gi];
+    delete[] gapm[gi];
+  }
+  delete[] gapp;
+  delete[] gappm;
+  delete[] gap;
+  delete[] gapm;
+  delete[] gapv;
+  delete[] u;
+  delete[] us;
+  delete[] un;
+  delete[] uns;
+  delete[] up;
+  delete[] ups;
+  delete[] unmp;
+  delete[] unsmp;
+  delete[] upmp;
+  delete[] upsmp;
+  delete[] argi;
+  delete[] args;
+  delete[] duk;
+  for (int ci = 3; ci > -1; ci--) {
+    delete[] cextlr[ci];
+    delete[] cext[ci];
+    delete[] cmullr[ci];
+    delete[] cmul[ci];
+  }
+  delete[] cextlr;
+  delete[] cext;
+  delete[] cmullr;
+  delete[] cmul;
+}
+
 ScatteringAngles::ScatteringAngles(GeometryConfiguration *gconf) {
   int isam = gconf->isam;
   _th = gconf->in_theta_start;

src/libnptm/clu_subs.cpp

@@ -9,6 +9,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif

src/libnptm/sph_subs.cpp

@@ -8,6 +8,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif

src/libnptm/tra_subs.cpp

@@ -11,6 +11,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif

src/trapping/cfrfme.cpp

@@ -18,6 +18,10 @@
 #include "../include/Parsers.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif

src/trapping/clffft.cpp

@@ -18,6 +18,10 @@
 #include "../include/Parsers.h"
 #endif
 
+#ifndef INCLUDE_LOGGING_H_
+#include "../include/logging.h"
+#endif
+
 #ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
 #endif