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)