I am dealing with a box of spheres, all having the same radius and inserted into a finite box size without any overlaps. Along a given direction $\mathbf v$, I am trying to find:

  • (a) The two spheres $i,j$ whose difference vectors $\mathbf r= \mathbf r_i-\mathbf r_j$ is parallel to $\mathbf v$ and are furthest apart. That is, among the pairs of spheres that satisfy $\mathbf r \parallel \mathbf v,$ which pair has the largest $|r|.$

  • (b) Similarly, the two spheres furthest apart whose difference vector is now perpendicular to $\mathbf v.$ That is, among the pairs of spheres that satisfy $\mathbf r \perp \mathbf v,$ which pair has the largest $|r|.$

Regarding the thresholds with which we deem two normalized vectors $\mathbf u, \mathbf v$ to be parallel or perpendicular:

  • $\mathbf u \parallel \mathbf v$ iff $\mathbf u\cdot \mathbf v \ge 0.9,$ let's call these cases type=1,
  • $\mathbf u \perp \mathbf v$ iff $\mathbf u\cdot \mathbf v \le 0.1.$, let's call these cases type=-1,
  • and anything else we'll call type=0.


What I have managed to do so far is:

  • First create a list spherespairs of all possible pairs using Subset, and a n-by-n matrix matdists of all zeros to be filled with pairwise distances, where n is the number of spheres.
  • Defined a function furthest that takes two spheres, the target vector $\mathbf v$ (targetvec in the code), and computes their Euclidean distance (and updates matdists), their normalized difference vector and its inner product with $\mathbf v$ and the obtained type.
  • Then I apply furthest to all the possible pairs.
  • With the former, to find (a), I'd then first filter pairs whose type=1, then find the one with largest dist. Similarly for (b).

Working example and code:


spheres = {};
n = 200; (*number of spheres*)
r = 0.5; (*radius*)
boxlen = 20; (*cubic box length*)
targetvec = {0., 0., 1.}; (*this is our direction vector v*) 

Inserting the spheres randomly and without overlap:

While[Length[spheres] < n, s = RandomReal[{r, boxlen - r}, 3];
  If[And @@ (Norm[# - s] > 2*r & /@ spheres), AppendTo[spheres, s]]];

Visualisation of spheres and in red arrow the given vector $\mathbf v:$

cube = {Opacity[0.1], Cuboid[{0, 0, 0}, {boxlen, boxlen, boxlen}]};
  Sphere[#, r] & /@ spheres, {Red, Arrowheads[0.1], 
   Arrow[Tube[{{0, 0, 0}, targetvec*boxlen}, r]]}}, Boxed -> False]

List of all sphere pairs and initialization of matrix of distances:

spherepairs = Subsets[spheres, {2}];
matdists = ConstantArray[0, {n, n}];

The described function furthest:

furthest[sphere1_, sphere2_, targetvec_, boxsize_] := 
  Module[{r1 = sphere1, r2 = sphere2, v = targetvec, l = boxsize},
   r = Normalize[r1 - r2];
   innerprod = r.v;
   dist = Norm[r1 - r2];
   indexsphere1 = Flatten@Position[spheres, r1];
   indexsphere2 = Flatten@Position[spheres, r2];
   matdists[[indexsphere1[[1]], indexsphere2[[1]]]] = dist;
   matdists[[indexsphere2[[1]], indexsphere1[[1]]]] = dist;
   type = 0;
   Which[Abs[innerprod] >= 0.9,
    type = 1;,
    Abs[innerprod] <= 0.1,
    type = -1;,
    type = 0;
   {type, innerprod, dist}

Applying the function to all the pairs:

result = furthest[#1, #2, targetvec, boxlen] & @@@ 
  spherepairs // AbsoluteTiming

which takes about $6$ seconds for $n=200$ spheres!

Problem and question:

  • My approach seems to scale very inefficiently, to obtain the type, innerprod and dist of all possible pairs (all stored in result) takes about $6$ seconds for only $200$ particles, and I haven't even proceeded to filtering result to solve (a) and (b). Later I will have to do these calculations for systems of $n\approx 2000,$ so efficiency in finding (a) and (b) is of essence. I know my approach is really naive because I calculate everything for all pairs, as opposed to targetting my search to cases likely to be furthest apart. But I don't know how I could achieve such a targeted search! Any hints would be really helpful.

  • In my approach, is there any part where I am doing something completely inefficiently considering Mathematica's capabilities? In other words, is there a simple change from which my approach would benefit a major speed-up? I have a feeling my approach is really over-killing it...


1 Answer 1


Follows a script which determines with assigned accuracy the furthers parallel and perpendicular set of points.

distmin = 0;
distminp = 0;
e1 = 0.01;
e2 = 0.01;
e3 = 0.01;
For[i = 1, i <= n, i++, 
 For[j = i + 1, j <= n, 
  j++, {dx, dy, dz} = spheres[[i]] - spheres[[j]]; 
  If[Abs[dx/dz] < e1 && Abs[dy/dz] < e2, dist = dx^2 + dy^2 + dz^2; 
   If[dist > distmin, distmin = dist; imax = i; jmax = j]];
  l0 = dx^2 + dy^2;
  lpro = Sqrt[l0];
  If[Abs[dz/lpro] < e3, len = l0 + dz^2; 
   If[len > distminp, distminp = len; iperp = i; jperp = j]]

If[distmin == 0 || distminp == 0, Print["Tolerances must be loosened"],
  v1 = spheres[[imax]] - spheres[[jmax]];
  v2 = spheres[[iperp]] - spheres[[jperp]];
  gr1 = ParametricPlot3D[spheres[[imax]] mu + spheres[[jmax]] (1 - mu),{mu, 0, 1}, PlotStyle -> {Thick, Blue}];
  gr2 = ParametricPlot3D[spheres[[iperp]] mu + spheres[[jperp]] (1 - mu),{mu, 0, 1}, PlotStyle -> {Thick, Blue}];
  Show[gr0, gr1, gr2]

enter image description here


Given a direction vector $\vec v$ the points should be previously rotated such that the $z$ axis direction coincides with $\vec v$ direction.

  • $\begingroup$ Many thanks Cesareo! Regarding the cases where the direction vector is not along the z-axis: if I understood correctly the same approach holds if we first rotate everything such that $\vec{v}$ becomes the new z-axis, so I defined a rotation transform: targetvec = Normalize[{1.,1.,1.}], then rot = RotationTransform[{{0.,0.,1.},targetvec}] and rotated spheres according to spheres = MapThread[rot,{spheres}] and the cube: cube = GeometricTransformation[cube,rot]; But the results when drawn seem wrong, am I mis-implemented the transformation? Thanks again! $\endgroup$
    – user52181
    Feb 27, 2020 at 12:22
  • $\begingroup$ You can implement a rotation formula like Rodrigue's. See the reference en.wikipedia.org/wiki/Rodrigues%27_rotation_formula Here the rotation axis is obtained by Cross[v, {0,0,1}] $\endgroup$
    – Cesareo
    Feb 27, 2020 at 12:32
  • $\begingroup$ To use the built-in rotation functions, I did the following to implement your idea: angle = VectorAngle[{0.,0.,1.},targetvec] and the transform applied to all particles defined as rot = RotationTransform[angle,Cross[targetvec,{0.,0.,1.}]]; $\endgroup$
    – user52181
    Feb 27, 2020 at 14:47
  • $\begingroup$ Did it work now? $\endgroup$
    – Cesareo
    Feb 27, 2020 at 15:03
  • $\begingroup$ Sure! I've tested with targetvec=Normalize[{1.,1.,1.}] and it produces sensible results. But surprisingly, for greater n's, e.g. for $2000$ spheres (n=2000), the two For loops perform more slowly than my original approach, and the loops seem to be the only bottleneck (about 20 seconds absolute time) since the rotation transformations are very quick. $\endgroup$
    – user52181
    Feb 27, 2020 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.