Generalize scatterer configuration to multi-layer case

src/scripts/model_maker.py

@@ -18,7 +18,6 @@
 ## @package pycompare
 #  \brief Script to build models from YAML configuration files
 
-import cmath
 import numpy as np
 import yaml
 
@@ -93,12 +92,6 @@ def interpolate_constants(sconf):
                     ry0 = rpart[wi - 1]
                     iy0 = ipart[wi - 1]
                     dx = w - x0
-                    if (dci == 1):
-                        print("DEBUG: wi = %d"%wi)
-                        print("DEBUG: w = %e"%w)
-                        print("DEBUG: x0 = %e"%x0)
-                        print("DEBUG: x1 = %e"%x1)
-                        print("DEBUG: dx = %e"%dx)
                     dry = (ry1 - ry0) / (x1 - x0) * dx
                     diy = (iy1 - iy0) / (x1 - x0) * dx
                     ry = ry0 + dry
@@ -140,7 +133,7 @@ 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['exri'] = float(model['material_settings']['extern_refr'])
+        sconf['exdc'] = float(model['material_settings']['extern_diel'])
         sconf['wp'] = float(model['radiation_settings']['wp'])
         sconf['xip'] = float(model['radiation_settings']['xip'])
         sconf['idfc'] = int(model['material_settings']['diel_flag'])
@@ -148,30 +141,65 @@ def load_model(model_file):
         sconf['xi_start'] = float(model['radiation_settings']['scale_start'])
         sconf['xi_end'] = float(model['radiation_settings']['scale_end'])
         sconf['xi_step'] = float(model['radiation_settings']['scale_step'])
-        if (model['material_settings']['match_mode'] != "GRID"):
+        sconf['configurations'] = int(model['particle_settings']['n_types'])
+        sconf['dielec_file'] = model['material_settings']['dielec_file']
+        num_dielec = len(model['particle_settings']['dielec_id'])
+        if (len(model['particle_settings']['n_layers']) != sconf['configurations']):
+            print("ERROR: Declared number of layers does not match number of types!")
+            return (None, None)
+        else:
+            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])
+        max_layers = max(sconf['nshl'])
+        if (num_dielec != sconf['configurations']):
+            print("ERROR: declared array of optical constants does not match configurations!")
+            return (None, None)
+        else:
+            sconf['dielec_id'] = [
+                [ 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)
+                    return (None, None)
+                else:
+                    for j in range(sconf['nshl'][i]):
+                        sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
+        if (model['material_settings']['match_mode'] == "INTERPOLATE"):
             sconf['nxi'] = 1 + int((sconf['xi_end'] - sconf['xi_start']) / sconf['xi_step'])
             sconf['vec_xi'] = [0.0 for i in range(sconf['nxi'])]
             for i in range(sconf['nxi']):
                 sconf['vec_xi'][i] = sconf['xi_start'] + i * sconf['xi_step']
+            # Set up the dielectric constants
+            if (sconf['idfc'] == 0):
+                sconf['rdc0'] = [
+                    [
+                        [0.0 for k in range(sconf['nxi'])] for j in range(sconf['configurations'])
+                    ] for i in range(max_layers)
+                ]
+                sconf['idc0'] = [
+                    [
+                        [0.0 for k in range(sconf['nxi'])] for j in range(sconf['configurations'])
+                    ] for i in range(max_layers)
+                ]
+                interpolate_constants(sconf)
+        elif (model['material_settings']['match_mode'] == "GRID"):
+            match_grid(sconf)
+        else:
+            print("ERROR: %s is not a recognized match mode!"%(model['material_settings']['match_mode']))
+            return (None, None)
         if (model['radiation_settings']['scale_name'] == "WAVELENGTH"):
             sconf['insn'] = 3
-        sconf['configurations'] = int(model['particle_settings']['n_types'])
+        else:
+            print("ERROR: unsupported scaling mode (only WAVELENGTH available)!")
+            return (None, None)
         sph_types = model['particle_settings']['sph_types']
         if (len(sph_types) != sconf['nsph']):
             print("ERROR: vector of sphere types does not match the declared number of spheres!")
             return (None, None)
         else:
             sconf['vec_types'] = [int(str_typ) for str_typ in sph_types]
-        max_layers = 0
-        if (len(model['particle_settings']['n_layers']) != sconf['configurations']):
-            print("ERROR: Declared number of layers does not match number of types!")
-            return (None, None)
-        else:
-            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['nshl'][i] > max_layers):
-                    max_layers = sconf['nshl'][i]
         if (len(model['particle_settings']['radii']) != sconf['configurations']):
             print("ERROR: Declared number of radii does not match number of types!")
             return (None, None)
@@ -193,38 +221,84 @@ def load_model(model_file):
                 else:
                     for j in range(sconf['nshl'][i]):
                         sconf['rcf'][i][j] = float(model['particle_settings']['rad_frac'][i][j])
-        # Set up the dielectric constants
-        sconf['dielec_file'] = model['material_settings']['dielec_file']
-        num_dielec = len(model['particle_settings']['dielec_id'])
-        if (num_dielec != sconf['configurations']):
-            print("ERROR: declared array of optical constants does not match configurations!")
-            return (None, None)
-        else:
-            sconf['dielec_id'] = [
-                [ 0 for j in range(max_layers)] for i in range(sconf['configurations'])
-            ]
+    else: # model is None
+        print("ERROR: could not parse " + model_file + "!")
+    return (sconf, gconf)
+
+## \brief Populate the dielectric constant data matching a grid.
+#
+#  Important note: if the configuration requests that more than one
+#  optical constants file should be used, all the files must provide
+#  their constants for the same vector of wavelengths.
+#
+#  \param sconf: `dict` Scatterer configuration dictionary.
+#  \return result: `int` An exit code (0 if successful).
+def match_grid(sconf):
+    result = 0
+    max_layers = 0
+    nxi = 0
+    sconf['vec_xi'] = []
     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)
-                    return (None, None)
-                else:
         for j in range(sconf['nshl'][i]):
-                        sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
-        if (sconf['idfc'] == 0):
+            file_idx = sconf['dielec_id'][i][j]
+            file_name = sconf['dielec_file'][int(file_idx) - 1]
+            dielec_file = open(file_name, 'r')
+            wavelengths = []
+            rpart = []
+            ipart = []
+            str_line = dielec_file.readline()
+            while (str_line != ""):
+                if (not str_line.startswith('#')):
+                    split_line = str_line.split(',')
+                    if (len(split_line) == 3):
+                        wavelengths.append(float(split_line[0]))
+                        rpart.append(float(split_line[1]))
+                        ipart.append(float(split_line[2]))
+                str_line = dielec_file.readline()
+            dielec_file.close()
+            if (max_layers == 0):
+                # This is executed only once
+                max_layers = max(sconf['nshl'])
+                w_start = sconf['xi_start']
+                w_end = sconf['xi_end']
+                for wi in range(len(wavelengths)):
+                    w = wavelengths[wi]
+                    if (w >= w_start and w <= w_end):
+                        sconf['vec_xi'].append(w)
+                        nxi += 1
                 sconf['rdc0'] = [
                     [
-                    [0.0 for k in range(sconf['nxi'])] for j in range(sconf['nshl'][i])
-                ] for i in range(sconf['configurations'])
+                        [
+                            0.0 for dk in range(nxi)
+                        ] for dj in range(sconf['configurations'])
+                    ] for di in range(max_layers)
                 ]
                 sconf['idc0'] = [
                     [
-                    [0.0 for k in range(sconf['nxi'])] for j in range(sconf['nshl'][i])
-                ] for i in range(sconf['configurations'])
+                        [
+                            0.0 for dk in range(nxi)
+                        ] for dj in range(sconf['configurations'])
+                    ] for di in range(max_layers)
                 ]
-            interpolate_constants(sconf)
-    else: # model is None
-        print("ERROR: could not parse " + model_file + "!")
-    return (sconf, gconf)
+                sconf['nxi'] = nxi
+            # This is executed for all layers in all configurations
+            wi = 0
+            x = wavelengths[wi]
+            ry = rpart[wi]
+            iy = ipart[wi]
+            for dci in range(sconf['nxi']):
+                w = sconf['vec_xi'][dci]
+                while (w > x):
+                    x = wavelengths[wi]
+                    ry = rpart[wi]
+                    iy = ipart[wi]
+                    if (wi == len(wavelengths)):
+                        print("ERROR: file %s does not cover requested wavelengths!"%file_name)
+                        return 1
+                    wi += 1
+                sconf['rdc0'][j][i][dci] = ry
+                sconf['idc0'][j][i][dci] = iy
+    return result
 
 ## \brief Print a command-line help summary.
 def print_help():
@@ -246,7 +320,7 @@ def write_legacy_sconf(sconf):
     out_file = str(sconf['out_file'])
     nsph = sconf['nsph']
     ies = sconf['ies']
-    exri = sconf['exri']
+    exdc = sconf['exdc']
     wp = sconf['wp']
     xip = sconf['xip']
     idfc = sconf['idfc']
@@ -258,7 +332,7 @@ def write_legacy_sconf(sconf):
     str_line = " {0:3d}{1:3d}\n".format(nsph, ies)
     output.write(str_line)
     str_line = " {0:12.7E} {1:12.7E} {2:12.7E} {3:2d} {4:7d} {5:4d} {6:3d}\n".format(
-        exri, wp, xip, idfc, nxi, instpc, xi_flag
+        exdc, wp, xip, idfc, nxi, instpc, xi_flag
     )
     output.write(str_line)
     if (instpc == 0):