diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 99b2315ccf6b12e144f5f2525a5de3063f9ff6cd..46c4433a84f0be4ccaae6df2f8518c5be53a4268 100755
--- a/src/scripts/model_maker.py
+++ b/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,36 +241,62 @@ 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':
- # Write legacy output
- #print("DEBUG: writing to file.")
- output = open(out_file, 'w')
- str_line = " {0:3d}{1:3d}\n".format(nsph, ies)
+ nxi = sconf['nxi']
+ # Write legacy output
+ output = open(out_file, 'w')
+ 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
+ )
+ 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)
- 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
- )
+ 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)
- output.close()
+ 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
def write_gconf(model, out_file):