diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index cfb0e28e6ae3258d6e90324a06311bb8154c31f9..fcfd0a61c9a43ea44e2dc09f9f306ea66e7d9ae4 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -176,11 +176,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]) != sconf['nshl'][i]):
-                        print("ERROR: Declared materials in type %d do not match number of layers!"%i)
+                    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))
                         return (None, None)
                     else:
-                        for j in range(sconf['nshl'][i]):
+                        for j in range(1 + int(sconf['nshl'][i] / 2)):
                             sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
         if (model['radiation_settings']['scale_name'] == "XI"):
             sconf['insn'] = 1
@@ -189,7 +189,8 @@ def load_model(model_file):
             for i in range(sconf['nxi']):
                 sconf['vec_xi'][i] = sconf['xi_start'] + i * sconf['xi_step']
             # Set up dielectric constants
-            if (num_dielec != sconf['configurations']):
+            allow_dielec = True # TODO: define logic to check dielectric constants
+            if (not allow_dielec):
                 print("ERROR: delcared dielectric constants do not match number of sphere types!")
                 return (None, None)
             else:
@@ -561,20 +562,19 @@ def write_legacy_sconf(conf):
             output.write(str_line)
     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)
+            # Write all layers in each configuration for every wavelength
+            for xk in range(conf['nxi']):
+                for xi in range(conf['configurations']):
+                    for xj in range(1 + int(conf['nshl'][xi] / 2)):
+                        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)
         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]):
+                for xj in range(1 + int(conf['nshl'][xi] / 2)):
                     rdc0 = conf['rdc0'][xj][xi][0]
                     idc0 = conf['idc0'][xj][xi][0]
                     if (rdc0 != 0.0 or idc0 != 0.0):
@@ -582,27 +582,26 @@ def write_legacy_sconf(conf):
                         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]
+            # 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
+                    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)
-                            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
+                layers = 1 + int(conf['nshl'][xi]) - 1
                 for xj in range(layers):
                     rdc0 = conf['rdc0'][xj][xi][0]
                     idc0 = conf['idc0'][xj][xi][0]