diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index fcfd0a61c9a43ea44e2dc09f9f306ea66e7d9ae4..21a2eb1bb710bf56dd9f1a8271bba565d4937634 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -164,9 +164,10 @@ 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 (sconf['application'] == "INCLUSION"):
+            if (max_layers < sconf['nshl'][0] + 1):
+                max_layers = sconf['nshl'][0] + 1
         if (num_dielec != sconf['configurations']):
             print("ERROR: declared array of optical constants does not match configurations!")
             return (None, None)
@@ -176,8 +177,11 @@ def load_model(model_file):
                     [ 0 for j in range(max_layers)] for i in range(sconf['configurations'])
                 ]
                 for i in range(sconf['configurations']):
-                    if (len(model['particle_settings']['dielec_id'][i]) != 1 + int(sconf['nshl'][i] / 2)):
-                        print("ERROR: Declared materials in type %d do not match number of layers!"%(i + 1))
+                    expected_layer_num = 1 + int(sconf['nshl'][i] / 2)
+                    if (sconf['application'] == "INCLUSION" and i == 0):
+                        expected_layer_num += 1
+                    if (len(model['particle_settings']['dielec_id'][i]) != expected_layer_num):
+                        print("ERROR: Declared materials in type %d do not match the number of layers!"%(i + 1))
                         return (None, None)
                     else:
                         for j in range(1 + int(sconf['nshl'][i] / 2)):
@@ -268,11 +272,17 @@ def load_model(model_file):
                 [0.0 for j in range(max_layers)] for i in range(sconf['configurations'])
             ]
             for i in range(sconf['configurations']):
-                if (len(model['particle_settings']['rad_frac'][i]) != sconf['nshl'][i]):
-                    print("ERROR: Declared transition radii in type %d do not match number of layers!"%i)
+                expected_radii = sconf['nshl'][i]
+                if (sconf['application'] == "INCLUSION" and i == 0):
+                    expected_radii += 1
+                if (len(model['particle_settings']['rad_frac'][i]) != expected_radii):
+                    print("ERROR: Declared transition radii in type %d do not match the number of layers!"%(i + 1))
                     return (None, None)
                 else:
-                    for j in range(sconf['nshl'][i]):
+                    expected_radii = sconf['nshl'][i]
+                    if (sconf['application'] == "INCLUSION" and i == 0):
+                        expected_radii += 1
+                    for j in range(expected_radii):
                         sconf['rcf'][i][j] = float(model['particle_settings']['rad_frac'][i][j])
         # Create the gconf dict
         str_polar = model['radiation_settings']['polarization']
@@ -356,7 +366,10 @@ def match_grid(sconf):
     nxi = 0
     sconf['vec_xi'] = []
     for i in range(sconf['configurations']):
-        for j in range(sconf['nshl'][i]):
+        layers = sconf['nshl'][i]
+        if (sconf['application'] == "INCLUSION" and i == 0):
+            layers += 1
+        for j in range(layers):
             file_idx = sconf['dielec_id'][i][j]
             dielec_path = PurePath(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
             file_name = str(dielec_path)
@@ -377,6 +390,8 @@ def match_grid(sconf):
             if (max_layers == 0):
                 # This is executed only once
                 max_layers = max(sconf['nshl'])
+                if (sconf['application'] == "INCLUSION" and max_layers < sconf['nshl'][0] + 1):
+                    max_layers = sconf['nshl'][0] + 1
                 w_start = sconf['xi_start']
                 w_end = sconf['xi_end']
                 for wi in range(len(wavelengths)):
@@ -553,11 +568,11 @@ def write_legacy_sconf(conf):
         output.write("\n")
     for ci in range(conf['configurations']):
         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]):
+        if (conf['application'] == "INCLUSION" and ci == 0):
+            layers += 1
+        for cj in range(layers):
             str_line = " {0:.7E}\n".format(conf['rcf'][ci][cj])
             output.write(str_line)
     if (conf['application'] != "INCLUSION"):
@@ -585,35 +600,27 @@ def write_legacy_sconf(conf):
             # Write all layers in each configuration for every wavelength
             for xk in range(conf['nxi']):
                 for xi in range(conf['configurations']):
-                    layers = 1 + int(conf['nshl'][xi]) - 1
+                    layers = int(conf['nshl'][xi])
+                    if (xi == 0):
+                        layers += 1
                     for xj in range(layers):
                         rdc0 = conf['rdc0'][xj][xi][xk]
                         idc0 = conf['idc0'][xj][xi][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 (xi == 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 = 1 + int(conf['nshl'][xi]) - 1
+                layers = int(conf['nshl'][xi])
+                if (xi == 0):
+                    layers += 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)
     output.write("0\n")
     output.close()
     return result