(*BB[*)(* Initialization *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
ClearAll[ball, pos, mom, bins, peg];

ball = <|
  r -> 0.1,         (* Radius of the ball *)
  n -> 50          (* Number of simultaneous balls *)
|>;

peg = <|
  r -> 0.25         (* Radius of the pegs *)
|>;

sim = <|
  \[CapitalDelta]t -> 0.05,    (* Higher Value need to correct collision positions *)
  \[Epsilon] -> 0.4,           (* 1.0 means all energy back in the ball *)
  g -> {0, -9.81},             (* You can tilt the board *)
  n -> 300                    (* Total Number of Balls *)
|>;

(* Initialize Position and Momentum Matrices *)
randomBall := {RandomReal[{-0.05, 0.05}], RandomReal[{0, 1}] + 10};
ball[pos] = Table[randomBall, ball[n]];
ball[mom] = Table[0, ball[n], 2];
ball[bins] = {};

(* positions of bin rectangles : only for visualizaiton *)
rectangleBins = Table[0.01, {i,-8,8,1}];

(* Peg position *)
peg[pos] = Flatten[Table[{ii + 0.5 jj, jj}, {ii, -15, 10}, {jj, 0, 8}], 1];

visPositions = ball[pos];

simulateNext[] := With[{},
  (* Step *)
  (* Scheme: Backward Euler *)

  (* Next Step Momentum *)
  ball[mom] += Threaded[sim[\[CapitalDelta]t] sim[g]];

  (* Next Position Prediction *)
  next = ball[pos] + sim[\[CapitalDelta]t] ball[mom];

  (* Distance: Do not parallelize, it's already in the algo. *)
  dist = DistanceMatrix[next, peg[pos]];

  (* Collisions between current and future step *)
  col = MapIndexed[
    {v, i} |-> {
      SelectFirst[v, (# < ball[r] + peg[r]) & -> "Index"], i[[1]]
    } /. ({Missing[_], _} -> Nothing),
    dist
  ];

  (* Correct Momentum to include collisions *)
  If[Length[col] > 1,
    {
      up, val
    } = Transpose[
       MapApply[
        {p, b} |-> {
          b,
          sim[\[Epsilon]] * Norm[ball[mom][[b]]] Normalize[next[[b]] - peg[pos][[p]]]
        },
        col
      ]
    ];
    momt = ball[mom];
    momt[[up]] = val;
    ball[mom] = momt;
  ];

  (* Update pos with corrected momentum *)
  ball[pos] += sim[\[CapitalDelta]t] ball[mom];

  (* Filter through the balls, remove/reset and register those exiting *)
  exited = Select[ball[pos], #[[2]] < 0 & -> "Index"];
  ingame = Select[ball[pos], #[[2]] > 0 & -> "Index"];
  
  If[Length[exited] > 0,
    ball[bins] = Join[ball[bins], ball[pos][[exited, 1]]]
  ];

  (* Remove them from the pos and mom matrices *) 
  ball[pos] = ball[pos][[ingame]];
  ball[mom] = ball[mom][[ingame]];

  (* Insert new ball *)
  left = sim[n] - Length[ball[bins]] - Length[ball[pos]];
  toAdd = Min[left, ball[n] - Length[ball[pos]]];
  Do[
    AppendTo[ball[pos], randomBall];
    AppendTo[ball[mom], {0, 0}],
    toAdd
  ];

  (* Update visible positions *)
  visPositions = ball[pos];

  (* Update bins *)
  If[Length[ball[bins]] > 5,  With[{h = HistogramList[ball[bins], {-8, 8, 1}]},
    rectangleBins = Rescale @ ReplacePart[
        ConstantArray[0.01, Length@h[[1]]],           (* 17 entries at .01 *)
        Thread[Most@h[[1]] + 9 -> h[[2]]]              (* now integers 1,2,…16 *)
      ]
  ]];
];