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 901f58dd2c081ee5e4f4923e0e73f3c266e53ff4..0000000000000000000000000000000000000000
--- a/src/scripts/generator.py
+++ /dev/null
@@ -1,105 +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_rand.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)]
-
-n_types = config['particle_settings']['n_types']
-if (len(config['particle_settings']['sph_types']) > 0):
- if (len(config['particle_settings']['sph_types']) != nsph):
- print("ERROR: declared number of types does not match the number of spheres!")
- exit(1)
-else:
- # Random type generation
- for ti in range(nsph):
- itype = 1 + int(n_types * random.random())
- config['particle_settings']['sph_types'].append(itype)
-sph_type_index = config['particle_settings']['sph_types'][0] - 1
-vec_spheres = [{'itype': sph_type_index + 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
-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_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 < 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_spheres.append({
- 'itype': sph_type_index + 1,
- 'x': x,
- 'y': y,
- 'z': z
- })
- is_placed = True
- placed_spheres += 1
- attempts = 0
-
-print(vec_spheres)
-print(sorted(vec_spheres, key=lambda item: item['itype']))
diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 70b4df231d689d9d03045c0c014e9a923a456023..555e3227a509c48015c5dcf358b075bd2e4bcc0f 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -1,4 +1,4 @@
-#!/bin/python3
+#!/usr/bin/env python3
# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
#
@@ -24,13 +24,15 @@
# The script requires python3.
import math
+import os
+import pdb
import random
import yaml
allow_3d = True
try:
import pyvista as pv
-except ModuleNotFound as ex:
+except ModuleNotFoundError as ex:
print("WARNING: pyvista not found!")
allow_3d = False
@@ -140,6 +142,7 @@ 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 if model['system_settings']['make_3D'] == "0" else True
if (make_3d and not allow_3d):
print("WARNING: 3D visualization of models is not available. Disabling.")
@@ -363,8 +366,8 @@ def load_model(model_file):
# Generate random cluster
rnd_seed = int(model['system_settings']['rnd_seed'])
# random_aggregate() checks internally whether application is INCLUSION
- print("DEBUG: make_3d is %s."%str(make_3d))
- random_aggregate(sconf, gconf, rnd_seed, 20.0 * max(sconf['ros']), make_3d)
+ max_rad = float(model['particle_settings']['max_rad'])
+ random_aggregate(sconf, gconf, rnd_seed, max_rad)
else:
if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
print("ERROR: coordinate vectors do not match the number of spheres!")
@@ -373,6 +376,11 @@ 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 (model['system_settings']['make_3D'] != "0" and allow_3d):
+ if (max_rad == 0.0):
+ max_rad = 20.0 * max(sconf['ros'])
+ write_obj(sconf, gconf, max_rad)
else: # model is None
print("ERROR: could not parse " + model_file + "!")
return (sconf, gconf)
@@ -500,8 +508,21 @@ def print_help():
print("--help Print this help and exit.")
print(" ")
-# random_aggregate()
-def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attempts=100):
+
+## \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
+def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
random.seed(seed)
nsph = scatterer['nsph']
vec_thetas = [0.0 for i in range(nsph)]
@@ -540,13 +561,14 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
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 < rr2):
+ 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])
@@ -560,7 +582,17 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
+ sinthi * sinphi * sinthj * sinphj
+ costhi * costhj
)
- rho += 2.0 * vec_rads[j] * cosalpha
+ 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])
@@ -586,17 +618,6 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
geometry['vec_sph_y'][sph_index] = sphere['y']
geometry['vec_sph_z'][sph_index] = sphere['z']
sph_index += 1
- print("DEBUG: make3d is %s"%str(make3d))
- if make3d:
- # Do something
- for vi in range(len(vec_rads)):
- radius = vec_rads[vi] / max_rad
- x = vec_spheres[vi]['x'] / max_rad
- y = vec_spheres[vi]['y'] / max_rad
- z = vec_spheres[vi]['z'] / max_rad
- mesh = pv.Sphere(radius, (x, y, z))
- mesh_name = "sphere_{0:04d}.obj".format(vi)
- mesh.save(mesh_name)
# end random_aggregate()
## \brief Write the geometry configuration dictionary to legacy format.
@@ -749,6 +770,67 @@ 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.
+#
+# The spherical units are saved as `sphere_XXXX.obj` files. The material
+# information is collected in a `model.mtl` file and the whole scene is written
+# to a `model.obj` file.
+#
+# \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):
+ color_strings = [
+ "1.0 1.0 1.0\n", # white
+ "1.0 0.0 0.0\n", # red
+ "0.0 0.0 1.0\n", # blue
+ "0.0 1.0 0.0\n", # green
+ ]
+ color_names = [
+ "white", "red", "blue", "green"
+ ]
+ mtl_file = open("model.mtl", "w")
+ for mi in range(len(color_strings)):
+ mtl_line = "newmtl mtl{0:d}\n".format(mi)
+ 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_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))
+ mesh.save("tmp_mesh.obj")
+ pl.add_mesh(mesh) #, color=color_by_name)
+ mesh_name = "sphere_{0:04d}.obj".format(si)
+ in_obj_file = open("tmp_mesh.obj", "r")
+ out_obj_file = open(mesh_name, "w")
+ in_line = in_obj_file.readline()
+ out_obj_file.write(in_line)
+ out_obj_file.write("mtllib model.mtl\n")
+ out_obj_file.write("usemtl mtl{0:d}\n".format(sph_type_index))
+ while (in_line != ""):
+ in_line = in_obj_file.readline()
+ out_obj_file.write(in_line)
+ in_obj_file.close()
+ out_obj_file.close()
+ pl.export_obj("model.obj")
+ os.remove("tmp_mesh.obj")
+
## \brief Exit code (0 for success)
exit_code = main()
exit(exit_code)