Implement model making for CLUSTER development case

src/scripts/model_maker.py

@@ -36,10 +36,11 @@ from sys import argv
 # \returns result: `int` Number of detected error-level inconsistencies.
 def main():
     result = 0
-    if (len(argv) < 2):
+    config = parse_arguments()
+    if (config['help_mode'] or config['yml_file_name'] == ""):
         print_help()
     else:
-        sconf, gconf = load_model(argv[1])
+        sconf, gconf = load_model(config['yml_file_name'])
         if (sconf is not None) and (gconf is not None):
             result = write_legacy_sconf(sconf)
             if (result == 0): result += write_legacy_gconf(gconf)
@@ -150,6 +151,10 @@ def load_model(model_file):
         sconf['configurations'] = int(model['particle_settings']['n_types'])
         sconf['dielec_path'] = model['material_settings']['dielec_path']
         sconf['dielec_file'] = model['material_settings']['dielec_file']
+        num_dielec = 0 # len(model['particle_settings']['dielec_id'])
+        if (sconf['idfc'] == -1):
+            num_dielec = len(model['material_settings']['diel_const'])
+        elif (sconf['idfc'] == 0):
             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!")
@@ -173,6 +178,41 @@ def load_model(model_file):
                 else:
                     for j in range(sconf['nshl'][i]):
                         sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
+        if (model['radiation_settings']['scale_name'] == "XI"):
+            sconf['insn'] = 1
+            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 dielectric constants
+            if (num_dielec != sconf['configurations']):
+                print("ERROR: delcared dielectric constants do not match number of sphere types!")
+                return (None, None)
+            else:
+                if (sconf['idfc'] == -1):
+                    sconf['rdc0'] = [
+                        [
+                            [0.0 for k in range(1)] for j in range(sconf['configurations'])
+                        ] for i in range(max_layers)
+                    ]
+                    sconf['idc0'] = [
+                        [
+                            [0.0 for k in range(1)] for j in range(sconf['configurations'])
+                        ] for i in range(max_layers)
+                    ]
+                    for di in range(max_layers):
+                        for dj in range(sconf['configurations']):
+                            if (len(model['material_settings']['diel_const'][dj]) / 2 != sconf['nshl'][dj]):
+                                print("ERROR: dielectric constants for type %dj do not match number of layers!"%(dj + 1))
+                                return (None, None)
+                            else:
+                                sconf['rdc0'][di][dj][0] = float(model['material_settings']['diel_const'][dj][2 * di])
+                                sconf['idc0'][di][dj][0] = float(model['material_settings']['diel_const'][dj][2 * di + 1])
+                else: # sconf[idfc] != -1 and scaling on XI
+                    print("ERROR: for scaling on XI, optical constants must be defined!")
+                    return (None, None)
+        elif (model['radiation_settings']['scale_name'] == "WAVELENGTH"):
+            sconf['insn'] = 3
             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'])]
@@ -191,15 +231,16 @@ def load_model(model_file):
                         ] for i in range(max_layers)
                     ]
                     interpolate_constants(sconf)
+                else: # sconf[idfc] != 0 and scaling on wavelength
+                    print("ERROR: for wavelength scaling, optical constants must be tabulated!")
+                    return (None, None)
             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
         else:
-            print("ERROR: unsupported scaling mode (only WAVELENGTH available)!")
+            print("ERROR: %s is not a supported scaling mode!"%(model['radiation_settings']['scale_name']))
             return (None, None)
         sph_types = model['particle_settings']['sph_types']
         if (len(sph_types) != sconf['nsph']):
@@ -250,22 +291,37 @@ def load_model(model_file):
         gconf['npntts'] = int(model['geometry_settings']['npntts'])
         gconf['iavm'] = int(model['geometry_settings']['iavm'])
         gconf['isam'] = int(model['geometry_settings']['isam'])
+        gconf['jwtm'] = int(model['output_settings']['jwtm'])
         gconf['th'] = float(model['geometry_settings']['in_th_start'])
         gconf['thstp'] = float(model['geometry_settings']['in_th_step'])
         gconf['thlst'] = float(model['geometry_settings']['in_th_end'])
         gconf['ph'] = float(model['geometry_settings']['in_ph_start'])
         gconf['phstp'] = float(model['geometry_settings']['in_ph_step'])
         gconf['phlst'] = float(model['geometry_settings']['in_ph_end'])
-        gconf['ths'] = float(model['geometry_settings']['in_th_start'])
-        gconf['thsstp'] = float(model['geometry_settings']['in_th_step'])
-        gconf['thslst'] = float(model['geometry_settings']['in_th_end'])
-        gconf['phs'] = float(model['geometry_settings']['in_ph_start'])
-        gconf['phsstp'] = float(model['geometry_settings']['in_ph_step'])
-        gconf['phslst'] = float(model['geometry_settings']['in_ph_end'])
+        gconf['ths'] = float(model['geometry_settings']['sc_th_start'])
+        gconf['thsstp'] = float(model['geometry_settings']['sc_th_step'])
+        gconf['thslst'] = float(model['geometry_settings']['sc_th_end'])
+        gconf['phs'] = float(model['geometry_settings']['sc_ph_start'])
+        gconf['phsstp'] = float(model['geometry_settings']['sc_ph_step'])
+        gconf['phslst'] = float(model['geometry_settings']['sc_ph_end'])
+        if (gconf['application'] != "SPHERE" or gconf['nsph'] != 1):
+            # TODO: put here a test to allow for empty vectors in case of random generation
+            if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
+                print("ERROR: X coordinates do not match the number of spheres!")
+                return (None, None)
+            if (len(model['geometry_settings']['y_coords']) != gconf['nsph']):
+                print("ERROR: Y coordinates do not match the number of spheres!")
+                return (None, None)
+            if (len(model['geometry_settings']['z_coords']) != gconf['nsph']):
+                print("ERROR: Z coordinates do not match the number of spheres!")
+                return (None, None)
         gconf['vec_sph_x'] = [0.0 for i in range(gconf['nsph'])]
         gconf['vec_sph_y'] = [0.0 for i in range(gconf['nsph'])]
         gconf['vec_sph_z'] = [0.0 for i in range(gconf['nsph'])]
-        gconf['jwtm'] = int(model['output_settings']['jwtm'])
+        for si in range(gconf['nsph']):
+            gconf['vec_sph_x'][si] = float(model['geometry_settings']['x_coords'][si])
+            gconf['vec_sph_y'][si] = float(model['geometry_settings']['y_coords'][si])
+            gconf['vec_sph_z'][si] = float(model['geometry_settings']['z_coords'][si])
     else: # model is None
         print("ERROR: could not parse " + model_file + "!")
     return (sconf, gconf)
@@ -346,6 +402,30 @@ def match_grid(sconf):
                 sconf['idc0'][j][i][dci] = iy
     return result
 
+## \brief Parse the command line arguments.
+#
+#  The script behaviour can be modified through a set of optional arguments.
+#  The purpose of this function is to parse the command line in search for
+#  such arguments and prepare the execution accordingly.
+#
+#  \returns config: `dict` A dictionary containing the script configuration.
+def parse_arguments():
+    config = {
+        'yml_file_name': "",
+        'help_mode': False,
+    }
+    yml_set = False
+    for arg in argv[1:]:
+        if (arg.startswith("--help")):
+            config['help_mode'] = True
+        elif (not yml_set):
+            if (not arg.startswith("--")):
+                config['yml_file_name'] = arg
+                yml_set = True
+        else:
+            raise Exception("Unrecognized argument \'{0:s}\'".format(arg))
+    return config
+
 ## \brief Print a command-line help summary.
 def print_help():
     print("###############################################")
@@ -382,7 +462,7 @@ def write_legacy_gconf(conf):
         output.write(str_line)
     else:
         mxndm = 2 * nsph * conf['li'] * (conf['li'] + 2)
-        str_line = " {0:4d} {1:4d} {2:4d} {3:4d} {4:4d} {5:4d} {6:4d} {7:4d}\n".format(
+        str_line = " {0:4d} {1:4d} {2:4d} {3:4d} {4:4d} {5:4d} {6:4d} {7:4d} {8:4d}\n".format(
             nsph, conf['li'], conf['le'], mxndm, conf['inpol'],
             conf['npnt'], conf['npntts'], conf['iavm'], conf['isam']
         )
@@ -462,7 +542,16 @@ def write_legacy_sconf(conf):
                     for xii in range(conf['configurations']):
                         rdc0 = conf['rdc0'][xj][xii][xk]
                         idc0 = conf['idc0'][xj][xii][xk]
-                        if (rdc0 != 0.0 and idc0 != 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")