The algorithm essentially starts with a set of spheres in random positions and in a relatively dense state, and runs an iterative process where the spheres are moved in each iteration in order to reduce the overlaps.<br>The core of the script is more or less like this<pre style="font-family:Monospaced;font-style:normal;font-size:12.0">
<span style="color:rgb(160,32,240)"></span>
superp<span style="color:rgb(74,85,219)"></span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(0,0,0)">calc_D_superp_dobcel</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">pos</span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span> //this line uses an external function to calculate a matrix with the overlap between paritcles, defined as superp(i,j)=max(R(i)+R(j)-D(i,j),0) . D(i,j) is the distance between the spheres and R is the radius. pos is a 3-colum matrix containging the coordinates of the position of each sphere<span style="color:rgb(160,32,240)"></span>
control<span style="color:rgb(92,92,92)">=</span><span style="color:rgb(188,143,143)">0</span><span style="color:rgb(0,0,0)">;</span><span style="color:rgb(0,0,0)"></span> //number of iterations
<span style="color:rgb(0,0,0)"></span>
<span style="color:rgb(0,0,0)">sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">int</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(0,0,0)">sumsupant</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(0,0,0)">;</span> // this defines a third control system as will be described bellow. ** <span style="color:rgb(0,0,0)"></span>
<span style="color:rgb(160,32,240)">while</span> <span style="color:rgb(50,185,185)">max</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">superp</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(0,0,0)">tol</span> // the maximum value of superp controls when the system is considered overlap-free and the algorithm is finished<br>
<br>...<br>// Several lines where the positions of each sphere are modified, acording to the values of the matrix superp<br>... <span style="color:rgb(100,174,100);font-style:italic"><br></span>
superp<span style="color:rgb(74,85,219)"></span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(0,0,0)">calc_D_superp_dobcel</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">pos</span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span> //matrix superp is re-calculated with the new positions<br>
<br>control=control+1;<br><br> <span style="color:rgb(160,32,240)">if</span> <span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(0,0,0)">max_mov</span> <span style="color:rgb(160,32,240)">then</span> <span style="color:rgb(95,158,160)">break</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(160,32,240)">end</span> // control by maximum number of iterations
<br>//the following belongs to the third control system:
<span style="color:rgb(160,32,240)">if</span> <span style="color:rgb(50,185,185)">length</span><span style="color:rgb(74,85,219)">(sup_2</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">==</span><span style="color:rgb(50,185,185)">int</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(160,32,240)">then</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">int</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(0,0,0)">sum(superp)</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(0,0,0)">;</span> //**
<span style="color:rgb(160,32,240)">else </span>//**
<span style="color:rgb(74,85,219)">sup_2</span><span style="color:rgb(74,85,219)"></span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">int</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">=sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">int</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">control</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">+</span> <span style="color:rgb(0,0,0)">sum(superp)</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">max_mov</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(188,143,143)">20</span><span style="color:rgb(0,0,0)">;</span> //**
<span style="color:rgb(160,32,240)">end</span> //**
<span style="color:rgb(160,32,240)">if</span> <span style="color:rgb(50,185,185)">length</span><span style="color:rgb(74,85,219)">(sup_2</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">4.5</span> <span style="color:rgb(160,32,240)">then</span> //**
<span style="color:rgb(160,32,240)">if</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">0.85</span><span style="color:rgb(92,92,92)">*sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">&</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">0.85</span><span style="color:rgb(92,92,92)">*sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">3</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(160,32,240)">then</span> <span style="color:rgb(95,158,160)">break</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(160,32,240)">end</span> //**
<span style="color:rgb(160,32,240)">if</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">0.7</span><span style="color:rgb(92,92,92)">*sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">&</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">0.7</span><span style="color:rgb(92,92,92)">*sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">3</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">&</span> sup_2<span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">3</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">></span><span style="color:rgb(188,143,143)">0.6</span><span style="color:rgb(92,92,92)">*sup_2</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(255,170,0)">$</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(188,143,143)">4</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(160,32,240)">then</span> <span style="color:rgb(95,158,160)">break</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(160,32,240)">end</span> //**
<span style="color:rgb(160,32,240)">end</span> //**
<span style="color:rgb(160,32,240)">end </span>//ends while loop</pre><br>The third control system essentially controls if the overlap is decreasing as the iterations proceed. It defines a new variable, sup_2, that stores the average value of sum(superp) in max_mov/20 iterations. If this average does not decrease fast enough (>15% decrease in two or >30%decrease in three consecutive groups of max_mov/20 iterations), it considers that it will not converge and breaks.<br>
Now, this third control system seems to be problem. As you can see, non of its lines (marked with ** for more clarity), affects the values of superp (at least, it shouldnīt). But still, the value of superp is modified when I include this control system.<br>
The other problem I had was with a similar algorithm, without the third control system, but just adding a "pause" rigth after the last end. This will also modify the values of superp and the number of iterations.<br>
As I said, the error in superp is not important (it is in the order of 1.d-15), but this difference propagates exponentially through the iterations and can significantly modify (like a 30%) the total number of iterations (I was trying to evaluate if the third control really worked, but I couldnīt do it because the algortihm CONVERGES with a different number of iterations...)<br>