clc, clear, close all

% Update this link with the output from §4 of Multistorey_examples.m
sims_data = dlmread('para_cases_output_new_20191220142149.txt','\t',2,0);
n_cases = 1000;
Pcr_FEA = zeros(1,n_cases);
Pcr_ME = zeros(1,n_cases);
for sim_num = 1:n_cases
    
    % Read the ABAQUS output
    case_ID = fopen(strcat('inp_para_case_',num2str(sim_num),'.dat'),'r');
    case_data = fscanf(case_ID,'%c');
    fclose('all');
    case_data = case_data(strfind(case_data,'MODE NO'):strfind(case_data,'THE ANALYSIS HAS BEEN COMPLETED')-1);
    case_data = case_data(strfind(case_data,'1')+1:length(case_data));
    case_data = case_data(find(~isspace(case_data)));
    Pcr_FEA(sim_num) = str2double(case_data);
    
    % Read the MATLAB output
    Pcr_ME(sim_num) = sims_data(sim_num,22);
end

%Diff = (Pcr_ME-Pcr_FEA)./max(Pcr_ME,Pcr_FEA);
Diff = (Pcr_ME-Pcr_FEA)./Pcr_FEA;
Diff(Pcr_ME<0) = NaN;
avg = mean(Diff(~isnan(Diff)));
stdev = std(Diff(~isnan(Diff)));
minim = min(Diff(~isnan(Diff)));
maxim = max(Diff(~isnan(Diff)));

figure('Renderer', 'painters', 'Position', [10 10 550 230])
Diff = Diff*100;
histogram(Diff,'FaceColor','k')
set(gca, 'FontName', 'Times New Roman');
xlabel('Diff (%)')
ylabel('Number of realizations')
grid on