Implement scatterer configuration for single-layer single sphere

src/scripts/model_maker.py

@@ -19,8 +19,10 @@
 #  \brief Script to build models from YAML configuration files
 
 import cmath
+import numpy as np
 import yaml
 
+from pathlib import PurePath
 from sys import argv
 
 ## \brief Main execution code
@@ -34,25 +36,195 @@ def main():
     if (len(argv) < 2):
         print_help()
     else:
-        model = load_model(argv[1])
-        if model is not None:
-            result = write_sconf(model)
+        sconf, gconf = load_model(argv[1])
+        if sconf is not None:
+            result = write_legacy_sconf(sconf)
         else:
             print("ERROR: could not create configuration.")
             result = 1
     return result
 
+## \brief Populate the dielectric constant data via interpolation
+#
+#  \param sconf: `dict` Scatterer configuration dictionary.
+#  \return result: `int` An exit code (0 if successful).
+def interpolate_constants(sconf):
+    result = 0
+    for i in range(sconf['configurations']):
+        for j in range(sconf['nshl'][i]):
+            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()
+            wi = 0
+            x0 = 0.0
+            x1 = 0.0
+            ry0 = 0.0
+            iy0 = 0.0
+            ry1 = 0.0
+            iy1 = 0.0
+            for dci in range(sconf['nxi']):
+                w = sconf['vec_xi'][dci]
+                while (w > x1):
+                    x0 = wavelengths[wi]
+                    ry0 = rpart[wi]
+                    iy0 = ipart[wi]
+                    if (wi == len(wavelengths)):
+                        print("ERROR: file %s does not cover requested wavelengths!"%file_name)
+                        return 1
+                    wi += 1
+                    x1 = wavelengths[wi]
+                    ry1 = rpart[wi]
+                    iy1 = ipart[wi]
+                if (wi > 0):
+                    x0 = wavelengths[wi - 1]
+                    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
+                    iy = iy0 + diy
+                    sconf['rdc0'][j][i][dci] = ry
+                    sconf['idc0'][j][i][dci] = iy
+                else:
+                    if (wavelengths[wi] == w):
+                        sconf['rdc0'][j][i][dci] = rpart[0]
+                        sconf['idc0'][j][i][dci] = ipart[0]
+                    else:
+                        print("ERROR: file %s does not cover requested wavelengths!"%file_name)
+                        return 2
+    return result
+
+## \brief Create tha calculation configuration structure from YAML input.
+#
+#  \param model_file: `str` Full path to the YAML input file.
+#  \return sconf, gconf: `tuple` A dictionary tuple for scatterer and
+#                        geometric configurations.
 def load_model(model_file):
-    result = None
+    sconf = None
+    gconf = None
+    model = None
     try:
         with open(model_file, 'r') as stream:
-            result = yaml.safe_load(stream)
+            model = yaml.safe_load(stream)
     except yaml.YAMLError:
-        print("ERROR: " + model_file + " is not a valid configuration!")
-        result = None
+        print("ERROR: " + model_file + " is not a valid YAML file!")
     except FileNotFoundError:
         print("ERROR: " + model_file + " was not found!")
-    return result
+    if model is not None:
+        # Create the sconf dict
+        sconf = {
+            'out_file': PurePath(
+                model['input_settings']['input_folder'],
+                model['input_settings']['spheres_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['wp'] = float(model['radiation_settings']['wp'])
+        sconf['xip'] = float(model['radiation_settings']['xip'])
+        sconf['idfc'] = int(model['material_settings']['diel_flag'])
+        sconf['instpc'] = int(model['radiation_settings']['step_flag'])
+        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['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']
+        if (model['radiation_settings']['scale_name'] == "WAVELENGTH"):
+            sconf['insn'] = 3
+        sconf['configurations'] = int(model['particle_settings']['n_types'])
+        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)
+        else:
+            sconf['ros'] = [0.0 for i in range(sconf['configurations'])]
+            for i in range(sconf['configurations']):
+                sconf['ros'][i] = float(model['particle_settings']['radii'][i])
+        if (len(model['particle_settings']['rad_frac']) != sconf['configurations']):
+            print("ERROR: Declared number of fractional radii does not match number of types!")
+            return (None, None)
+        else:
+            sconf['rcf'] = [
+                [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)
+                    return (None, None)
+                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'])
+            ]
+            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):
+            sconf['rdc0'] = [
+                [
+                    [0.0 for k in range(sconf['nxi'])] for j in range(sconf['nshl'][i])
+                ] for i in range(sconf['configurations'])
+            ]
+            sconf['idc0'] = [
+                [
+                    [0.0 for k in range(sconf['nxi'])] for j in range(sconf['nshl'][i])
+                ] for i in range(sconf['configurations'])
+            ]
+            interpolate_constants(sconf)
+    else: # model is None
+        print("ERROR: could not parse " + model_file + "!")
+    return (sconf, gconf)
 
 ## \brief Print a command-line help summary.
 def print_help():
@@ -69,28 +241,19 @@ def print_help():
     print("--help                Print this help and exit.")
     print("                                               ")
 
-def write_sconf(model, form='legacy'):
+def write_legacy_sconf(sconf):
     result = 0
-    out_file = (
-        model['input_settings']['input_folder'] + "/" +
-        model['input_settings']['spheres_file']
-    )
-    #print("DEBUG: out_file = " + out_file)
-    nsph = model['particle_settings']['n_spheres']
-    ies = 1 if model['particle_settings']['application'] == "inclusion" else 0
-    exri= float(model['material_settings']['extern_refr'])
-    wp = float(model['radiation_settings']['wp'])
-    xip = float(model['radiation_settings']['xip'])
-    scale_start = float(model['radiation_settings']['scale_start'])
-    scale_end = float(model['radiation_settings']['scale_end'])
-    scale_step = float(model['radiation_settings']['scale_step'])
-    idfc = int(model['radiation_settings']['diel_flag'])
-    instpc = int(model['radiation_settings']['step_flag'])
+    out_file = str(sconf['out_file'])
+    nsph = sconf['nsph']
+    ies = sconf['ies']
+    exri = sconf['exri']
+    wp = sconf['wp']
+    xip = sconf['xip']
+    idfc = sconf['idfc']
+    instpc = sconf['instpc']
     xi_flag = 3
-    nxi = 1 + int((scale_end - scale_start) / scale_step)
-    if form == 'legacy':
+    nxi = sconf['nxi']
     # Write legacy output
-        #print("DEBUG: writing to file.")
     output = open(out_file, 'w')
     str_line = " {0:3d}{1:3d}\n".format(nsph, ies)
     output.write(str_line)
@@ -98,6 +261,41 @@ def write_sconf(model, form='legacy'):
         exri, wp, xip, idfc, nxi, instpc, xi_flag
     )
     output.write(str_line)
+    if (instpc == 0):
+        for ixi in range(nxi):
+            str_line = "{0:.3E}\n".format(sconf['vec_xi'][ixi])
+            output.write(str_line)
+    else:
+        str_line = "{0:.3E}  {1:.3E}\n".format(sconf['xi_start'], sconf['xi_step'])
+        output.write(str_line)
+    sphere_count = 0
+    for si in range(nsph):
+        str_line = "{0:5d}".format(sconf['vec_types'][si])
+        output.write(str_line)
+        sphere_count += 1
+        if (sphere_count == 16):
+            output.write("\n")
+            sphere_count = 0
+    if (sphere_count != 0):
+        output.write("\n")
+    for ci in range(sconf['configurations']):
+        str_line = "{0:3d}   {1:15.7E}\n".format(sconf['nshl'][ci], sconf['ros'][ci])
+        output.write(str_line)
+        for cj in range(sconf['nshl'][ci]):
+            str_line = " {0:.7E}\n".format(sconf['rcf'][ci][cj])
+            output.write(str_line)
+    if (sconf['idfc'] == 0):
+        # Write all wavelength dependent constants for each layer in each configuration
+        for xi in range(sconf['configurations']):
+            for xj in range(sconf['nshl'][xi]):
+                for xk in range(sconf['nxi']):
+                    for xii in range(sconf['configurations']):
+                        rdc0 = sconf['rdc0'][xj][xii][xk]
+                        idc0 = sconf['idc0'][xj][xii][xk]
+                        if (rdc0 != 0.0 and 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