From ed31ec867514d0b2a19d9c5a0b9ac058c5e26c3e Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Wed, 5 Mar 2025 17:24:05 +0100
Subject: [PATCH] Implement geometry configuration from script

---
 src/scripts/model_maker.py | 153 +++++++++++++++++++++++++++++--------
 1 file changed, 121 insertions(+), 32 deletions(-)

diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 55a00d62..d7e031e9 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -36,8 +36,9 @@ def main():
         print_help()
     else:
         sconf, gconf = load_model(argv[1])
-        if sconf is not None:
+        if (sconf is not None) and (gconf is not None):
             result = write_legacy_sconf(sconf)
+            if (result == 0): result += write_legacy_gconf(gconf)
         else:
             print("ERROR: could not create configuration.")
             result = 1
@@ -52,7 +53,8 @@ def interpolate_constants(sconf):
     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_path = PurePath(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
+            file_name = str(dielec_path)
             dielec_file = open(file_name, 'r')
             wavelengths = []
             rpart = []
@@ -142,6 +144,7 @@ def load_model(model_file):
         sconf['xi_end'] = float(model['radiation_settings']['scale_end'])
         sconf['xi_step'] = float(model['radiation_settings']['scale_step'])
         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 = len(model['particle_settings']['dielec_id'])
         if (len(model['particle_settings']['n_layers']) != sconf['configurations']):
@@ -221,6 +224,44 @@ 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])
+        # Create the gconf dict
+        str_polar = model['radiation_settings']['polarization']
+        if (str_polar not in ["LINEAR", "CIRCULAR"]):
+            print("ERROR: %s is not a recognized polarization state."%str_polar)
+            return (None, None)
+        gconf = {
+            'out_file': PurePath(
+                model['input_settings']['input_folder'],
+                model['input_settings']['geometry_file']
+            )
+        }
+        gconf['nsph'] = sconf['nsph']
+        gconf['application'] = model['particle_settings']['application']
+        gconf['li'] = int(model['geometry_settings']['li'])
+        gconf['le'] = int(
+            gconf['li'] if gconf['application'] == "SPHERE" else model['geometry_settings']['le']
+        )
+        gconf['inpol'] = 0 if str_polar == "LINEAR" else 1
+        gconf['npnt'] = int(model['geometry_settings']['npnt'])
+        gconf['npntts'] = int(model['geometry_settings']['npntts'])
+        gconf['iavm'] = int(model['geometry_settings']['iavm'])
+        gconf['isam'] = int(model['geometry_settings']['isam'])
+        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['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'])
     else: # model is None
         print("ERROR: could not parse " + model_file + "!")
     return (sconf, gconf)
@@ -241,7 +282,8 @@ def match_grid(sconf):
     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_path = PurePath(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
+            file_name = str(dielec_path)
             dielec_file = open(file_name, 'r')
             wavelengths = []
             rpart = []
@@ -302,8 +344,11 @@ def match_grid(sconf):
 
 ## \brief Print a command-line help summary.
 def print_help():
-    print("                                               ")
-    print("***              MODEL_MAKER                ***")
+    print("###############################################")
+    print("#                                             #")
+    print("#           NPtm_code MODEL_MAKER             #")
+    print("#                                             #")
+    print("###############################################")
     print("                                               ")
     print("Create input files for FORTRAN and C++ code.   ")
     print("                                               ")
@@ -315,18 +360,65 @@ def print_help():
     print("--help                Print this help and exit.")
     print("                                               ")
 
-def write_legacy_sconf(sconf):
+## \brief Write the geometry configuration dictionary to legacy format.
+#
+#  \param conf: `dict` Geometry configuration dictionary.
+#  \return result: `int` An exit code (0 if successful).
+def write_legacy_gconf(conf):
     result = 0
-    out_file = str(sconf['out_file'])
-    nsph = sconf['nsph']
-    ies = sconf['ies']
-    exdc = sconf['exdc']
-    wp = sconf['wp']
-    xip = sconf['xip']
-    idfc = sconf['idfc']
-    instpc = sconf['instpc']
+    out_file = str(conf['out_file'])
+    nsph = conf['nsph']
+    str_line = "INIT"
+    # Write legacy output
+    output = open(out_file, 'w')
+    if (conf['application'] == "SPHERE"):
+        str_line = " {0:4d} {1:4d} {2:4d} {3:4d} {4:4d} {5:4d}\n".format(
+            nsph, conf['li'], conf['inpol'], conf['npnt'], conf['npntts'], conf['isam']
+        )
+        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(
+            nsph, conf['li'], conf['le'], mxndm, conf['inpol'],
+            conf['npnt'], conf['npntts'], conf['iavm'], conf['isam']
+        )
+        output.write(str_line)
+        for si in range(nsph):
+            str_line = " {0:15.8E} {1:15.8E} {2:15.8E}\n".format(
+                conf['vec_sph_x'][si], conf['vec_sph_y'][si], conf['vec_sph_z'][si]
+            )
+            output.write(str_line)
+    str_line = " {0:7.2E}  {1:7.2E}  {2:7.2E}  {3:7.2E}  {4:7.2E}  {5:7.2E}\n".format(
+        conf['th'], conf['thstp'], conf['thlst'],
+        conf['ph'], conf['phstp'], conf['phlst']
+    )
+    output.write(str_line)
+    str_line = " {0:7.2E}  {1:7.2E}  {2:7.2E}  {3:7.2E}  {4:7.2E}  {5:7.2E}\n".format(
+        conf['ths'], conf['thsstp'], conf['thslst'],
+        conf['phs'], conf['phsstp'], conf['phslst']
+    )
+    output.write(str_line)
+    str_line = " {0:d}\n0\n".format(conf['jwtm'])
+    output.write(str_line)
+    output.close()
+    return result
+
+## \brief Write the scatterer configuration dictionary to legacy format.
+#
+#  \param conf: `dict` Scatterer configuration dictionary.
+#  \return result: `int` An exit code (0 if successful).
+def write_legacy_sconf(conf):
+    result = 0
+    out_file = str(conf['out_file'])
+    nsph = conf['nsph']
+    ies = conf['ies']
+    exdc = conf['exdc']
+    wp = conf['wp']
+    xip = conf['xip']
+    idfc = conf['idfc']
+    instpc = conf['instpc']
     xi_flag = 3
-    nxi = sconf['nxi']
+    nxi = conf['nxi']
     # Write legacy output
     output = open(out_file, 'w')
     str_line = " {0:3d}{1:3d}\n".format(nsph, ies)
@@ -337,14 +429,14 @@ def write_legacy_sconf(sconf):
     output.write(str_line)
     if (instpc == 0):
         for ixi in range(nxi):
-            str_line = "{0:.3E}\n".format(sconf['vec_xi'][ixi])
+            str_line = "{0:.3E}\n".format(conf['vec_xi'][ixi])
             output.write(str_line)
     else:
-        str_line = "{0:.3E}  {1:.3E}\n".format(sconf['xi_start'], sconf['xi_step'])
+        str_line = "{0:.3E}  {1:.3E}\n".format(conf['xi_start'], conf['xi_step'])
         output.write(str_line)
     sphere_count = 0
     for si in range(nsph):
-        str_line = "{0:5d}".format(sconf['vec_types'][si])
+        str_line = "{0:5d}".format(conf['vec_types'][si])
         output.write(str_line)
         sphere_count += 1
         if (sphere_count == 16):
@@ -352,20 +444,20 @@ def write_legacy_sconf(sconf):
             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])
+    for ci in range(conf['configurations']):
+        str_line = "{0:3d}   {1:15.7E}\n".format(conf['nshl'][ci], conf['ros'][ci])
         output.write(str_line)
-        for cj in range(sconf['nshl'][ci]):
-            str_line = " {0:.7E}\n".format(sconf['rcf'][ci][cj])
+        for cj in range(conf['nshl'][ci]):
+            str_line = " {0:.7E}\n".format(conf['rcf'][ci][cj])
             output.write(str_line)
-    if (sconf['idfc'] == 0):
+    if (conf['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]
+        for xi in range(conf['configurations']):
+            for xj in range(conf['nshl'][xi]):
+                for xk in range(conf['nxi']):
+                    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):
                             str_line = "({0:11.5E},{1:11.5E})\n".format(rdc0, idc0)
                             output.write(str_line)
@@ -373,9 +465,6 @@ def write_legacy_sconf(sconf):
     output.close()
     return result
 
-def write_gconf(model, out_file):
-    return
-
 ## \brief Exit code (0 for success)
 exit_code = main()
 exit(exit_code)
-- 
GitLab