diff --git a/src/include/Commons.h b/src/include/Commons.h
index 34254bf3d189fc60719e8568d2eb911d81e82248..173c21d916bac072d6a6d497f221e8ffce2841a7 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -654,6 +654,8 @@ public:
 class ParticleDescriptor {
 protected:
   // >>> COMMON TO ALL DESCRIPTOR TYPES <<< //
+  //! \brief Sub-class identification code.
+  short _class_type;
   //! \brief Number of spheres composing the model.
   int _nsph;
   //! \brief Maximum internal field expansion order.
@@ -716,6 +718,17 @@ protected:
 
 public:
   // >>> COMMON TO ALL DESCRIPTOR TYPES <<< //
+  //! \brief Base sub-class identification code.
+  const static short BASE_TYPE = 0;
+  //! \brief Sphere sub-class identification code.
+  const static short SPHERE_TYPE = 1;
+  //! \brief Cluster sub-class identification code.
+  const static short CLUSTER_TYPE = 2;
+  //! \brief Inclusion sub-class identification code.
+  const static short INCLUSION_TYPE = 3;
+  
+  //! \brief Sub-class identification code.
+  const short& class_type = _class_type;
   //! \brief Read-only view of number of spheres composing the model.
   const int &nsph = _nsph;
   //! \brief Read-only view of maximum internal field expansion order.
@@ -840,6 +853,25 @@ public:
   //! \brief J-vector components index matrix.
   int **ind3j;
   // >>> END OF SECTION NEEDED BY CLUSTER AND INCLU <<< //
+
+  // >>> NEEDED BY INCLU <<< //
+  //! \brief TBD
+  dcomplex *rm0;
+  //! \brief TBD
+  dcomplex *re0;
+  //! \brief TBD
+  dcomplex *rmw;
+  //! \brief TBD
+  dcomplex *rew;
+  //! \brief TBD
+  dcomplex *tm;
+  //! \brief TBD
+  dcomplex *te;
+  //! \brief TBD
+  dcomplex *tm0;
+  //! \brief TBD
+  dcomplex *te0;
+  // >>> END OF SECTION NEEDED BY INCLU <<< //
   
   /*! \brief ParticleDescriptor instance constructor.
    *
@@ -874,10 +906,6 @@ public:
    */
   ParticleDescriptorCluster(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
 
-  /*! \brief ParticleDescriptorSphere instance destroyer.
-   */
-  ~ParticleDescriptorCluster();
-  
   /*! \brief Interface function to return the descriptor type as string.
    *
    * \return descriptor_type: `string` The descriptor type name.
@@ -885,6 +913,28 @@ public:
   std::string get_descriptor_type() override { return "cluster descriptor"; }
 };
 
+/*! \brief The data structure describing a particle model for a sphere with inclusions.
+ *
+ * This class is used to solve the problem of a spherical particle with a
+ * cluster of inclusions. It replaces the C1 and C1_AddOns class implementation
+ * of the C1 FORTRAN common block.
+ */
+class ParticleDescriptorInclusion: public ParticleDescriptor {
+public:
+  /*! \brief ParticleDescriptorInclusion instance constructor.
+   *
+   * \param gconf: `const GeometryConfiguration *` Pointer to GeometryConfiguration instance.
+   * \param sconf: `const ScattererConfiguration *` Pointer to ScattererConfiguration instance.
+   */
+  ParticleDescriptorInclusion(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
+
+  /*! \brief Interface function to return the descriptor type as string.
+   *
+   * \return descriptor_type: `string` The descriptor type name.
+   */
+  std::string get_descriptor_type() override { return "inclusion descriptor"; }
+};
+
 /*! \brief The data structure describing a spherical particle model.
  *
  * This class is used to solve the problem of a single spherical particle. It
@@ -900,10 +950,6 @@ public:
    */
   ParticleDescriptorSphere(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
 
-  /*! \brief ParticleDescriptorSphere instance destroyer.
-   */
-  ~ParticleDescriptorSphere();
-  
   /*! \brief Interface function to return the descriptor type as string.
    *
    * \return descriptor_type: `string` The descriptor type name.
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 5f4c73668f21c54716176b048271fc241a7849e5..05875b1fe8e056df953c8ed9b6b2b69f2a57fd87 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1519,6 +1519,7 @@ ClusterIterationData::~ClusterIterationData() {
 
 // >>> ParticleDescriptor class implementation. <<< //
 ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
+  _class_type = BASE_TYPE;
   _nsph = gconf->number_of_spheres;
   _li = (nsph == 1) ? gconf->l_max : gconf->li;
   _le = 0;
@@ -1555,7 +1556,7 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
     rc[rci] = vec_rc + last_layer_index;
     last_layer_index += sconf->get_nshl(rci);
     if (rci == 0 && sconf->use_external_sphere) last_layer_index++;
-    for (int rcj = cur_layer_index; rcj < last_layer_index; rcj++) {
+    for (int rcj = 0; rcj < last_layer_index - cur_layer_index; rcj++) {
       rc[rci][rcj] = sconf->get_rcf(rci, rcj);
     }
     cur_layer_index = last_layer_index;
@@ -1616,9 +1617,19 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
   scscm = NULL;
   ecscm = NULL;
   ind3j = NULL;
+  // >>> NEEDED BY INCLU <<<
+  rm0 = NULL;
+  re0 = NULL;
+  rmw = NULL;
+  rew = NULL;
+  tm = NULL;
+  te = NULL;
+  tm0 = NULL;
+  te0 = NULL;
 }
 
 ParticleDescriptor::~ParticleDescriptor() {
+  // Base class members, always destroyed.
   delete[] nshl;
   delete[] ros;
   delete[] iog;
@@ -1634,12 +1645,73 @@ ParticleDescriptor::~ParticleDescriptor() {
   delete[] vec_rei;
   delete[] rmi;
   delete[] vec_rmi;
+  // Inclusion class members, destroyed only if sub-class is INCLUSION
+  if (_class_type == INCLUSION_TYPE) {
+    delete[] rm0;
+    delete[] re0;
+    delete[] rmw;
+    delete[] rew;
+    delete[] tm;
+    delete[] te;
+    delete[] tm0;
+    delete[] te0;
+  }
+  // Inclusion/cluster class members, destroyed if sub-class is INCLUSION or CLUSTER
+  if (_class_type == INCLUSION_TYPE || _class_type == CLUSTER_TYPE) {
+    delete[] vec_am0m;
+    delete[] vec_fsac;
+    delete[] vec_sac;
+    delete[] vec_fsacm;
+    delete[] vec_ind3j;
+  
+    delete[] vh;
+    delete[] vj0;
+    delete[] vyhj;
+    delete[] vyj0;
+    delete[] am0m;
+    delete[] fsac;
+    delete[] fsacm;
+    delete[] vintm;
+    delete[] scscp;
+    delete[] ecscp;
+    delete[] scscpm;
+    delete[] ecscpm;
+    delete[] v3j0;
+    delete[] scsc;
+    delete[] ecsc;
+    delete[] scscm;
+    delete[] ecscm;
+    delete[] ind3j;
+  }
+  // Cluster/sphere class members, destroyed if sub-class is CLUSTER or SPHERE
+  if (_class_type == CLUSTER_TYPE || _class_type == SPHERE_TYPE) {
+    for (int vi = 0; vi < _nsph; vi++) {
+      delete[] sas[vi];
+    }
+    delete[] vints;
+    delete[] sas;
+    delete[] gcsv;
+    delete[] sqabs;
+    delete[] sqexs;
+    delete[] sqscs;
+    delete[] sabs;
+    delete[] sexs;
+    delete[] sscs;
+    delete[] fsas;
+    delete[] vec_vints;
+    delete[] vec_sas;
+  }
+  // Cluster class members, destroyed only if sub-class is CLUSTER
+  if (_class_type == CLUSTER_TYPE) {
+    delete[] vintt;
+  }
 }
 
 // >>> End of ParticleDescriptor class implementation. <<< //
 
 // >>> ParticleDescriptorCluster class implementation. <<< //
 ParticleDescriptorCluster::ParticleDescriptorCluster(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {
+  _class_type = CLUSTER_TYPE;
   // Needed by SPHERE and CLUSTER
   vec_sas = new dcomplex[4 * nsph]();
   vec_vints = new dcomplex[16 * nsph]();
@@ -1708,58 +1780,70 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(GeometryConfiguration *gcon
   for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
 }
 
-ParticleDescriptorCluster::~ParticleDescriptorCluster() {
-  // Needed by SPHERE and CLUSTER
-  for (int vi = 0; vi < nsph; vi++) {
-    delete[] sas[vi];
-  }
-  delete[] vints;
-  delete[] sas;
-  delete[] gcsv;
-  delete[] sqabs;
-  delete[] sqexs;
-  delete[] sqscs;
-  delete[] sabs;
-  delete[] sexs;
-  delete[] sscs;
-  delete[] fsas;
-  delete[] vec_vints;
-  delete[] vec_sas;
-  
-  // Needed by CLUSTER
-  delete[] vintt;
+// >>> End of ParticleDescriptorCluster class implementation. <<< //
 
+// >>> ParticleDescriptorInclusion class implementation. <<< //
+ParticleDescriptorInclusion::ParticleDescriptorInclusion(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {
+  _class_type = INCLUSION_TYPE;
   // Needed by CLUSTER and INCLU
-  delete[] vec_am0m;
-  delete[] vec_fsac;
-  delete[] vec_sac;
-  delete[] vec_fsacm;
-  delete[] vec_ind3j;
+  _le = gconf->le;
+  _lm = (_li > _le) ? _li : _le;
+  _nlim = _nsph * _li * (_li + 2);
+  _nlem = _nsph * _le * (_le + 2);
+  _nlemt = 2 * _nlem;
+  _ncou = _nsph * _nsph - 1;
+  _litpo = _li + _li + 1;
+  _litpos = _litpo * _litpo;
+  _lmtpo = _li + _le + 1;
+  _lmtpos = _lmtpo * _lmtpo;
+  _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
+  vec_am0m = new dcomplex[_nlemt * _nlemt]();
+  vec_fsac = new dcomplex[4]();
+  vec_sac = new dcomplex[4]();
+  vec_fsacm = new dcomplex[4]();
+  vec_ind3j = new int[(_lm + 1) * _lm]();
   
-  delete[] vh;
-  delete[] vj0;
-  delete[] vyhj;
-  delete[] vyj0;
-  delete[] am0m;
-  delete[] fsac;
-  delete[] fsacm;
-  delete[] vintm;
-  delete[] scscp;
-  delete[] ecscp;
-  delete[] scscpm;
-  delete[] ecscpm;
-  delete[] v3j0;
-  delete[] scsc;
-  delete[] ecsc;
-  delete[] scscm;
-  delete[] ecscm;
-  delete[] ind3j;
+  vh = new dcomplex[_ncou * _litpo]();
+  vj0 = new dcomplex[_nsph * _lmtpo]();
+  vyhj = new dcomplex[_ncou * _litpos]();
+  vyj0 = new dcomplex[_nsph * _lmtpos]();
+  am0m = new dcomplex*[_nlemt];
+  for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
+  fsac = new dcomplex*[2];
+  fsacm = new dcomplex*[2];
+  for (int fi = 0; fi < 2; fi++) {
+    fsac[fi] = vec_fsac + (fi * 2);
+    fsacm[fi] = vec_fsacm + (fi * 2);
+  }
+  vintm = new dcomplex[16]();
+  scscp = new dcomplex[2]();
+  ecscp = new dcomplex[2]();
+  scscpm = new dcomplex[2]();
+  ecscpm = new dcomplex[2]();
+  v3j0 = new double[_nv3j]();
+  scsc = new double[2]();
+  ecsc = new double[2]();
+  scscm = new double[2]();
+  ecscm = new double[2]();
+  ind3j = new int*[_lm + 1];
+  for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+
+  // Needed by INCLU
+  rm0 = new dcomplex[_le]();
+  re0 = new dcomplex[_le]();
+  rmw = new dcomplex[_le]();
+  rew = new dcomplex[_le]();
+  tm = new dcomplex[_le]();
+  te = new dcomplex[_le]();
+  tm0 = new dcomplex[_le]();
+  te0 = new dcomplex[_le]();
 }
 
-// >>> End of ParticleDescriptorCluster class implementation. <<< //
+// >>> End of ParticleDescriptorInclusion class implementation. <<< //
 
 // >>> ParticleDescriptorSphere class implementation. <<< //
 ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {
+  _class_type = SPHERE_TYPE;
   vec_sas = new dcomplex[4 * nsph]();
   vec_vints = new dcomplex[16 * nsph]();
 
@@ -1781,23 +1865,6 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf,
   }
 }
 
-ParticleDescriptorSphere::~ParticleDescriptorSphere() {
-  for (int vi = 0; vi < nsph; vi++) {
-    delete[] sas[vi];
-  }
-  delete[] vints;
-  delete[] sas;
-  delete[] gcsv;
-  delete[] sqabs;
-  delete[] sqexs;
-  delete[] sqscs;
-  delete[] sabs;
-  delete[] sexs;
-  delete[] sscs;
-  delete[] fsas;
-  delete[] vec_vints;
-  delete[] vec_sas;
-}
 // >>> End of ParticleDescriptorSphere class implementation. <<< //
 
 // >>> ScatteringAngles class implementation. <<< //