%% Pushing the limit of all diamond anvils
% Edoweiss
% 2018-06-15

close all; clear all; clc

%% simulations parameters
trials = 100000;
available_potential = 13;
initial_level = 5;

%% Setup
pTable = [  1
            1
            1
            1
            1
            0.88
            0.78
            0.68
            0.59
            0.51
            0.50
            ];

d_levels = [];

%% Diamond Anvils
% we keep trying until we run out of diamond anvils or potential
for i = 1:trials;
    pot = available_potential; % potential remaining
    lvl = initial_level; % current enhancement level
    while pot > 0 && lvl < 40
        if lvl < length(pTable)
            check = pTable(lvl+1);
        else
            check = 0.50;
        end
        if rand < check % success
            lvl = lvl + 1;
        else %failure
            pot = pot - 1;
        end
    end
    d_levels = [d_levels; lvl];
end

%%
fig1 = figure;
[h_d,e_d]  = histcounts(d_levels,'Normalization','pdf');
for i = 1:length(e_d)-1 % get centers of bins
    c_d(i) = (e_d(i)+e_d(i+1))/2;
end

area(c_d,h_d)

ylabel('Probability Density')
xlabel('enhancement level')
saveas(fig1,'dia-maxpdf.png')

%%
fig2 = figure;
[h_d,e_d]  = histcounts(d_levels,'Normalization','cdf');
for i = 1:length(e_d)-1 % get centers of bins
    c_d(i) = (e_d(i)+e_d(i+1))/2;
end
c_d = [5 c_d];
h_d = [1 1-h_d];
area(c_d,h_d)

ylabel('Probability')
xlabel('enhancement level')
saveas(fig2,'dia-maxcdf.png')