diff --git a/src/scripts/config_gen.yml b/src/scripts/config_gen.yml
deleted file mode 100644
index f1742a04c8a8e56f6faf64c5333c238b8a194fe7..0000000000000000000000000000000000000000
--- a/src/scripts/config_gen.yml
+++ /dev/null
@@ -1,93 +0,0 @@
-system_settings:
- max_host_ram : 0
- max_gpu_ram : 0
-
-input_settings:
- input_folder : "."
- spheres_file : "DEDFB"
- geometry_file: "DCLU"
-
-output_settings:
- output_folder: "."
- output_name : "c_OCLU"
- formats : [ "LEGACY", "HDF5" ]
- jwtm : 1
-
-particle_settings:
- application : "CLUSTER"
- n_spheres : 8
- n_types : 2
- sph_types : [ 1, 1, 1, 1, 2, 2, 2, 2 ]
- n_layers : [ 1, 1 ]
- radii : [ 6.000e-08, 4.000e-8 ]
- rad_frac : [ [ 1.0 ], [ 1.0 ] ]
- dielec_id : [ [ 1 ], [ 1 ] ]
-
-material_settings:
- diel_flag : 0
- extern_diel : 1.00e+00
- dielec_path : "../../ref_data"
- dielec_file : [ "eps_ashok_C.csv" ]
- dielec_fmt : [ "CSV" ]
- match_mode : "GRID"
- diel_const : [ ]
-
-radiation_settings:
- polarization: "LINEAR"
- scale_start : 1.00e-06
- scale_end : 2.00e-05
- scale_step : 5.00e-09
- wp : 2.99792e+08
- xip : 1.00e+00
- step_flag : 0
- scale_name : "WAVELENGTH"
-
-geometry_settings:
- li : 4
- le : 8
- npnt : 149
- npntts : 300
- iavm : 0
- isam : 0
- in_th_start : 0.0
- in_th_step : 0.0
- in_th_end : 0.0
- in_ph_start : 0.0
- in_ph_step : 0.0
- in_ph_end : 0.0
- sc_th_start : 0.0
- sc_th_step : 0.0
- sc_th_end : 0.0
- sc_ph_start : 0.0
- sc_ph_step : 0.0
- sc_ph_end : 0.0
- x_coords : [
- 8.00e-08,
- -8.00e-08,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00
- ]
- y_coords : [
- 0.00e+00,
- 0.00e+00,
- 8.00e-08,
- -8.00e-08,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00
- ]
- z_coords : [
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 0.00e+00,
- 8.00e-08,
- -8.00e-08,
- 16.00e-08,
- -16.00e-08
- ]
diff --git a/src/scripts/generator.py b/src/scripts/generator.py
deleted file mode 100644
index 97491a5f92f33e029b1691d70228fedd790c3cfd..0000000000000000000000000000000000000000
--- a/src/scripts/generator.py
+++ /dev/null
@@ -1,95 +0,0 @@
-import math
-import random
-import yaml
-
-#import pdb
-
-from pathlib import PurePath
-from sys import argv
-
-seed = int(argv[1])
-
-random.seed(seed)
-
-config = {}
-with open('config.yml', 'r') as stream:
- config = yaml.safe_load(stream)
-
-nsph = int(config['particle_settings']['n_spheres'])
-
-vec_thetas = [0.0 for i in range(nsph)]
-vec_phis = [0.0 for i in range(nsph)]
-vec_rads = [0.0 for i in range(nsph)]
-vec_sph_x = [0.0 for i in range(nsph)]
-vec_sph_y = [0.0 for i in range(nsph)]
-vec_sph_z = [0.0 for i in range(nsph)]
-
-sph_type_index = config['particle_settings']['sph_types'][0] - 1
-vec_rads[0] = config['particle_settings']['radii'][sph_type_index]
-max_rad = 20.0 * vec_rads[0]
-placed_spheres = 1
-attempts = 0
-max_attempts = 100
-for i in range(1, nsph):
- sph_type_index = config['particle_settings']['sph_types'][i] - 1
- vec_rads[i] = config['particle_settings']['radii'][sph_type_index]
- is_placed = False
- #breakpoint()
- while (not is_placed):
- if (attempts > max_attempts):
- print("WARNING: could not place sphere %d in allowed radius!"%i)
- break # while(not is_placed)
- vec_thetas[i] = math.pi * random.random()
- vec_phis[i] = 2.0 * math.pi * random.random()
- rho = vec_rads[0] + vec_rads[i]
- z = rho * math.cos(vec_thetas[i])
- y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
- x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
- j = 0
- while (j < i - 1):
- j += 1
- dx2 = (x - vec_sph_x[j]) * (x - vec_sph_x[j])
- dy2 = (y - vec_sph_y[j]) * (y - vec_sph_y[j])
- dz2 = (z - vec_sph_z[j]) * (z - vec_sph_z[j])
- dist2 = dx2 + dy2 + dz2
- rr2 = (vec_rads[i] + vec_rads[j]) * (vec_rads[i] + vec_rads[j])
- if (dist2 < rr2):
- # Spheres i and j are compenetrating.
- # Sphere i is moved out radially until it becomes externally
- # tangent to sphere j. Then the check is repeated, to verify
- # that no other sphere was penetrated. The process is iterated
- # until sphere i is placed or the maximum allowed radius is
- # reached.
- sinthi = math.sin(vec_thetas[i])
- sinthj = math.sin(vec_thetas[j])
- costhi = math.cos(vec_thetas[i])
- costhj = math.cos(vec_thetas[j])
- sinphi = math.sin(vec_phis[i])
- sinphj = math.sin(vec_phis[j])
- cosphi = math.cos(vec_phis[i])
- cosphj = math.cos(vec_phis[j])
- cosalpha = (
- sinthi * cosphi * sinthj * cosphj
- + sinthi * sinphi * sinthj * sinphj
- + costhi * costhj
- )
- rho += 2.0 * vec_rads[j] * cosalpha
- z = rho * math.cos(vec_thetas[i])
- y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
- x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
- j = 0
- continue # while(j < i - 1)
- if (rho + vec_rads[i] > max_rad):
- # The current direction is filled. Try another one.
- attempts += 1
- continue # while(not is_placed)
- vec_sph_x[i] = x
- vec_sph_y[i] = y
- vec_sph_z[i] = z
- is_placed = True
- placed_spheres += 1
- attempts = 0
-
-print(vec_sph_x)
-print(vec_sph_y)
-print(vec_sph_z)
diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 21a2eb1bb710bf56dd9f1a8271bba565d4937634..6e2d3d50e20f538a528ba64cf1850ebda9ab21b2 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -1,6 +1,6 @@
-#!/bin/python3
+#!/usr/bin/env python3
-# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
+# Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -23,9 +23,22 @@
#
# The script requires python3.
+import math
+import multiprocessing
+import numpy as np
+import os
+import pdb
+import random
import yaml
-from pathlib import PurePath
+allow_3d = True
+try:
+ import pyvista as pv
+except ModuleNotFoundError as ex:
+ print("WARNING: pyvista not found!")
+ allow_3d = False
+
+from pathlib import Path
from sys import argv
## \brief Main execution code
@@ -58,7 +71,7 @@ def interpolate_constants(sconf):
for i in range(sconf['configurations']):
for j in range(sconf['nshl'][i]):
file_idx = sconf['dielec_id'][i][j]
- dielec_path = PurePath(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
+ dielec_path = Path(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
file_name = str(dielec_path)
dielec_file = open(file_name, 'r')
wavelengths = []
@@ -131,9 +144,18 @@ def load_model(model_file):
except FileNotFoundError:
print("ERROR: " + model_file + " was not found!")
if model is not None:
+ max_rad = 0.0
+ make_3d = False
+ try:
+ make_3d = False if model['system_settings']['make_3D'] == "0" else True
+ except KeyError:
+ make_3d = False
+ if (make_3d and not allow_3d):
+ print("WARNING: 3D visualization of models is not available. Disabling.")
+ make_3d = False
# Create the sconf dict
sconf = {
- 'out_file': PurePath(
+ 'out_file': Path(
model['input_settings']['input_folder'],
model['input_settings']['spheres_file']
)
@@ -252,11 +274,21 @@ def load_model(model_file):
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'] and len(sph_types) != 0):
- print("ERROR: vector of sphere types does not match the declared number of spheres!")
- return (None, None)
- else:
+ if (len(sph_types) == sconf['nsph']):
sconf['vec_types'] = [int(str_typ) for str_typ in sph_types]
+ else:
+ if (len(sph_types) != 0):
+ print("ERROR: vector of sphere types does not match the declared number of spheres!")
+ return (None, None)
+ else: # len(sph_types) = 0
+ len_vec_x = len(model['geometry_settings']['x_coords'])
+ len_vec_y = len(model['geometry_settings']['y_coords'])
+ len_vec_z = len(model['geometry_settings']['z_coords'])
+ if (len_vec_x != 0 or len_vec_y != 0 or len_vec_z != 0):
+ print("ERROR: cannot assign random types with explicit sphere positions!")
+ return (None, None)
+ else:
+ sconf['vec_types'] = [0 for ti in range(sconf['nsph'])]
if (len(model['particle_settings']['radii']) != sconf['configurations']):
print("ERROR: Declared number of radii does not match number of types!")
return (None, None)
@@ -290,7 +322,7 @@ def load_model(model_file):
print("ERROR: %s is not a recognized polarization state."%str_polar)
return (None, None)
gconf = {
- 'out_file': PurePath(
+ 'out_file': Path(
model['input_settings']['input_folder'],
model['input_settings']['geometry_file']
)
@@ -324,22 +356,45 @@ def load_model(model_file):
gconf['vec_sph_y'] = [0.0 for i in range(gconf['nsph'])]
gconf['vec_sph_z'] = [0.0 for i in range(gconf['nsph'])]
if (gconf['application'] != "SPHERE" or gconf['nsph'] != 1):
- if (len(model['geometry_settings']['x_coords']) != len(model['geometry_settings']['y_coords'])):
+ len_vec_x = len(model['geometry_settings']['x_coords'])
+ len_vec_y = len(model['geometry_settings']['y_coords'])
+ len_vec_z = len(model['geometry_settings']['z_coords'])
+ if (len_vec_x != len_vec_y):
print("ERROR: X and Y coordinate vectors have different lengths!")
return (None, None)
- if (len(model['geometry_settings']['x_coords']) != len(model['geometry_settings']['z_coords'])):
+ if (len_vec_x != len_vec_z):
print("ERROR: X and Z coordinate vectors have different lengths!")
return (None, None)
- if (len(model['geometry_settings']['y_coords']) != len(model['geometry_settings']['z_coords'])):
+ if (len_vec_y != len_vec_z):
print("ERROR: Y and Z coordinate vectors have different lengths!")
return (None, None)
- if (len(model['particle_settings']['sph_types']) == 0 and len(model['geometry_settings']['x_coords']) != 0):
- print("ERROR: random type assignment is not allowed with assigned coordinates!")
- return (None, None)
- # TODO: put here a test to allow for empty vectors in case of random generation
- if (len(model['geometry_settings']['x_coords']) == 0):
+ if (len_vec_x == 0):
# Generate random cluster
- print("DEBUG: would generate random cluster.")
+ rnd_seed = int(model['system_settings']['rnd_seed'])
+ try:
+ max_rad = float(model['particle_settings']['max_rad'])
+ except KeyError:
+ print("ERROR: random model generation requires specification of particle_settings:max_rad.")
+ return (None, None)
+ rnd_engine = "COMPACT"
+ try:
+ rnd_engine = model['system_settings']['rnd_engine']
+ except KeyError:
+ # use compact generator, if no specification is given
+ rnd_engine = "COMPACT"
+ if (rnd_engine == "COMPACT"):
+ check = random_compact(sconf, gconf, rnd_seed, max_rad)
+ if (check == 1):
+ print("ERROR: compact random generator works only when all sphere types have the same radius.")
+ return (None, None)
+ elif (rnd_engine == "LOOSE"):
+ # random_aggregate() checks internally whether application is INCLUSION
+ check = random_aggregate(sconf, gconf, rnd_seed, max_rad)
+ else:
+ print("ERROR: unrecognized random generator engine.")
+ return (None, None)
+ if (check != 0):
+ print("WARNING: %d sphere(s) could not be placed."%check)
else:
if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
print("ERROR: coordinate vectors do not match the number of spheres!")
@@ -348,6 +403,40 @@ def load_model(model_file):
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])
+ if (make_3d and allow_3d):
+ if (max_rad == 0.0):
+ max_rad = 20.0 * max(sconf['ros'])
+ write_obj(sconf, gconf, max_rad)
+ try:
+ max_gpu_ram = int(model['system_settings']['max_gpu_ram'])
+ if (max_gpu_ram > 0):
+ max_gpu_ram_bytes = max_gpu_ram * 1024 * 1024 * 1024
+ matrix_dim = 2 * gconf['nsph'] * gconf['li'] * (gconf['li'] + 2)
+ matrix_size_bytes = 16 * matrix_dim * matrix_dim
+ if (matrix_size_bytes < max_gpu_ram_bytes):
+ max_gpu_processes = int(max_gpu_ram_bytes / matrix_size_bytes)
+ print("INFO: system supports up to %d simultaneous processes on GPU."%max_gpu_processes)
+ else:
+ print("WARNING: estimated matrix size is larger than available GPU memory!")
+ else:
+ print("INFO: no GPU RAM declared.")
+ max_host_ram = int(model['system_settings']['max_host_ram'])
+ if (max_host_ram > 0):
+ max_host_ram_bytes = max_host_ram * 1024 * 1024 * 1024
+ matrix_dim = 2 * gconf['nsph'] * gconf['li'] * (gconf['li'] + 2)
+ matrix_size_bytes = 16 * matrix_dim * matrix_dim
+ if (matrix_size_bytes < max_host_ram_bytes):
+ max_host_processes = int(max_host_ram_bytes / matrix_size_bytes / 2)
+ print("INFO: system supports up to %d simultaneous processes."%max_host_processes)
+ else:
+ print("WARNING: estimated matrix size is larger than available host memory!")
+ else:
+ print("WARNING: no host RAM declared!")
+ except KeyError as ex:
+ print(ex)
+ print("WARNING: missing system description! Cannot estimate recommended execution.")
+ cpu_count = multiprocessing.cpu_count()
+ print("INFO: the number of detected CPUs is %d."%cpu_count)
else: # model is None
print("ERROR: could not parse " + model_file + "!")
return (sconf, gconf)
@@ -371,7 +460,7 @@ def match_grid(sconf):
layers += 1
for j in range(layers):
file_idx = sconf['dielec_id'][i][j]
- dielec_path = PurePath(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
+ dielec_path = Path(sconf['dielec_path'], sconf['dielec_file'][int(file_idx) - 1])
file_name = str(dielec_path)
dielec_file = open(file_name, 'r')
wavelengths = []
@@ -475,6 +564,242 @@ def print_help():
print("--help Print this help and exit.")
print(" ")
+
+## \brief Generate a random cluster aggregate from YAML configuration options.
+#
+# This function generates a random aggregate of spheres using radial ejection
+# in random directions of new spheres until they become tangent to the
+# outermost sphere existing in that direction. The result of the generated
+# model is directly saved in the parameters of the scatterer and geometry
+# configuration dictionaries.
+#
+# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
+# \param geometry: `dict` Geometry configuration dictionary (gets modified)
+# \param seed: `int` Seed for the random sequence generation
+# \param max_rad: `float` Maximum allowed radial extension of the aggregate
+# \param max_attempts: `int` Maximum number of attempts to place a particle in any direction
+# \return result: `int` Function exit code (0 for success, otherwise number of
+# spheres that could not be placed)
+def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
+ result = 0
+ random.seed(seed)
+ nsph = scatterer['nsph']
+ vec_thetas = [0.0 for i in range(nsph)]
+ vec_phis = [0.0 for i in range(nsph)]
+ vec_rads = [0.0 for i in range(nsph)]
+ n_types = scatterer['configurations']
+ if (0 in scatterer['vec_types']):
+ tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
+ for ti in range(nsph):
+ itype = tincrement + int(n_types * random.random())
+ scatterer['vec_types'][ti] = itype
+ sph_type_index = scatterer['vec_types'][0] - 1
+ vec_spheres = [{'itype': sph_type_index + 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
+ vec_rads[0] = scatterer['ros'][sph_type_index]
+ placed_spheres = 1
+ attempts = 0
+ for i in range(1, nsph):
+ sph_type_index = scatterer['vec_types'][i] - 1
+ vec_rads[i] = scatterer['ros'][sph_type_index]
+ is_placed = False
+ while (not is_placed):
+ if (attempts > max_attempts):
+ result += 1
+ break # while(not is_placed)
+ vec_thetas[i] = math.pi * random.random()
+ vec_phis[i] = 2.0 * math.pi * random.random()
+ rho = vec_rads[0] + vec_rads[i]
+ z = rho * math.cos(vec_thetas[i])
+ y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
+ x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
+ j = 0
+ while (j < i - 1):
+ j += 1
+ dx2 = (x - vec_spheres[j]['x']) * (x - vec_spheres[j]['x'])
+ dy2 = (y - vec_spheres[j]['y']) * (y - vec_spheres[j]['y'])
+ dz2 = (z - vec_spheres[j]['z']) * (z - vec_spheres[j]['z'])
+ dist2 = dx2 + dy2 + dz2
+ rr2 = (vec_rads[i] + vec_rads[j]) * (vec_rads[i] + vec_rads[j])
+ if (dist2 < 0.9999 * rr2):
+ # Spheres i and j are compenetrating.
+ # Sphere i is moved out radially until it becomes externally
+ # tangent to sphere j. Then the check is repeated, to verify
+ # that no other sphere was penetrated. The process is iterated
+ # until sphere i is placed or the maximum allowed radius is
+ # reached.
+ # breakpoint()
+ sinthi = math.sin(vec_thetas[i])
+ sinthj = math.sin(vec_thetas[j])
+ costhi = math.cos(vec_thetas[i])
+ costhj = math.cos(vec_thetas[j])
+ sinphi = math.sin(vec_phis[i])
+ sinphj = math.sin(vec_phis[j])
+ cosphi = math.cos(vec_phis[i])
+ cosphj = math.cos(vec_phis[j])
+ cosalpha = (
+ sinthi * cosphi * sinthj * cosphj
+ + sinthi * sinphi * sinthj * sinphj
+ + costhi * costhj
+ )
+ D12 = math.sqrt(
+ vec_spheres[j]['x'] * vec_spheres[j]['x']
+ + vec_spheres[j]['y'] * vec_spheres[j]['y']
+ + vec_spheres[j]['z'] * vec_spheres[j]['z']
+ )
+ O1K = D12 * cosalpha
+ sinalpha = math.sqrt(1.0 - cosalpha * cosalpha)
+ sinbetaprime = D12 / (vec_rads[i] + vec_rads[j]) * sinalpha
+ cosbetaprime = math.sqrt(1.0 - sinbetaprime * sinbetaprime)
+ Op3K = (vec_rads[i] + vec_rads[j]) * cosbetaprime
+ rho = O1K + Op3K
+ z = rho * math.cos(vec_thetas[i])
+ y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
+ x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
+ j = 0
+ continue # while(j < i - 1)
+ if (rho + vec_rads[i] > max_rad):
+ # The current direction is filled. Try another one.
+ attempts += 1
+ continue # while(not is_placed)
+ vec_spheres.append({
+ 'itype': sph_type_index + 1,
+ 'x': x,
+ 'y': y,
+ 'z': z
+ })
+ is_placed = True
+ placed_spheres += 1
+ attempts = 0
+ sph_index = 0
+ for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
+ scatterer['vec_types'][sph_index] = sphere['itype']
+ geometry['vec_sph_x'][sph_index] = sphere['x']
+ geometry['vec_sph_y'][sph_index] = sphere['y']
+ geometry['vec_sph_z'][sph_index] = sphere['z']
+ sph_index += 1
+ return result
+
+## \brief Generate a random compact cluster from YAML configuration options.
+#
+# This function generates a random aggregate of spheres using the maximum
+# compactness packaging to fill a spherical volume with given maximum radius,
+# then it proceeds by subtracting random spheres from the outer layers, until
+# the aggregate is reduced to the desired number of spheres. The function
+# can only be used if all sphere types have the same radius. The result of the
+# generated model is directly saved in the parameters of the scatterer and
+# geometry configuration dictionaries.
+#
+# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
+# \param geometry: `dict` Geometry configuration dictionary (gets modified)
+# \param seed: `int` Seed for the random sequence generation
+# \param max_rad: `float` Maximum allowed radial extension of the aggregate
+# \return result: `int` Function exit code (0 for success, otherwise error code)
+def random_compact(scatterer, geometry, seed, max_rad):
+ result = 0
+ random.seed(seed)
+ nsph = scatterer['nsph']
+ n_types = scatterer['configurations']
+ if (0 in scatterer['vec_types']):
+ tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
+ for ti in range(nsph):
+ itype = tincrement + int(n_types * random.random())
+ scatterer['vec_types'][ti] = itype
+ if (max(scatterer['ros']) != min(scatterer['ros'])):
+ result = 1
+ else:
+ radius = scatterer['ros'][0]
+ x_centers = np.arange(-1.0 * max_rad + radius, max_rad, 2.0 * radius)
+ y_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
+ z_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
+ x_offset = radius
+ y_offset = radius
+ x_layer_offset = radius
+ y_layer_offset = radius / math.sqrt(3.0)
+ tmp_spheres = []
+ n_cells = len(x_centers) * len(y_centers) * len(z_centers)
+ print("INFO: the cubic space would contain %d spheres."%n_cells)
+ n_max_spheres = int((max_rad / radius) * (max_rad / radius) * (max_rad / radius) * 0.74)
+ print("INFO: the maximum radius allows for %d spheres."%n_max_spheres)
+ for zi in range(len(z_centers)):
+ if (x_layer_offset == 0.0):
+ x_layer_offset = radius
+ else:
+ x_layer_offset = 0.0
+ if (y_offset == 0.0):
+ y_offset = radius
+ else:
+ y_offset = 0.0
+ for yi in range(len(y_centers)):
+ if (x_offset == 0.0):
+ x_offset = radius
+ else:
+ x_offset = 0.0
+ for xi in range(len(x_centers)):
+ x = x_centers[xi] + x_offset + x_layer_offset
+ y = y_centers[yi] + y_offset
+ z = z_centers[zi]
+ extent = radius + math.sqrt(x * x + y * y + z * z)
+ if (extent < max_rad):
+ tmp_spheres.append({
+ 'itype': 1,
+ 'x': x,
+ 'y': y,
+ 'z': z
+ })
+ #tmp_spheres = [{'itype': 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
+ current_n = len(tmp_spheres)
+ print("INFO: before erosion there are %d spheres in use."%current_n)
+ rho = 10.0 * max_rad
+ discard_rad = 100.0 * max_rad
+ while (current_n > nsph):
+ theta = math.pi * random.random()
+ phi = 2.0 * math.pi * random.random()
+ x0 = rho * math.sin(theta) * math.cos(phi)
+ y0 = rho * math.sin(theta) * math.sin(phi)
+ z0 = rho * math.cos(theta)
+ closest_index = 0
+ minimum_distance = 1000.0 * max_rad
+ for di in range(len(tmp_spheres)):
+ x1 = tmp_spheres[di]['x']
+ if (x1 == discard_rad):
+ continue
+ y1 = tmp_spheres[di]['y']
+ z1 = tmp_spheres[di]['z']
+ distance = math.sqrt(
+ (x1 - x0) * (x1 - x0)
+ + (y1 - y0) * (y1 - y0)
+ + (z1 - z0) * (z1 - z0)
+ )
+ if (distance < minimum_distance):
+ closest_index = di
+ minimum_distance = distance
+ tmp_spheres[closest_index]['x'] = discard_rad
+ current_n -= 1
+ vec_spheres = []
+ sph_index = 0
+ for ti in range(len(tmp_spheres)):
+ sphere = tmp_spheres[ti]
+ if (sphere['x'] < max_rad):
+ sphere['itype'] = scatterer['vec_types'][sph_index]
+ sph_index += 1
+ vec_spheres.append(sphere)
+ #pl = pv.Plotter()
+ #for si in range(len(vec_spheres)):
+ # x = vec_spheres[si]['x'] / max_rad
+ # y = vec_spheres[si]['y'] / max_rad
+ # z = vec_spheres[si]['z'] / max_rad
+ # mesh = pv.Sphere(radius / max_rad, (x, y, z))
+ # pl.add_mesh(mesh)
+ #pl.export_obj("scene.obj")
+ sph_index = 0
+ for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
+ scatterer['vec_types'][sph_index] = sphere['itype']
+ geometry['vec_sph_x'][sph_index] = sphere['x']
+ geometry['vec_sph_y'][sph_index] = sphere['y']
+ geometry['vec_sph_z'][sph_index] = sphere['z']
+ sph_index += 1
+ return result
+
## \brief Write the geometry configuration dictionary to legacy format.
#
# \param conf: `dict` Geometry configuration dictionary.
@@ -625,6 +950,82 @@ def write_legacy_sconf(conf):
output.close()
return result
+## \brief Export the model to a set of OBJ files for 3D visualization.
+#
+# This function exports the model as a set of OBJ files (one for every
+# spherical unit, plus a single scene file) to allow for model visualization
+# with 3D software tools.
+#
+# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
+# \param geometry: `dict` Geometry configuration dictionary (gets modified)
+# \param max_rad: `float` Maximum allowed radial extension of the aggregate
+def write_obj(scatterer, geometry, max_rad):
+ out_dir = scatterer['out_file'].absolute().parent
+ out_model_path = Path(str(geometry['out_file']) + ".obj")
+ out_material_path = Path(str(geometry['out_file']) + ".mtl")
+ color_strings = [
+ "0.5 0.5 0.5\n", # white
+ "0.3 0.0 0.0\n", # red
+ "0.0 0.0 0.3\n", # blue
+ "0.0 0.3 0.0\n", # green
+ ]
+ color_names = [
+ "white", "red", "blue", "green"
+ ]
+ mtl_file = open(str(out_material_path), "w")
+ for mi in range(len(color_strings)):
+ mtl_line = "newmtl " + color_names[mi] + "\n"
+ mtl_file.write(mtl_line)
+ color_line = color_strings[mi]
+ mtl_file.write(" Ka " + color_line)
+ mtl_file.write(" Ks " + color_line)
+ mtl_file.write(" Kd " + color_line)
+ mtl_file.write(" Ns 100.0\n")
+ mtl_file.write(" Tr 0.0\n")
+ mtl_file.write(" illum 2\n\n")
+ mtl_file.close()
+ pl = pv.Plotter()
+ for si in range(scatterer['nsph']):
+ sph_type_index = scatterer['vec_types'][si]
+ # color_index = 1 + (sph_type_index % (len(color_strings) - 1))
+ # color_by_name = color_names[sph_type_index]
+ radius = scatterer['ros'][sph_type_index - 1] / max_rad
+ x = geometry['vec_sph_x'][si] / max_rad
+ y = geometry['vec_sph_y'][si] / max_rad
+ z = geometry['vec_sph_z'][si] / max_rad
+ mesh = pv.Sphere(radius, (x, y, z))
+ pl.add_mesh(mesh, color=None)
+ pl.export_obj(str(Path(str(out_dir), "TMP_MODEL.obj")))
+ tmp_model_file = open(str(Path(str(out_dir), "TMP_MODEL.obj")), "r")
+ out_model_file = open(str(out_model_path), "w")
+ mtl_line = "mtllib {0:s}\n".format(out_material_path.name)
+ sph_index = 0
+ sph_type_index = 0
+ old_sph_type_index = 0
+ str_line = tmp_model_file.readline()
+ while (str_line != ""):
+ if (str_line.startswith("mtllib")):
+ str_line = mtl_line
+ elif (str_line.startswith("g ")):
+ sph_index += 1
+ sph_type_index = scatterer['vec_types'][sph_index - 1]
+ if (sph_type_index == old_sph_type_index):
+ str_line = tmp_model_file.readline()
+ str_line = tmp_model_file.readline()
+ else:
+ old_sph_type_index = sph_type_index
+ color_index = sph_type_index % (len(color_names) - 1)
+ str_line = "g grp{0:04d}\n".format(sph_type_index)
+ out_model_file.write(str_line)
+ str_line = tmp_model_file.readline()
+ str_line = "usemtl {0:s}\n".format(color_names[color_index])
+ out_model_file.write(str_line)
+ str_line = tmp_model_file.readline()
+ out_model_file.close()
+ tmp_model_file.close()
+ os.remove(str(Path(str(out_dir), "TMP_MODEL.obj")))
+ os.remove(str(Path(str(out_dir), "TMP_MODEL.mtl")))
+
## \brief Exit code (0 for success)
exit_code = main()
exit(exit_code)
diff --git a/src/scripts/pycompare.py b/src/scripts/pycompare.py
index b4b0a0c39702116396203f7d401778ce134e28d1..778e40554789e73f1231893f9c4110fb0e4a75b3 100755
--- a/src/scripts/pycompare.py
+++ b/src/scripts/pycompare.py
@@ -1,6 +1,6 @@
-#!/bin/python3
+#!/usr/bin/env python3
-# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
+# Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
diff --git a/src/scripts/pydynrange.py b/src/scripts/pydynrange.py
index 491f68f80e83ac7542ce3b9eb873e1fbe813f6a2..5dd250a3b6f5dc0d6272b38b2c08c25873b8b9a0 100755
--- a/src/scripts/pydynrange.py
+++ b/src/scripts/pydynrange.py
@@ -1,6 +1,6 @@
-#!/bin/python3
+#!/usr/bin/env python3
-# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
+# Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
diff --git a/src/scripts/pytiming.py b/src/scripts/pytiming.py
index bc1ab97ae595ea8af25197e031f4befc50c38fce..c2d2e99e4acaf54250350612954df659f86a3ef9 100755
--- a/src/scripts/pytiming.py
+++ b/src/scripts/pytiming.py
@@ -1,6 +1,6 @@
-#!/bin/python3
+#!/usr/bin/env python3
-# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
+# Copyright (C) 2025 INAF - Osservatorio Astronomico di Cagliari
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
diff --git a/test_data/config_random_example.yml b/test_data/config_random_example.yml
new file mode 100644
index 0000000000000000000000000000000000000000..deb73c39dfb0ecda1c28f72c6253e16247001311
--- /dev/null
+++ b/test_data/config_random_example.yml
@@ -0,0 +1,140 @@
+# EXAMPLE CONFIGURATION FILE TO BUILD A RANDOM CLUSTER WITH 100 SPHERES.
+
+system_settings:
+ # Limit on host RAM use in Gb (0 for no configuration limit)
+ max_host_ram : 24
+ # Limit on GPU RAM use in Gb (0 for no configuration limit)
+ max_gpu_ram : 8
+ # Random seed (a positive integer number)
+ rnd_seed : 105
+ # Random engine (COMPACT or LOOSE)
+ rnd_engine : "COMPACT"
+ # OBJ model export flag (requires pyvista; 0 is FALSE)
+ make_3D : 0
+
+input_settings:
+ # Folder to write the code input configuration files
+ input_folder : "."
+ # Name of the scatterer description file
+ spheres_file : "DEDFB_comp"
+ # Name of the geometry description file
+ geometry_file: "DCLU_comp"
+
+output_settings:
+ # Folder for the code output storage
+ output_folder: "."
+ # Name of the main output file
+ output_name : "c_OCLU"
+ # Requested output formats
+ formats : [ "LEGACY", "HDF5" ]
+ # Index of the scale for transition matrix output
+ jwtm : 1
+
+particle_settings:
+ # What application to use (SPHERE | CLUSTER | INCLUSION)
+ application : "CLUSTER"
+ # Number of spheres
+ n_spheres : 100
+ # Number of sphere types
+ n_types : 2
+ # Vector of sphere type identifiers (what type is each sphere)
+ # For random genration it can be empty (random sphere types),
+ # or contain "n_spheres" elements (random position with assigned
+ # types).
+ sph_types : [ ]
+ # Vector of layers in types (how many layers in each type)
+ n_layers : [ 1, 1 ]
+ # Spherical monomer radii in m (one size for each type)
+ radii : [ 4.000e-8, 4.000e-8 ]
+ # Layer fractional radii (one per layer in each type)
+ rad_frac : [ [ 1.0 ], [ 1.0 ] ]
+ # Index of the dielectric constants (one per layer in each type)
+ #
+ # 1 is first file in `dielec_file`, 2 is second ...
+ dielec_id : [ [ 1 ], [ 1 ] ]
+ # Maximum radius for random particle size in m
+ max_rad : 3.0e-7
+
+material_settings:
+ # Dielectric properties definition (-1 = scaled, 0 = tabulated)
+ diel_flag : 0
+ # External medium dielectric constant
+ extern_diel : 1.00e+00
+ # Dielectric constant files folder
+ dielec_path : "../../ref_data"
+ # List of dielectric constant files (used if diel_flag = 0)
+ dielec_file : [ "eps_ashok_C.csv" ]
+ # Dielectric constant files format (same for all files)
+ dielec_fmt : [ "CSV" ]
+ # Matching method between optical constants and radiation wavelengths
+ #
+ # INTERPOLATE: the constants are interpolated on wavelengths
+ # GRID: only the wavelengths with defined constants are computed
+ #
+ match_mode : "GRID"
+ # Reference dielectric constants (used if diel_flag = -1)
+ #
+ # One real and one imaginary part for each layer in each type.
+ diel_const : [ ]
+
+radiation_settings:
+ # Radiation field polarization (LINEAR | CIRCULAR)
+ polarization: "LINEAR"
+ # First scale to be used
+ scale_start : 1.00e-06
+ # Last scale to be used
+ scale_end : 2.00e-05
+ # Calculation step (overridden if `match_mode` is GRID)
+ scale_step : 5.00e-09
+ # Peak Omega
+ wp : 2.99792e+08
+ # Peak scale
+ xip : 1.00e+00
+ # Define scale explicitly (0) or in equal steps (1)
+ step_flag : 0
+ # Type of scaling variable
+ scale_name : "WAVELENGTH"
+
+geometry_settings:
+ # Maximum internal field expansion
+ li : 6
+ # Maximum external field expansion (not used by SPHERE)
+ le : 8
+ # Number of transition layer integration points
+ npnt : 149
+ # Number of non transition layer integration points
+ npntts : 300
+ # Averaging mode
+ iavm : 0
+ # Meridional plane flag
+ isam : 0
+ # Starting incidence azimuth angle
+ in_th_start : 0.0
+ # Incidence azimuth angle incremental step
+ in_th_step : 0.0
+ # Ending incidence azimuth angle
+ in_th_end : 0.0
+ # Starting incidence elevation angle
+ in_ph_start : 0.0
+ # Incidence elevation angle incremental step
+ in_ph_step : 0.0
+ # Ending incidence elevation angle
+ in_ph_end : 0.0
+ # Starting scattered azimuth angle
+ sc_th_start : 0.0
+ # Scattered azimuth angle incremental step
+ sc_th_step : 0.0
+ # Ending scattered azimuth angle
+ sc_th_end : 0.0
+ # Starting scattered elevation angle
+ sc_ph_start : 0.0
+ # Scattered elevation angle incremental step
+ sc_ph_step : 0.0
+ # Ending scattered elevation angle
+ sc_ph_end : 0.0
+ # Vector of sphere X coordinates (one per sphere or empty for random)
+ x_coords : [ ]
+ # Vector of sphere Y coordinates (one per sphere or empty for random)
+ y_coords : [ ]
+ # Vector of sphere Z coordinates (one per sphere or empty for random)
+ z_coords : [ ]