cylaxis[t_, {r1_,t1_,p1_,a1_,b1_}] := {t Cos[a1] Cos[b1]+r1 Cos[p1] Cos[t1],t Cos[b1] Sin[a1]+r1 Cos[p1] Sin[t1],t Sin[b1]+r1 Sin[p1]}

axialdist[{r1_, t1_, p1_,a1_, b1_} ,{r2_, t2_, p2_,a2_, b2_}]:=Chop[√((((Cos[b1] Sin[b2] (-r1 Cos[p1] Sin[a1-t1]+r2 Cos[p2] Sin[a1-t2])+Cos[b2] (Cos[a2] Cos[b1] Sin[a1] (r1 Sin[p1]-r2 Sin[p2])-Cos[a1] Cos[b1] Sin[a2] (r1 Sin[p1]-r2 Sin[p2])+Sin[b1] (r1 Cos[p1] Sin[a2-t1]-r2 Cos[p2] Sin[a2-t2]))))^2)/(((Cos[a2])^2) ((Cos[b2])^2) (((Cos[b1])^2) ((Sin[a1])^2)+((Sin[b1])^2))+(((Cos[b2] Sin[a2] Sin[b1]-Cos[b1] Sin[a1] Sin[b2]))^2)-2 Cos[a1] Cos[a2] Cos[b1] Cos[b2] (Cos[b1] Cos[b2] Sin[a1] Sin[a2]+Sin[b1] Sin[b2])+((Cos[a1])^2) ((Cos[b1])^2) (((Cos[b2])^2) ((Sin[a2])^2)+((Sin[b2])^2))))]

randomcylaxis:= {RandomReal[],RandomReal[{0,2 π}],RandomReal[{-π,π}],RandomReal[{0,2 π}],RandomReal[{-π,π}]}

makesphere[rad_, cylstruct_] := Block[{vec, time}, vec= cylaxis[time,cylstruct]; 
time/.Solve[vec.vec==rad^2, time]]

cylist = {{{1,0,0,0,0},0.25}};
Length[cylist]

plotcyls[clist_] := Module[{t1,t2,curax,curad,maxrad},
maxrad=Max[Table[clist[[i,2]],{i,1,Length[clist]}]];
Graphics3D[Table[{curax,curad}=clist[[i]];{t1,t2}=makesphere[10,curax];{ColorData["RoseColors"][0.66curad/(maxrad)],
Specularity[White,50],EdgeForm[],Cylinder[{cylaxis[t1,curax],cylaxis[t2,curax]},curad]},{i,1,Length[clist]}],Boxed->False,Background->Lighter[ Lighter[Gray]],Lighting->{{"Ambient",Darker[Gray]},{"Directional",White,{10,0,0},{-10,0,0}}}]]

cylist = {randomcyl}
Length[cylist]

minpos[a_List]:= Module[{min=Min[a]},{min,Flatten[Position[a,min]][[1]]}]

icount=0;

While[Length[cylist] ≤200 && icount++ < 1000000,test = {RandomReal[{0,5}],RandomReal[{0,2 π}],RandomReal[{-π,π}],RandomReal[{0,2 π}],RandomReal[{-π,π}]};distlist = Table[axialdist[test,cylist[[i,1]]]-cylist[[i,2]],{i,1,Length[cylist]}];
{mindist,nbrcyl}= minpos[distlist];
If[mindist> 0, AppendTo[cylist,{test,Min[{mindist,RandomReal[{0.9,1.1}]}]}]]]

cylist

plotcyls[cylist]