Poisson Branching Process (Critical)

M = 10;
T = zeros(M, 1);
S = zeros(M, 1);

for i=1:M
  current_generation = 1;
  t = 0;
  size = 1;
  while current_generation > 0
    %% current_generation = sum(2 * round(rand(current_generation, 1)));
    current_generation = sum(poissrnd(1, [1, M]));
    size = size + current_generation;
    t = t + 1;
  end
  T(i) = t;
  S(i) = size;
end

fig = figure('visible', 'off');
scatter(T, S);
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');
grid on;
print(fig, "critical.png")