0

I am trying to create a Monte Carlo simulation for a physical situation. The situation is a layer of disks are lying on a plane. The disks have a known mean radius and standard deviation, and a known mean center-to-center distance (measured compared to it's nearest $n$ neighbors) and standard deviation. Both follow a normal distribution.

I would like to create a simulation that, given the mean values and standard deviations produces a list of disk radii and center locations that match the statistics. The disks should not overlap and fit within a given boundary.

David G. Stork
  • 41,180
  • 3
  • 34
  • 96
bicarlsen
  • 295
  • 2
  • 8
  • What have you tried? – David G. Stork Oct 04 '17 at 16:52
  • I haven't tried anything yet. My main problem is I don't know how to check if the disks are overlapping. – bicarlsen Oct 04 '17 at 16:53
  • 1
    Take the centers, two at a time, and check that the sum of the associated two radii is less than the separation between the centers. Note: You cannot have true normal distribution under your constraints. – David G. Stork Oct 04 '17 at 16:58
  • So simple, thanks. However, if two particles do overlap how can I reassign either the radii or the center position so the distributions are not disturbed too much? – bicarlsen Oct 04 '17 at 17:00
  • I would just randomly delete one of the overlapping disks. Again, any such change guarantees your distribution with not be Gaussian. – David G. Stork Oct 04 '17 at 17:06
  • 1
    start here. https://mathematica.stackexchange.com/questions/126486/generate-randomly-sized-non-overlapping-disks. I guess this is not quite a dup due to the matching statistics part, but not really on topic with that either. – george2079 Oct 04 '17 at 17:14
  • 2
    your target packing density is a big factor in this problem. If you are going after a fairly dense packing the problem becomes nearly intractable. At low density you can readily pack with your target size and mean spacing statistics. Capturing the spacing distribution can be done but it is a challenge. – george2079 Oct 04 '17 at 17:34
  • Thanks @george2079, indeed the packing density is quite high (61% area coverage). It seems the best option is to add particles by hand (luckily this is do-able) and adjust parameters until they somewhat match all relevant values. – bicarlsen Oct 04 '17 at 18:20
  • I would be happy to post my solution. What is the etiquette for this? Should I post it as an answer? – bicarlsen Oct 04 '17 at 18:20
  • Yes. Just post your result as an answer. – bbgodfrey Oct 04 '17 at 23:32
  • You actually should post it, otherwise the question sits on the unanswered list. Plus maybe folks who see it will come up with better approaches. – george2079 Oct 05 '17 at 16:36
  • 2
    Could you provide some explicit data to play with? Also you could try incorporating RegionDisjoint. See the 'Neat Example' in the documentation here. – Greg Hurst Oct 06 '17 at 17:42
  • @ChipHurst, I didn't know about that function. Looks like it could definitely make things cleaner for checking disk overlaps. I like using RegionUnion for checking with the boundary intersection because it also works if a disk is completely outside the region.

    I provided a minimum example of data as a comment in the answer.

    – bicarlsen Oct 06 '17 at 20:24
  • @bicarlsen , Next time post the work of your own effort and have this opportunity to debug and learn what is that is keeping you from completing the problem. – Jose Enrique Calderon Oct 07 '17 at 12:55

1 Answers1

1

To solve the problem I created some function to help visualize the placement and statistics of the disks. For each functions disks is a list with elements {radius, {xCenter, yCenter}} and boundary is a vector {width, height} describing the bounding area, which is centered at the origin.

areaCoverage calculates the percent of the enclosed area covered by particles.

raiusStats calculates relevant statistics of the disk radii.

visualizeLayer creates a graphic with the disks and bounding area, and gives a list of any disks that intersect each other or the boundary.

areaCoverage[disks_, boundary_] := Module[
   {diskArea, boundaryArea},
   diskArea = Pi*Total[disks[[All, 1]]^2];
   boundaryArea = Times @@ boundary;

   N[diskArea/boundaryArea*100]
];

radiusStats[disks_] := Module[
   {radii},
   radii = disks[[All, 1]] ;

   N[ {Mean[radii], StandardDeviation[radii] } ]
];

visualizeLayer[disks_, boundary_] := Module[
   {d, b, pairs, intersecting, bInt},
   b = Rectangle[ -boundary/2, boundary/2 ];
   d = Table[Disk[disk[[2]], disk[[1]] ], {disk, disks}];
   pairs = Subsets[d, {2} ];
   pairs = Table[With[{
      intersect = (RegionMeasure@RegionIntersection@pair != 0)
     },
     {pair, intersect}
    ],
   {pair, pairs}
   ];
  intersecting = Select[pairs, #[[2]] == True &][[All, 1]];

  bInt = Table[With[{
      intersect = ! (RegionUnion@{b, disk} === b)
      },
     {disk, intersect}
     ],
    {disk, d}
    ];
  bInt = Select[bInt, #[[2]] == True &][[All, 1]];

 {Show[Graphics[ {FaceForm[], EdgeForm[Black], b} ], d // Graphics, 
   ImageSize -> 500], intersecting, bInt}
];

I started with a list of three disks, adding new ones one-by-one playing with the parameters and placement to match the statistics I wanted.

While not automated, these functions certainly helped create scenarios more quickly than would otherwise be possible. The visualization of the particles made adding new ones quite easy, and most of the time was spent making small adjustments to match the statistics to given values.

Here is a minimum set of data:

boundary = {200, 200};

disks={ 
    { 10, {0, 0} },
    { 10, { 50, 50 } },
    { 10, {-50, -50} },
    { 20, {-40, -50} },
    { 30, {-80, 0} } 
};
Greg Hurst
  • 35,921
  • 1
  • 90
  • 136
bicarlsen
  • 295
  • 2
  • 8
  • 2
    I think what people would need to work this problem is an example of the statistics you are after. – george2079 Oct 06 '17 at 20:49
  • @george2079 that's a good point. I was mentioned in the question that I was interested in the distribution of the particle radii, and the center-to-center distance between the n nearest neighbors. The nearest neighbors calculation ended up being overkill, so I switched to coverage percentage instead. – bicarlsen Oct 07 '17 at 20:37