#
# FreeCAD Python Macro to generate a Gaussian random pattern of cylinders in a cuboid volume
#
# Adrian Bowyer
# reprapltd.com
#
# 23 November 2018
#
# Licence: GPL
#

import random
from FreeCAD import Base

# Parameters

radius = 0.5         # The radius of the cylinders
meanDensity = 0.0005 # The expected number of cylinder centres per unit volume (1 mm^3)
sdDensity = 0.0001   # The standard deviation in the number of cylinder centres per unit volume 
meanLength = 20      # The mean length of the cylinders (mm)
sdLength = 5         # The standard deviation of the cylinder lengths (mm)

# The box from the origin in which all the centres of the cylinders will lie

xBox = 100
yBox = 100
zBox = 100

#********************************************************************************

# Generate a point from the uniform distribution over the defined box

def RanBoxPoint():
 x = random.uniform(0, xBox)
 y = random.uniform(0, yBox)
 z = random.uniform(0, zBox)
 return Base.Vector(x,y,z)


# Wrapper function for the Gaussian distribution
# (prevents negative and zero values)

def RanPositiveGaussian(mean, sd):
 g = -1
 while (g <= 0.0):
  g = random.gauss(mean, sd)
 return g


# Generate a point from the uniform distribution
# on the surface of the unit ball by rejection sampling
# from the surrounding cube.

def RanDirection():
 r2 = 2
 while (r2 > 1):
  x = random.uniform(-1, 1)
  y = random.uniform(-1, 1)
  z = random.uniform(-1, 1)
  r2 = x*x + y*y + z*z
 v = Base.Vector(x,y,z)
 v = v.multiply(1.0/v.Length)
 return v


# Make a single cylinder
# This assumes that the axis is a unit vector

def Cylinder(r, l, centre, axis):
 c = Part.makeCylinder(r,l,centre,axis,360)
 c.translate(axis.multiply(-l/2))
 return c

# The pattern is the union of all the cylinders

volume = xBox*yBox*zBox
number = int(RanPositiveGaussian(volume*meanDensity, volume*sdDensity))
c = Cylinder(radius, RanPositiveGaussian(meanLength, sdLength),RanBoxPoint(), RanDirection())
for i in range(0,number):
 c = c.fuse(Cylinder(radius, RanPositiveGaussian(meanLength, sdLength), RanBoxPoint(), RanDirection()))
Part.show(c)