Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giacomo.mulas/np_tmcode
1 result
Select Git revision
Show changes
Commits on Source (14)
Show whitespace changes
Inline Side-by-side
src/scripts/config_gen.yml deleted 100644 → 0
View file @ 3fcc1996
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
]
src/scripts/generator.py deleted 100644 → 0
View file @ 3fcc1996
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)
src/scripts/model_maker.py
View file @ 25477404
#!/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):
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'] = [int(str_typ) for str_typ in sph_types]
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)
src/scripts/pycompare.py
View file @ 25477404
#!/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
......
src/scripts/pydynrange.py
View file @ 25477404
#!/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
......
src/scripts/pytiming.py
View file @ 25477404
#!/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
......
test_data/config_random_example.yml 0 → 100644
View file @ 25477404
# 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 : [ ]