Prepare a random generator script and a test configuration set

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

    ]

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)