generator.py

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)