format long g;

% define data processing constants
SET_SIZE = 16;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%load data file
%data = load('chsh_0_20141124.dat');
%data = load('chsh_180_20141125');
data_1 = load('chsh_h_optimise_20141129.dat');
data_2 = load('chsh_h_optimise_20141201.dat');
data_3 = load('chsh_h_optimise_20141202.dat');
data = [data_1; data_2; data_3];

% evaluate the number of complete set
NO_OF_SET = floor(size(data, 1)/SET_SIZE);

% truncating incomplete set
data(NO_OF_SET*SET_SIZE + 1 :end, :) = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data_singles = data(:, [2, 3]);
data_2_fold = reshape(data(:, 4), 4, 4*NO_OF_SET);

data_singles = data_singles./[data(:, 1), data(:, 1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a0_1 = data_singles(1:16:end, 1);
a0_2 = data_singles(3:16:end, 1);
a0_3 = data_singles(5:16:end, 1);
a0_4 = data_singles(7:16:end, 1);
a0 = [a0_1; a0_2; a0_3; a0_4];

a0_90_1 = data_singles(2:16:end, 1);
a0_90_2 = data_singles(4:16:end, 1);
a0_90_3 = data_singles(6:16:end, 1);
a0_90_4 = data_singles(8:16:end, 1);
a0_90 = [a0_90_1; a0_90_2; a0_90_3; a0_90_4];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a1_1 = data_singles(9:16:end, 1);
a1_2 = data_singles(11:16:end, 1);
a1_3 = data_singles(13:16:end, 1);
a1_4 = data_singles(15:16:end, 1);
a1 = [a1_1; a1_2; a1_3; a1_4];

a1_90_1 = data_singles(10:16:end, 1);
a1_90_2 = data_singles(12:16:end, 1);
a1_90_3 = data_singles(14:16:end, 1);
a1_90_4 = data_singles(16:16:end, 1);
a1_90 = [a1_90_1; a1_90_2; a1_90_3; a1_90_4];

plot(a1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b0_1 = data_singles(1:16:end, 2);
b0_2 = data_singles(2:16:end, 2);
b0_3 = data_singles(9:16:end, 2);
b0_4 = data_singles(10:16:end, 2);
b0 = [b0_1; b0_2; b0_3; b0_4];

b0_90_1 = data_singles(3:16:end, 2);
b0_90_2 = data_singles(4:16:end, 2);
b0_90_3 = data_singles(11:16:end, 2);
b0_90_4 = data_singles(12:16:end, 2);
b0_90 = [b0_90_1; b0_90_2; b0_90_3; b0_90_4];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b1_1 = data_singles(5:16:end, 2);
b1_2 = data_singles(6:16:end, 2);
b1_3 = data_singles(13:16:end, 2);
b1_4 = data_singles(14:16:end, 2);
b1 = [b1_1; b1_2; b1_3; b1_4];

b1_90_1 = data_singles(7:16:end, 2);
b1_90_2 = data_singles(8:16:end, 2);
b1_90_3 = data_singles(15:16:end, 2);
b1_90_4 = data_singles(16:16:end, 2);
b1_90 = [b1_90_1; b1_90_2; b1_90_3; b1_90_4];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
index = 1:length(a0);
plot(index, a0, 'b', index, a0_90, 'g', index, a1, 'r', index, a1_90, 'c',...
    index, b0, 'm', index, b0_90, 'y', index, b1, 'k', index, b1_90, 'b-.'),...
    xlabel('No. of Measurments'),...
    ylabel('Normalized Counts'),...
    title('Normalized Counts versus No. of Measurments'),...
    grid;

legend('a0', 'a0_9_0', 'a1', 'a1_9_0', 'b0', 'b0_9_0', 'b1', 'b1_9_0');

mean(a0)
std(a0)

mean(a0_90)
std(a0_90)

mean(b0)
std(b0)

mean(b0_90)
std(b0_90)


mean(a1)
std(a1)

mean(a1_90)
std(a1_90)

mean(b1)
std(b1)

mean(b1_90)
std(b1_90)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

std(a0)/mean(a0)
std(a0_90)/mean(a0_90)

std(b0)/mean(b0)
std(b0_90)/mean(b0_90)

std(a1)/mean(a1)
std(a1_90)/mean(a1_90)

std(b1)/mean(b1)
std(b1_90)/mean(b1_90)

test = [mean(a0), mean(a0_90), mean(a1), mean(a1_90)];
test_a = test/max(test) % max is a0

test = [mean(b0), mean(b0_90), mean(b1), mean(b1_90)];
test_b = test/max(test) % max is b0_90

delta_a0 = std(a0);
delta_a0_90 = std(a0_90);
delta_a1 = std(a1);
delta_a1_90 = std(a1_90);
delta_a = [std(a0), std(a0_90), std(a1), std(a1_90)];

delta_b0 = std(b0);
delta_b0_90 = std(b0_90);
delta_b1 = std(b1);
delta_b1_90 = std(b1_90);
delta_b = [std(b0), std(b0_90), std(b1), std(b1_90)];

a = [mean(a0), mean(a0_90), mean(a1), mean(a1_90)];
a_norm = a/max(a);
delta_a_norm = sqrt((1/mean(a0))^2*delta_a.^2 + (a/mean(a0)^2).^2*delta_a0^2);

b = [mean(b0), mean(b0_90), mean(b1), mean(b1_90)];
b_norm = b/max(b);
delta_b_norm = sqrt((1/mean(b0_90))^2*delta_b.^2 + (b/mean(b0_90)^2).^2*delta_b0_90^2);

term_1_ele = powerset(b_norm, delta_a_norm);
term_1 = term_1_ele(:, 1).*term_1_ele(:, 2);

term_2_ele = powerset(a_norm, delta_b_norm);
term_2 = term_2_ele(:, 1).*term_2_ele(:, 2);

ab_norm_ele = powerset(a_norm, b_norm);
ab_norm = ab_norm_ele(:, 1).*ab_norm_ele(:, 2);
delta_ab_norm = sqrt(term_1.^2 + term_2.^2);

correction_factor = 1./ab_norm;
delta_correction_factor = delta_ab_norm./ab_norm.^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


test_ab = powerset(test_a, test_b)
test_ab(:, 1).*test_ab(:, 2)
powerset(1:4, 1:4)
plot(test_ab(:, 1).*test_ab(:, 2), 'x-')

correction_factor = 1./(test_ab(:, 1).*test_ab(:, 2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
