% Simulates irreversible edge tile binding
% B     B
%  B   B
%   B B
%    C
corner_tiles = 500;
simulations = 10;
corner_assembly_sizes = zeros(corner_tiles*2,simulations);

for i=1:simulations
  edge_tiles = 1000;
  corner_tiles = 500;

  corner_tile_left_sizes = zeros(1,corner_tiles);
  number_of_left_additions = zeros(1,corner_tiles);
  corner_tile_right_sizes = zeros(1,corner_tiles);
  number_of_right_additions = zeros(1,corner_tiles);

  i
  edge_assemblies = ones(1,edge_tiles);
  number_of_edge_assemblies = edge_tiles;
%  edge_bond = zeros(1,edge_tiles);

  percentage_completion = 1;
  circular_assemblies = 0;

  while number_of_edge_assemblies > ...
	((1-percentage_completion) * edge_tiles)
    % Find a random possible bond
    corner_bond_ratio = corner_tiles / ...
	(2 * corner_tiles + number_of_edge_assemblies * 2);
    corner_bond_ratios(number_of_edge_assemblies) = ...
	corner_bond_ratio * 2;

    bond = random ('Uniform', 0, 1);
    
    % And a random possible other bond:
    lucky_other_assembly = random('Discrete Uniform', ...
				  number_of_edge_assemblies);
    
    if bond < corner_bond_ratio
      % Bond is from a left corner tile: pick the assembly
      lucky_assembly = ceil(bond / corner_bond_ratio * corner_tiles);
      
      corner_tile_left_sizes(lucky_assembly) = ...
	  corner_tile_left_sizes(lucky_assembly) + ...
	  edge_assemblies(lucky_other_assembly);
      number_of_left_additions(lucky_assembly) = ...
	  number_of_left_additions(lucky_assembly)+1;
      edge_assemblies(lucky_other_assembly) = ...
	  edge_assemblies(number_of_edge_assemblies);
      number_of_edge_assemblies = number_of_edge_assemblies - 1;
%      edge_bond(number_of_edge_assemblies) = 1;
    elseif bond < 2 * corner_bond_ratio
      % Bond is from a right corner tile
      lucky_assembly = ceil((bond - corner_bond_ratio)/corner_bond_ratio ...
			    * corner_tiles);
      
      corner_tile_right_sizes(lucky_assembly) = ...
	  corner_tile_right_sizes(lucky_assembly) + ...
	  edge_assemblies(lucky_other_assembly);
      number_of_right_additions(lucky_assembly) = ...
	  number_of_right_additions(lucky_assembly)+1;
      edge_assemblies(lucky_other_assembly) = ...
	  edge_assemblies(number_of_edge_assemblies);
      number_of_edge_assemblies = number_of_edge_assemblies - 1;
%      edge_bond(number_of_edge_assemblies) = 1;
    else
      % Bond is between two edge assemblies
      lucky_assembly = ceil((bond - 2 * corner_bond_ratio) / ...
			    (1 - 2 * corner_bond_ratio) * ...
			    number_of_edge_assemblies);
      if lucky_assembly == lucky_other_assembly 
	'Circular assembly; skipping';
	circular_assemblies = circular_assemblies + 1;
      else
	edge_assemblies(lucky_assembly) = edge_assemblies(lucky_assembly) ...
	    + edge_assemblies(lucky_other_assembly);
	edge_assemblies(lucky_other_assembly) = ...
	    edge_assemblies(number_of_edge_assemblies);
	number_of_edge_assemblies = number_of_edge_assemblies - 1;
      end;
    end;
  end;
  corner_assembly_sizes(:,i) = [corner_tile_left_sizes,corner_tile_right_sizes]';
end;
  
% Display the current assemblies
  
%max_corner_size = max(max(corner_assembly_sizes));
%min_corner_size = min(min(corner_assembly_sizes));
max_corner_size = 8;
min_corner_size = 0;
bins = 9;

for i=1:simulations
  corner_assembly_sizes_hist(:,i) = histc(corner_assembly_sizes(:, ...
						  i),linspace(min_corner_size, max_corner_size,bins));
end;

corner_assembly_sizes_mean = mean(corner_assembly_sizes_hist,2);
corner_assembly_sizes_std = std(corner_assembly_sizes_hist,0,2);

figure(3);

hold off;
bar(linspace(min_corner_size,max_corner_size,bins),corner_assembly_sizes_mean,'y');
hold on;
axis([-1 9 0 1000]);
errorbar(linspace(min_corner_size,max_corner_size,bins), ...
	 corner_assembly_sizes_mean,corner_assembly_sizes_std,'.');
hold off;
xlabel('Corner Assembly Length, in Number of Tiles','FontSize',20);
ylabel('Number of Assemblies','FontSize',20);
title('Completed Irreversible Aggregation 1:1', ...
      'FontSize',20);
set(gca,'FontSize',20);
saveas(gcf, 'irreversible_aggregation_complete_length_distribution_one_to_one', ...
       'fig');
print -depsc2 irreversible_aggregation_complete_length_distribution_one_to_one.eps