Update CLUSTER to use the new C1 constructor

src/cluster/cluster.cpp

@@ -113,16 +113,8 @@ void cluster(string config_file, string data_file, string output_path) {
     int lm = gconf->l_max;
     if (gconf->li > lm) lm = gconf->li;
     if (gconf->le > lm) lm = gconf->le;
-    C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
+    C1 *c1 = new C1(gconf, sconf);
     C3 *c3 = new C3();
-    for (int c1i = 0; c1i < nsph; c1i++) {
-      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(gconf->li, gconf->le, nsph);
     C1_AddOns *c1ao = new C1_AddOns(c4);
     // End of add-ons initialization

src/include/Commons.h

@@ -191,25 +191,25 @@ public:
  */
 class C4 {
 public:
-  //! \brief QUESTION: definition?
+  //! \brief LITPO = 2 * LI + 1.
   int litpo;
-  //! \brief QUESTION: definition?
+  //! \brief LITPOS = LITPO * LITPO
   int litpos;
-  //! \brief Maximum field expansion order plus one. QUESTION: correct?
+  //! \brief LMPO = LM + 1.
   int lmpo;
-  //! \brief Twice maximum field expansion order plus one. QUESTION: correct?
+  //! \brief LMTPO = 2 * LM + 1.
   int lmtpo;
-  //! \brief Square of `lmtpo`.
+  //! \brief LMTPOS = LMTPO * LMTPO.
   int lmtpos;
-  //! \brief QUESTION: definition?
+  //! \brief Internal field expansion order.
   int li;
   //! \brief QUESTION: definition?
   int nlim;
-  //! \brief QUESTION: definition?
+  //! \brief External field expansion order.
   int le;
   //! \brief QUESTION: definition?
   int nlem;
-  //! \brief Maximum field expansion order. QUESTION: correct?
+  //! \brief Maximum field expansion order.
   int lm;
   //! \brief Number of spheres.
   int nsph;

src/libnptm/Commons.cpp

@@ -26,6 +26,11 @@
 
 C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
   lm = (int)gconf->get_param("l_max");
+  int li = (int)gconf->get_param("li");
+  int le = (int)gconf->get_param("le");
+  if (lm == 0) {
+    lm = (li > le) ? li : le;
+  }
   nsph = (int)gconf->get_param("nsph");
   nlmmt = 2 * (lm * (lm + 2));
 
@@ -68,8 +73,15 @@ C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
   rxx = new double[nsph]();
   ryy = new double[nsph]();
   rzz = new double[nsph]();
-  ros = new double[nsph]();
-  for (int ri = 0; ri < nsph; ri++) {
+  if (nsph > 1) {
+    for (int c1i = 0; c1i < nsph; c1i++) {
+      rxx[c1i] = gconf->get_sph_x(c1i);
+      ryy[c1i] = gconf->get_sph_y(c1i);
+      rzz[c1i] = gconf->get_sph_z(c1i);
+    }
+  }
+  ros = new double[configurations]();
+  for (int ri = 0; ri < configurations; ri++) {
     ros[ri] = sconf->get_radius(ri);
   }