diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 46c4433a84f0be4ccaae6df2f8518c5be53a4268..55a00d62b3460be4b50d9838ebf9546a4aa18c13 100755
--- a/src/scripts/model_maker.py
+++ b/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,39 +221,85 @@ 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'])
- ]
- 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 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']):
+ 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()
+ 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 dk in range(nxi)
+ ] for dj in range(sconf['configurations'])
+ ] for di in range(max_layers)
+ ]
+ sconf['idc0'] = [
+ [
+ [
+ 0.0 for dk in range(nxi)
+ ] for dj in range(sconf['configurations'])
+ ] for di in range(max_layers)
+ ]
+ 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():
print(" ")
@@ -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):