Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

For simplicity's sake, let the polymer's composing particle be integer 3-tuples (on 3D Cartesian grid), such that:

  1. Polymer-ness -- for all particle there is an directly adjacent neighbor. i.e in directions of {{0,0,1}, {0,0,-1}, {0,1,0}, {0,-1,0}, {1,0,0}, {-1,0,0}}

  2. Of loose spacial density -- there's no more than $d_n$ particles in a radius $n$ sphere around each particle

The plot looks interesting. DLA-esque. In fact, if we let $ d_n = <3^n> $, it appears rather dimensionless i.e. factal looking

However my code was slow and it's a pain to generate a set larger than 400 particles (~7mins)

200 particles

Moreover, how can I spin the plot independently around 2 spatial axis? (like that asteroid in Armageddon IYNWIM)

I tried ViewAngle and use 2 RotationTransform to rotate the vector. However, that's not how it seems to works.

Clear["Global`*"];
colorPalette = (#/3 + 0.2) & /@ {1, 2, 3};
radiusPlaette = {0.63, 0.53, 0.46};
adjLimit = {3, 9, 27, 81};
idx = Range[Length@adjLimit];
dir = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 
 0, -1}};
pos = {{0, 0, 0}}; test = {};

adj[x_, d_] := Cases[pos, _?(EuclideanDistance[x, #] <= d &)];

Dynamic[Length[pos]]

(*Polymer, nice-looking spheres*)
Dynamic[
Graphics3D[{Hue[Take[colorPalette, {Mod[#2, 3] + 1}], 1, 0.1],
  Sphere[#1 - {0.5, 0.5, 0.5} - Mean[pos], 
   Take[radiusPlaette, {Mod[#2, 3] + 1}]]} & @@@ 
   MapThread[List, {pos, Range[Length@pos]}], Axes -> False, 
  Boxed -> False, Lighting -> Automatic, SphericalRegion -> True ]]

(*Alternative layout
(*Colored Boxes*) 
Dynamic[Graphics3D[{Hue[1.*#2/3+0.15,0.7,1],Cuboid[#1-{0.5,0.5,0.5}-
Mean[pos]]}&@@@MapThread[List,{pos,Range[Length@pos]}],Axes->
False,Boxed->False,Lighting->{{"Ambient",White}}, SphericalRegion->True]]
*)

Timing[
 Do[
  out = Catch[While[True,
     stem = RandomChoice[pos];  
     rnd = stem + RandomChoice[dir];
     listadj1 = Length /@ (adj[rnd, #] & /@ idx);
         adj1 = And @@ MapThread[(#1 <= #2) &, {listadj1, adjLimit}];
     listadj2 = (Function[x, Max[Length /@ (adj[#, x] & /@ adj[rnd, x])]] /@ idx);
         adj2 = And @@ MapThread[(#1 <= #2) &, {listadj2, adjLimit}];

     If[! MemberQ[pos, rnd] && adj1 && adj2,(*Debug*)
      AppendTo[test, listadj2]; Throw[rnd]];
     ]];
  AppendTo[pos, out];
  , {399}]
 ]
(*Debug*)
Last@Sort[test]
Tally[Length /@ (adj[#, 2] & /@ pos)]

.

share|improve this question

1 Answer 1

I don't understand your original code very well (particularly the definitions of listadj1 and listadj2), but here's my naive translation for your description, I only rewrote the simulation part:

densitytest = 
  With[{dis2 = idx^2, limit = adjLimit}, 
   Compile[{{pos, _Real, 2}}, 
    Times @@ 
     UnitStep@(limit + 1 - 
        Max /@ Transpose[Total /@ Map[UnitStep[dis2 - #] &,
            Map[Total, Outer[Plus, pos, -pos, 1]^2, {-2}], {-1}]])]];

AbsoluteTiming[Do[
  Module[{pos2}, While[pos2 =!= pos,
    If[! MemberQ[pos, rnd = RandomChoice[pos] + RandomChoice[dir]], 
     pos2 = Flatten[{pos, {rnd}}, 1]; 
     If[1 == densitytest@pos2, pos = pos2]]]], {399}]]

enter image description here

I believe the code can be further optimized, but now I'd like to stop here and go to bed :)


Update:

I managed to further optimize the code, only about 40 seconds are needed for 400 particles now:

densitytest = 
  With[{dis2 = idx^2, limit = adjLimit}, 
   Compile[{{pos, _Real, 2}}, 
    Times @@ UnitStep@(limit + 1 - Max /@ Total /@ UnitStep@Outer[Plus, dis2,
            -Total@Transpose[Outer[Plus, pos, -pos, 1]^2, {2, 3, 1}]])]];

AbsoluteTiming[Do[
  Module[{pos2, toselect = Complement[Flatten[Outer[Plus, pos, dir, 1], 1], pos]},
    While[pos2 =!= pos,
    pos2 = Flatten[{pos, {RandomChoice@toselect}}, 1]; 
    If[1 == densitytest@pos2, pos = pos2]]], {399}]]
share|improve this answer
    
96 secs! That's virtuoso to me \^O^/ –  秦紀維 May 20 at 18:23
    
@秦紀維 And it's not the end, see my edit. –  xzczd May 21 at 4:29

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.