diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index bd31e186678b3c14a8191df724e04e4be7ee2ab9..9705290efc8b30736c79c18bc1a62d4cbea5b304 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -139,7 +139,8 @@ def load_model(model_file):
             )
         }
         sconf['nsph'] = int(model['particle_settings']['n_spheres'])
-        sconf['ies'] = 1 if model['particle_settings']['application'] == "INCLUSION" else 0
+        sconf['application'] = model['particle_settings']['application']
+        sconf['ies'] = 1 if sconf['application'] == "INCLUSION" else 0
         sconf['exdc'] = float(model['material_settings']['extern_diel'])
         sconf['wp'] = float(model['radiation_settings']['wp'])
         sconf['xip'] = float(model['radiation_settings']['xip'])
@@ -163,6 +164,8 @@ def load_model(model_file):
             sconf['nshl'] = [0 for i in range(sconf['configurations'])]
             for i in range(sconf['configurations']):
                 sconf['nshl'][i] = int(model['particle_settings']['n_layers'][i])
+            if (sconf['application'] == "INCLUSION"):
+                sconf['nshl'][0] += 1
         max_layers = max(sconf['nshl'])
         if (num_dielec != sconf['configurations']):
             print("ERROR: declared array of optical constants does not match configurations!")
@@ -531,31 +534,70 @@ def write_legacy_sconf(conf):
     if (sphere_count != 0):
         output.write("\n")
     for ci in range(conf['configurations']):
-        str_line = "{0:3d}   {1:15.7E}\n".format(conf['nshl'][ci], conf['ros'][ci])
+        layers = conf['nshl'][ci]
+        if (conf['application'] == "INCLUSION" and ci == 0):
+            layers -= 1
+        str_line = "{0:3d}   {1:15.7E}\n".format(layers, conf['ros'][ci])
         output.write(str_line)
         for cj in range(conf['nshl'][ci]):
             str_line = " {0:.7E}\n".format(conf['rcf'][ci][cj])
             output.write(str_line)
-    if (conf['idfc'] == 0):
-        # Write all wavelength dependent constants for each layer in each configuration
-        for xi in range(conf['configurations']):
-            for xj in range(conf['nshl'][xi]):
-                for xk in range(conf['nxi']):
-                    for xii in range(conf['configurations']):
-                        rdc0 = conf['rdc0'][xj][xii][xk]
-                        idc0 = conf['idc0'][xj][xii][xk]
+    if (conf['application'] != "INCLUSION"):
+        if (conf['idfc'] == 0):
+            # Write all wavelength dependent constants for each layer in each configuration
+            for xi in range(conf['configurations']):
+                for xj in range(conf['nshl'][xi]):
+                    for xk in range(conf['nxi']):
+                        for xii in range(conf['configurations']):
+                            rdc0 = conf['rdc0'][xj][xii][xk]
+                            idc0 = conf['idc0'][xj][xii][xk]
+                            if (rdc0 != 0.0 or idc0 != 0.0):
+                                str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
+                                output.write(str_line)
+        elif (conf['idfc'] == -1):
+            # Write reference scale constants for each layer in each configuration
+            for xi in range(conf['configurations']):
+                for xj in range(conf['nshl'][xi]):
+                    rdc0 = conf['rdc0'][xj][xi][0]
+                    idc0 = conf['idc0'][xj][xi][0]
+                    if (rdc0 != 0.0 or idc0 != 0.0):
+                        str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
+                        output.write(str_line)
+    else: # specialized output for INCLUSION
+        if (conf['idfc'] == 0):
+            # Write all wavelength dependent constants for each layer in each configuration
+            for xi in range(conf['configurations']):
+                layers = conf['nshl'][xi] - 1
+                for xj in range(layers):
+                    for xk in range(conf['nxi']):
+                        for xii in range(conf['configurations']):
+                            rdc0 = conf['rdc0'][xj][xii][xk]
+                            idc0 = conf['idc0'][xj][xii][xk]
+                            if (rdc0 != 0.0 or idc0 != 0.0):
+                                str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
+                                output.write(str_line)
+                            if (xii == conf['configurations'] - 1):
+                                rdc0 = conf['rdc0'][layers][0][xk]
+                                idc0 = conf['idc0'][layers][0][xk]
+                                if (rdc0 != 0.0 or idc0 != 0.0):
+                                    str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
+                                    output.write(str_line)
+        elif (conf['idfc'] == -1):
+            # Write reference scale constants for each layer in each configuration
+            for xi in range(conf['configurations']):
+                layers = conf['nshl'][xi] - 1
+                for xj in range(layers):
+                    rdc0 = conf['rdc0'][xj][xi][0]
+                    idc0 = conf['idc0'][xj][xi][0]
+                    if (rdc0 != 0.0 or idc0 != 0.0):
+                        str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
+                        output.write(str_line)
+                    if (xi == conf['configurations'] - 1):
+                        rdc0 = conf['rdc0'][layers][0][0]
+                        idc0 = conf['idc0'][layers][0][0]
                         if (rdc0 != 0.0 or idc0 != 0.0):
                             str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
                             output.write(str_line)
-    elif (conf['idfc'] == -1):
-        # Write reference scale constants for each layer in each configuration
-        for xi in range(conf['configurations']):
-            for xj in range(conf['nshl'][xi]):
-                rdc0 = conf['rdc0'][xj][xi][0]
-                idc0 = conf['idc0'][xj][xi][0]
-                if (rdc0 != 0.0 or idc0 != 0.0):
-                    str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
-                    output.write(str_line)
     output.write("0\n")
     output.close()
     return result
diff --git a/test_data/inclusion/config_dev.yml b/test_data/inclusion/config_dev.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e3269ef130a8d8214b416bfff0a93a1cbe29705c
--- /dev/null
+++ b/test_data/inclusion/config_dev.yml
@@ -0,0 +1,135 @@
+system_settings:
+  # Limit on host RAM use in Gb (0 for no configuration limit)
+  max_host_ram : 0
+  # Limit on GPU RAM use in Gb ( 0 for no configuration limit)
+  max_gpu_ram  : 0
+
+input_settings:
+  # Folder to write the code input configuration files
+  input_folder : "."
+  # Name of the scatterer description file
+  spheres_file : "DEDFB"
+  # Name of the geometry description file
+  geometry_file: "DINCLU"
+  
+output_settings:
+  # Folder for the code output storage
+  output_folder: "."
+  # Name of the main output file
+  output_name  : "c_OCLU"
+  # Requested output formats
+  formats      : [ "LEGACY", "HDF5" ]
+  # Index of the scale for transition matrix output
+  jwtm         : 2
+
+particle_settings:
+  # What application to use (SPHERE | CLUSTER | INCLUSION)
+  application : "INCLUSION"
+  # Number of spheres
+  n_spheres   : 2
+  # Number of sphere types
+  n_types     : 2
+  # Vector of sphere type identifiers (what type is each sphere)
+  sph_types   : [ 1, 2 ]
+  # Vector of layers in types (how many layers in each type)
+  n_layers    : [ 1, 1 ]
+  # Spherical monomer radii in m (one size for each type)
+  radii       : [ 2.0e-8, 2.0e-8 ]
+  # Layer fractional radii (one per layer in each type)
+  rad_frac    : [ [ 1.0, 3.0 ], [ 1.0 ] ]
+  # Index of the dielectric constants (one per layer in each type)
+  #
+  # 1 is first file in `dielec_file`, 2 is second ...
+  dielec_id   : [ [ 1, 3 ], [ 2 ] ]
+
+material_settings:
+  diel_flag   : 0
+  # External medium dielectric constant
+  extern_diel : 1.0e+00
+  # Dielectric constant files folder
+  dielec_path : "../../test_data/inclusion"
+  # List of dielectric constant files (used if diel_flag = 0)
+  dielec_file : [ "eps_sp1.csv", "eps_sp2.csv", "eps_spe.csv" ]
+  # Dielectric constant files format (same for all files)
+  dielec_fmt  : [ "CSV" ]
+  # Matching method between optical constants and radiation wavelengths
+  #
+  # INTERPOLATE: the constants are interpolated on wavelengths
+  # GRID: only the wavelengths with defined constants are computed
+  #
+  match_mode  : "GRID"
+  # Reference dielectric constants (used if diel_flag = -1)
+  #
+  # One real and one imaginary part for each layer in each type.
+  diel_const  : [ ]
+
+radiation_settings:
+  # Radiation field polarization (LINEAR | CIRCULAR)
+  polarization: "LINEAR"
+  # First scale to be used
+  scale_start : 1.0e-07
+  # Last scale to be used
+  scale_end   : 1.5e-07
+  # Calculation step (overridden if `match_mode` is GRID)
+  scale_step  : 5.0e-08
+  # Peak Omega
+  wp          : 3.00e+08
+  # Peak scale
+  xip         : 1.00e+00
+  # Define scale explicitly (0) or in equal steps (1)
+  step_flag   : 0
+  # Type of scaling variable
+  scale_name  : "WAVELENGTH"
+
+geometry_settings:
+  # Maximum internal field expansion
+  li          : 3
+  # Maximum external field expansion (not used by SPHERE)
+  le          : 3
+  # Number of transition layer integration points
+  npnt        : 149
+  # Number of non transition layer integration points
+  npntts      : 300
+  # Averaging mode
+  iavm        : 0
+  # Meridional plane flag
+  isam        : 0
+  # Starting incidence azimuth angle
+  in_th_start : 0.0
+  # Incidence azimuth angle incremental step
+  in_th_step  : 0.0
+  # Ending incidence azimuth angle
+  in_th_end   : 0.0
+  # Starting incidence elevation angle
+  in_ph_start : 0.0
+  # Incidence elevation angle incremental step
+  in_ph_step  : 0.0
+  # Ending incidence elevation angle
+  in_ph_end   : 0.0
+  # Starting scattered azimuth angle
+  sc_th_start : 0.0
+  # Scattered azimuth angle incremental step
+  sc_th_step  : 0.0
+  # Ending scattered azimuth angle
+  sc_th_end   : 0.0
+  # Starting scattered elevation angle
+  sc_ph_start : 0.0
+  # Scattered elevation angle incremental step
+  sc_ph_step  : 0.0
+  # Ending scattered elevation angle
+  sc_ph_end   : 0.0
+  # Vector of sphere X coordinates (one per sphere)
+  x_coords    : [
+    2.10e-08,
+    -2.10e-08
+    ]
+  # Vector of sphere Y coordinates (one per sphere)
+  y_coords    : [
+    0.00e+00,
+    0.00e+00
+    ]
+  # Vector of sphere Z coordinates (one per sphere)
+  z_coords    : [
+    0.00e+00,
+    0.00e+00
+    ]
diff --git a/test_data/inclusion/eps_sp1.csv b/test_data/inclusion/eps_sp1.csv
new file mode 100644
index 0000000000000000000000000000000000000000..9c6fb94a9ad6d8ed3b40dac6cf101e7b139cb1e1
--- /dev/null
+++ b/test_data/inclusion/eps_sp1.csv
@@ -0,0 +1,7 @@
+# Dielectric constants of SPHERE 1 for INCLUSION
+#
+# Wavelengths are in meters.
+#
+#Wavelength,Re(EPS),Im(EPS)
+1.0000000e-07,2.560000e00,0.000000e-01
+1.5000000e-07,4.000000e00,0.000000e-01
diff --git a/test_data/inclusion/eps_sp2.csv b/test_data/inclusion/eps_sp2.csv
new file mode 100644
index 0000000000000000000000000000000000000000..348a540efa3fbc7fe76a8b5ca18aaa32af6cb13d
--- /dev/null
+++ b/test_data/inclusion/eps_sp2.csv
@@ -0,0 +1,7 @@
+# Dielectric constants of SPHERE 2 for INCLUSION
+#
+# Wavelengths are in meters.
+#
+#Wavelength,Re(EPS),Im(EPS)
+1.0000000e-07,1.690000e00,0.000000e-01
+1.5000000e-07,2.000000e00,0.000000e-01
diff --git a/test_data/inclusion/eps_spe.csv b/test_data/inclusion/eps_spe.csv
new file mode 100644
index 0000000000000000000000000000000000000000..6218cb420bfc6a2ccf49bac1952745d6212c707a
--- /dev/null
+++ b/test_data/inclusion/eps_spe.csv
@@ -0,0 +1,7 @@
+# Dielectric constants of EXTERNAL SPHERE for INCLUSION
+#
+# Wavelengths are in meters.
+#
+#Wavelength,Re(EPS),Im(EPS)
+1.0000000e-07,2.000000e00,0.010000e-01
+1.5000000e-07,3.200000e00,0.030000e-01