% Representation of the experimental dataset for 
% "Short Communication: A simple and accurate method of measuring the
% Zeta-Potential of microfluidic channels"

close all;clear;clc

% General parameters
eta = 1e-3; % Pa s
L = 1e-2;  % m
eps = 80*8.85e-12; % C/V/m
E = 40/L; % V/m

% Definitions of color for the plots
colorR = [225,75,75]/255;
colorB = [56 117 176]/255;

%% Experimental obtention of the pre-factor

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 1.5 mS/m
% Applied voltage
V0 = [0.1 .2 .3 .4 .5 .75 1 1.25 1.5 1.75 2]*200; % V
% Measurements with sigma 1
V_1 = [8.304 21.31 29.46 38.72 47.46 70.47 94.62 120.1 143.4 167.5...
    190.5]*10^(-3); % V
% Measurement with sigma 2
V_2 = [7.321 17.30 25.74 33.30 41.29 61.32 81.87 102.8 123.8 145.8...
    166.1]*10^(-3); % V

figure('Name','Voltage Divider. 1.5 mS/m KCl Electrolyte')
box on
hold on
ax=gca;
ax.FontSize = 20;
ax.LineWidth = 1.5;
ax.FontName = 'Times New Roman';
plot(V0,V_1*1e3,'or','color',colorR,'MarkerSize',10,'LineWidth',2)
plot(V0,V_2*1e3,'ob','color',colorB,'MarkerSize',10,'LineWidth',2)
legend('\sigma_1','\sigma_2','Location','best')
xlabel('V_0 (Volts)'); ylabel('V (mV)')
title('Voltage Divider. 1.5 mS/m KCl Electrolyte')
p_1 = polyfit(V0,V_1,1);
p_2 = polyfit(V0,V_2,1);

ratio_1 = (1-p_1(1))/p_1(1); % R_ch(sigma_1)/R
ratio_2 = (1-p_2(1))/p_2(1); % R_ch(sigma_2)/R
twenty = ratio_2/(ratio_2-ratio_1);
factor_1mSm = twenty*ratio_1;

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 5.2 mS/m
% In this case it was assumed a 5% in conductivity
% Applied voltage
V0 = [0.05 0.1 .2 .3 .4 .5 .75 1 1.25 1.5 2]*200; % V
% Measurement with sigma 1
V = [11.55 22.7 46.0 69 92 114.4 172 230.5 285 327 441]*10^(-3); % V

figure('Name','Voltage Divider. 5.2 mS/m KCl Electrolyte')
box on
hold on
ax=gca;
ax.FontSize = 20;
ax.LineWidth = 1.5;
ax.FontName = 'Times New Roman';
plot(V0,V*1e3,'ok','MarkerSize',10,'LineWidth',2)
legend('\sigma_1','Location','best')
xlabel('V_0 (Volts)'); ylabel('V (mV)')
title('Voltage Divider. 5.2 mS/m KCl Electrolyte')
p = polyfit(V0,V,1);

ratio_5mSm = (1-p(1))/p(1); % R_ch/R

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 11.4 mS/m
% Applied voltage
V0 = [0.05 0.1 .2 .3 .4 .5 .75 1 1.25 1.5 2]*200; % V
% Measurements with sigma 1
V_1 = [24.7 50.4 102.2 152.8 204.5 256.2 382.8 512.1 641.5 784.1 1051]...
    *10^(-3); % V
% Measurement with sigma 2
V_2 = [24.05 48.75 98.00 145.5 194.2 244.2 364.5 489.9 611.8 746.5 996]...
    *10^(-3); % V

figure('Name','Voltage Divider. 11.4 mS/m KCl Electrolyte')
box on
hold on
ax=gca;
ax.FontSize = 20;
ax.LineWidth = 1.5;
ax.FontName = 'Times New Roman';
plot(V0,V_1*1e3,'or','color',colorR,'MarkerSize',10,'LineWidth',2)
plot(V0,V_2*1e3,'ob','color',colorB,'MarkerSize',10,'LineWidth',2)
legend('\sigma_1','\sigma_2','Location','best')
xlabel('V_0 (Volts)'); ylabel('V (mV)')
title('Voltage Divider. 11.4 mS/m KCl Electrolyte')
p_1 = polyfit(V0,V_1,1);
p_2 = polyfit(V0,V_2,1);

ratio_1 = (1-p_1(1))/p_1(1); % R_ch(sigma_1)/R
ratio_2 = (1-p_2(1))/p_2(1); % R_ch(sigma_2)/R
twenty = ratio_2/(ratio_2-ratio_1);
factor_11mSm = twenty*ratio_1;

%% Experimental zeta potential estimation data for non-treated surfaces

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 1.5 mS/m
V = [24.40 24.02 23.47 22.82 22.07 21.56 21.07 20.88 20.71 20.70 20.52...
    20.55 20.42 20.37 20.41 20.30 20.20 20.18 20.16 20.12 20.00 20.00...
    20.01 19.86 19.93 19.78 19.74 19.58 19.61 19.69 19.54 19.55 19.46...
    19.35 19.24 19.13 19.06 18.95 18.97 18.85 18.72 18.70 18.59 18.68...
    18.45 18.57 18.42 18.38 18.36 18.31 18.30 18.41 18.20 18.18 18.18...
    18.17 18.15 18.20 18.20 18.17 18.03 17.99 17.94 18.00 17.96 17.91...
    17.92 17.92 17.90 17.88 17.87 17.81 17.77 17.74 17.73 17.71 17.71...
    17.71 17.59 17.75 17.68 17.54 17.63]*1e-3; % V

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_1mSm_NonTreated = t;
V_1mSm_NonTreated = V;

figure('Name', 'Voltage Monitoring. 1.5 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t(5:59),V(5:59)*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',...
    colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 7;
t02 = 40;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 1.5 mS/m KCl Non-Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
factor_1mSm = 1.733608523735045e+04; 
zetaP_1mSm_NonTreated = eta/(eps*E^2)*factor_1mSm*sp;
Error_zetaP_1mSm_NonTreated = eta/(eps*E^2)*factor_1mSm*Esp;

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 5.2 mS/m
V = [4169 4153 4111 4096 4056 4071 4051 4024 4046 4035 4038 4037 4031 ... 
    4005 4005 4014 4013 3981 3982 3991 3995 3990 3969 3980 3950 3957 ...
    3956 3961 3958 3939 3949 3921 3939 3913 3912 3906 3902 3903 3887 ...
    3897 3876 3855 3880 3851 3863 3860 3867 3854 3848 3845 3833 3858 ...
    3837 3837 3836 3833 3846 3832 3839 3834 3833 3828 3817 3836 3801 ...
    3826 3834 3820 3808 3821 3795 3805 3807 3822 3803 3818 3813 3808 ...
    3777 3797 3799 3803 3793 3782 3780 3784 3779 3775 3796 3771 3765 ...
    3789 3766 3768 3783 3771 3775 3763 3769 3759 3774 3754 3755 3753 ...
    3762 3747 3745 3748 3739 3758 3753 3740 3764 3758 3745 3746 3743 ...
    3738 3730 3729 3750 3747 3733]'*1e-5;

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_5mSm_NonTreated = t;
V_5mSm_NonTreated = V;

figure('Name', 'Voltage Monitoring. 5.2 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t(1:60),V(1:60)*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',...
    colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 6;
t02 = 45;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 5.2 mS/m KCl Non-Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
ratio = 906.4916;
zetaP_5mSm_NonTreated = 20*eta*sp/(eps*E^2)*ratio;
Error_zetaP_5mSm_NonTreated = 20*eta*Esp/(eps*E^2)*ratio;

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 11.4 mS/m
V = [10.5 10.2 10.4 9.9 9.8 9.5 9.4 9.2 9.1 9.4 8.9 8.8 8.9 8.8 8.6 8.8...
    8.5 8.5 8.5 8.1 8.2 8.1 8.0 7.7 7.6 7.8 7.4 7.4 7.0 7.1 7.0 6.6 6.5...
    6.3 6.2 6.4 6.0 6.0 5.7 5.6 5.9 5.5 5.6 5.4 5.5 5.2 5.5 5.2 5.3 5.4...
    4.9 5.3 5.0 5.3 4.9 4.9 5.2 4.9 4.8 5.0 4.7 4.7 4.5 4.8 4.7 4.6 4.8...
    4.5 4.3 4.3 4.5 4.6 4.2 4.2 4.3 4.0 4.0 4.3 4.3 4.4 3.9 4.1 4.1 3.9...
    4.2 3.8 3.8 4.0 4.0 3.7]*1e-3;

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_11mSm_NonTreated = t;
V_11mSm_NonTreated = V;

figure('Name', 'Voltage Monitoring. 11.4 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t(1:53),V(1:53)*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',...
    colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 3;
t02 = 37;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 11.4 mS/m KCl Non-Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
factor_11mSm = 7.489392706815978e+03;
zetaP_11mSm_NonTreated = eta/(eps*E^2)*factor_11mSm*sp;
Error_zetaP_11mSm_NonTreated = eta/(eps*E^2)*factor_11mSm*Esp;

%% Experimental zeta potential estimation data for treated surfaces

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 1.5 mS/m
V = [19.03 18.72 18.62 18.78 18.54 18.37 18.62 18.26 18.24 18.53 18.61...
    18.14 18.18 18.51 18.28 18.23 18.41 18.26 18.04 18.06 18.35 18.36...
    18.02 18.23 18.34 18.29 18.05 18.03 18.09 18.14 18.04 18.20 17.93...
    18.05 18.14 18.06 17.91 18.10 18.05 18.03 17.67 18.20 17.82 18.05...
    17.84 17.82 17.94 18.03 17.78 17.94 17.73 17.86 17.89 18.00 18.11...
    17.52 17.56 17.74 17.84 17.79 17.74 17.63 17.85 17.87 17.96 17.86...
    17.65 17.74 17.64 17.80 17.56 17.78 17.65 17.64 17.63 17.86 17.63...
    17.68 17.92 17.56 17.77 17.77 17.87 17.85 17.72 17.63 18.04 17.86...
    17.93]*1e-3;

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_1mSm_Treated = t;
V_1mSm_Treated = V;

figure('Name', 'Voltage Monitoring. 1.5 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t(5:end),V(5:end)*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',...
    colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 1;
t02 = 50;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 1.5 mS/m KCl Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
factor_1mSm = 1.733608523735045e+04; 
zetaP_1mSm_Treated = eta/(eps*E^2)*factor_1mSm*sp;
Error_zetaP_1mSm_Treated = eta/(eps*E^2)*factor_1mSm*Esp;

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 5.2 mS/m
V = [45.86 45.63 45.66 45.48 45.31 45.17 45.19 45.20 45.10 45.07 44.98...
    45.11 44.93 44.92 44.90 45.13 44.83 45.02 44.76 44.87 44.77 44.78...
    44.93 44.77 44.67 44.74 44.62 44.64 44.77 44.67 44.67 44.80 44.69...
    44.61 44.55 44.51 44.62 44.64 44.44 44.49 44.59 44.38 44.59 44.56...
    44.38 44.55 44.25 44.39 44.31 44.46 44.40 44.11 44.27 44.35 44.15...
    44.32 44.22 44.10 44.04 44.16 44.14 43.98 44.13 43.99 44.00 43.91...
    44.05 44.08 43.90 43.99 43.93 43.90 44.05 43.97 43.90 43.90 43.88...
    43.96 43.81 43.90 43.87 44.15 44.11 43.99 44.02 43.91 44.01 43.99...
    43.85 43.87 44.04 43.95 44.09 43.86 43.99 43.88 43.82 43.94 43.92]*...
    1e-3;

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_5mSm_Treated = t;
V_5mSm_Treated = V;

figure('Name', 'Voltage Monitoring. 5.2 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t,V*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 7;
t02 = 60;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 5.2 mS/m KCl Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
ratio = 906.4916;
zetaP_5mSm_Treated = (1/(1-0.9369))*eta*sp/(eps*E^2)*ratio;
Error_zetaP_5mSm_Treated = (1/(1-0.9369))*eta*Esp/(eps*E^2)*ratio;

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 11.4 mS/m
V = [11.1 11.3 11.3 11.5 11.5 11.3 11.5 11.4 11.7 11.5 11.5 11.7 11.7...
    11.3 11.6 11.5 11.7 11.7 11.4 11.5 11.7 11.8 11.4 11.4 11.6 11.1...
    11.4 11.4 11.5 11.5 11.2 11.4 11.5 11.1 11.2 11.4 11.3 11.0 10.9...
    11.1 11.3 11.3 10.9 11.2 11.0 11.3 11.2 11.1 10.8 11.0 11.1 10.8...
    11.1 10.7 10.6 11.0 10.8 10.6 11.0 10.5 10.7 10.5 10.7 10.7 10.6...
    10.8 10.8 10.4 10.6 10.6 10.3 10.5 10.7 10.3 10.6 10.3 10.5 10.2...
    10.3 10.4 10.5 10.2 10.3 10.3 10.3 10.1 10.3 9.9 10.2 10.0 10.0...
    10.1 9.7 10 10 10 9.9 9.8 9.6 9.7 9.8 9.8 9.6 9.7 9.6 9.5 9.2 9.8...
    9.4 9.4 9.7 9.4 9.6 9.6 9.2 9.5 9.5 9.6 9.3 9.3 9.4 8.9 9.2 9.1 9.3...
    9.3 9.0 9.2 9.2 9.0 9.1 9.2 8.9 9.2 9.0 8.8 9 9 8.8 8.7 8.9 8.8 8.7...
    9.0 8.7 8.8 8.8 8.9 8.5 8.5 8.5 8.7 8.5 8.6 8.6 8.7 8.4 8.7 8.5 8.4...
    8.5 8.6 8.5 8.5 8.4 8.3 8.5 8.5 8.2 8.4 8.4 8.0 8.3 8.5 8.3 8.1 8.0...
    8.1 8.0 8.2 8 7.8 8 7.9 8 8.2 8.1 7.8 7.9 7.7 7.7 7.9 8 8 8 7.6 7.7...
    7.6 8.0 7.7 7.9 7.9 7.8 7.5 7.7 7.7 7.8 7.3 7.6 7.5 7.4 7.7 7.2 7.5...
    7.4 7.6 7.3 7.2 7.5 7.5 7.3 7.6 7.1 7.3 7.1 7.2 7.3 7.4 7.3 7.3 7.3...
    7.5 7.3 7.2 7.2 7.3 7 7.3 7 7.1 7.1 7.3 7.2 7 7 7.2 7.3 6.9 7.1 7.1...
    6.8 7 6.8 6.9 7.0 6.7 7.0 6.8 6.9 6.8 7.1 7 6.7 6.8 7 6.5 6.8 6.7...
    6.7 6.8 6.3 6.8 6.8 6.7 6.8 6.6 6.7 6.6 6.5 6.4 6.6 6.5 6.6 6.5 6.7...
    6.4 6.6 6.7 6.3 6.7 6.5 6.3 6.7 6.3 6.5 6.3 6.5 6.4 6.4 6.4 6.2 6.3...
    6.4 6.4 6.5 6.2 6.3 6.1 6.3 6.5 6.0 6.4 6.1 6.2 6.4 6.3 6.0 6.4 5.9...
    6.2 6.1 6.2 6 6.1 6.3 6.3 5.8 5.9 5.9 5.9 6.2 6.1 5.9 6.2 5.8 6.0...
    6.2 6 6 6.1 5.8 6.1 5.6 5.8 5.9 5.7 6 6 5.8 5.8 6.2 5.7 5.8 6.1 6 ...
    6.1 6 5.7 5.7 5.9 5.6 5.8 5.7 5.7 5.9 5.5 5.8]*1e-3;

N = length(V);
t = 0:1:(N-1); % Measurements taken every 1 second
t = t';
t_11mSm_Treated = t;
V_11mSm_Treated = V;

figure('Name', 'Voltage Monitoring. 11.4 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
plot(t,V*1e3,'ob','MarkerSize',10,'LineWidth',2, 'Color',colorB)

% Start and end time used for interpolation based on the error estimation
% (shown in last section)
t01 = 20;
t02 = 120;
pV = polyfit(t(t01:t02),V(t01:t02),1);
z = [t01,t02];
y = z*pV(1)+pV(2);
hold on
plot(z,y*1e3,'r','LineWidth',3, 'Color',colorR)
xlabel('t(s)'); ylabel('(mV)')
title('Voltage Monitoring. 11.4 mS/m KCl Treated')

sp = pV(1);
rP = corr(t(t01:t02),V(t01:t02));
Esp = Deltap(sp,rP,length(t(t01:t02)));
factor_11mSm = 7.489392706815978e+03;
zetaP_11mSm_Treated = eta/(eps*E^2)*factor_11mSm*sp;
Error_zetaP_11mSm_Treated = eta/(eps*E^2)*factor_11mSm*Esp;

%% Error Estimation for the obtained values. Non-Treated Surfaces

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 1.5 mS/m
Ni = 22; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 23-jj;
    t2 = 24+jj;
    pVtemp = polyfit(t_1mSm_NonTreated(t1:t2),V_1mSm_NonTreated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_1mSm_NonTreated(t1:t2),V_1mSm_NonTreated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_1mSm_NonTreated(t1:t2)));
    zeta(jj) = eta/(eps*E^2)*factor_1mSm*pV(jj);
    Ezeta(jj) = eta/(eps*E^2)*factor_1mSm*DpV(jj);
end

figure('Name','Error estimation. 1.5 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 7;
t02 = 40;
h1=errorbar(t02-t01,zetaP_1mSm_NonTreated*1e3,...
    Error_zetaP_1mSm_NonTreated*1e3,'or','MarkerFaceColor',colorR,...
    'Color',colorR);
axis([0,55,-120,-60])
xlabel('Number of points taken for interpolation around t = 23.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 1.5 mS/m KCl Non-Treated')

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 5.2 mS/m
Ni = 24; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 25-jj;
    t2 = 26+jj;
    pVtemp = polyfit(t_5mSm_NonTreated(t1:t2),V_5mSm_NonTreated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_5mSm_NonTreated(t1:t2),V_5mSm_NonTreated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_5mSm_NonTreated(t1:t2)));
    zeta(jj) = 20*eta*pVtemp(1)/(eps*E^2)*ratio;
    Ezeta(jj) = 20*eta*DpV(jj)/(eps*E^2)*ratio;
end

figure('Name','Error estimation. 5.2 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 6;
t02 = 45;
h1=errorbar(t02-t01,zetaP_5mSm_NonTreated*1e3,...
    Error_zetaP_5mSm_NonTreated*1e3,'or','MarkerFaceColor',colorR,...
    'Color',colorR);
axis([0,50,-90,-50])
xlabel('Number of points taken for interpolation around t = 25.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 5.2 mS/m KCl Non-Treated')

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 11.4 mS/m
Ni = 19; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 20-jj;
    t2 = 21+jj;
    pVtemp = polyfit(t_11mSm_NonTreated(t1:t2),...
        V_11mSm_NonTreated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_11mSm_NonTreated(t1:t2),V_11mSm_NonTreated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_11mSm_NonTreated(t1:t2)));
    zeta(jj) = eta/(eps*E^2)*factor_11mSm*pV(jj);
    Ezeta(jj) = eta/(eps*E^2)*factor_11mSm*DpV(jj);
end

figure('Name','Error estimation. 11.4 mS/m KCl Non-Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 3;
t02 = 37;
h1=errorbar(t02-t01,zetaP_11mSm_NonTreated*1e3,...
    Error_zetaP_11mSm_NonTreated*1e3,'or','MarkerFaceColor',colorR,...
    'Color',colorR);
axis([0,45,-90,-60])
xlabel('Number of points taken for interpolation around t = 20.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 11.4 mS/m KCl Non-Treated')

%% Error Estimation for the obtained values. Treated Surfaces

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 1.5 mS/m
Ni = 23; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 24-jj;
    t2 = 25+jj;
    pVtemp = polyfit(t_1mSm_Treated(t1:t2),V_1mSm_Treated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_1mSm_Treated(t1:t2),V_1mSm_Treated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_1mSm_Treated(t1:t2)));
    zeta(jj) = eta/(eps*E^2)*factor_1mSm*pV(jj);
    Ezeta(jj) = eta/(eps*E^2)*factor_1mSm*DpV(jj);
end

figure('Name','Error estimation. 1.5 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 1;
t02 = 50;
h1=errorbar(t02-t01,zetaP_1mSm_Treated*1e3,Error_zetaP_1mSm_Treated...
    *1e3,'or','MarkerFaceColor',colorR,'Color',colorR);
axis([0,55,-40,0])
xlabel('Number of points taken for interpolation around t = 24.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 1.5 mS/m KCl Treated')

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 5.2 mS/m
Ni = 27; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 33-jj;
    t2 = 34+jj;
    pVtemp = polyfit(t_5mSm_Treated(t1:t2),V_5mSm_Treated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_5mSm_Treated(t1:t2),V_5mSm_Treated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_5mSm_Treated(t1:t2)));
    zeta(jj) = (1/(1-0.9369))*eta*pVtemp(1)/(eps*E^2)*ratio;
    Ezeta(jj) = (1/(1-0.9369))*eta*DpV(jj)/(eps*E^2)*ratio;
end

figure('Name','Error estimation. 5.2 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 7;
t02 = 60;
h1=errorbar(t02-t01,zetaP_5mSm_Treated*1e3,Error_zetaP_5mSm_Treated*1e3,...
    'or','MarkerFaceColor',colorR,'Color',colorR);
axis([0,60,-30,-10])
xlabel('Number of points taken for interpolation around t = 33.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 5.2 mS/m KCl Treated')

% -------------------------------------------------------------------------
% KCl Electrolyte conductivity of 11.4 mS/m
Ni = 65; % Number of different slopes
pV = zeros(Ni,1); % Vector containing slopes
DpV = zeros(Ni,1); % Vector containing error of slopes
zeta = zeros(Ni,1); % Vector containing results of zeta potential
Ezeta = zeros(Ni,1); % Vector containing results of zeta potential error

for jj=1:Ni
    t1 = 70-jj;
    t2 = 71+jj;
    pVtemp = polyfit(t_11mSm_Treated(t1:t2),V_11mSm_Treated(t1:t2),1);
    pV(jj) = pVtemp(1);
    r = abs(corr(t_11mSm_Treated(t1:t2),V_11mSm_Treated(t1:t2)));
    DpV(jj) = Deltap(pV(jj),r,length(V_11mSm_Treated(t1:t2)));
    zeta(jj) = eta/(eps*E^2)*factor_11mSm*pV(jj);
    Ezeta(jj) = eta/(eps*E^2)*factor_11mSm*DpV(jj);
end

figure('Name','Error estimation. 11.4 mS/m KCl Treated')
hold on
box on
ax=gca;
ax.FontSize = 20;
ax.FontName = 'Times New Roman';
NEl = 3:2:2*Ni+2;
errorbar(NEl,zeta*1e3, Ezeta*1e3,'ob','MarkerFaceColor',colorB,'Color',...
    colorB)
t01 = 20;
t02 = 120;
h1=errorbar(t02-t01,zetaP_11mSm_Treated*1e3,Error_zetaP_11mSm_Treated...
    *1e3,'or','MarkerFaceColor',colorR,'Color',colorR);
axis([0,140,-25,-5])
xlabel('Number of points taken for interpolation around t = 70.5 s')
ylabel('zeta potential (mV)')
legend(h1,'Value used in Report')
title('Error estimation. 11.4 mS/m KCl Treated')

%% Functions for error estimation
% Error estimation of the slope
function [Dp]=Deltap(p,r,N)
Dp = sqrt(p^2*(r^(-2)-1)/(N-2));
end

% Correlation coefficient
function [r]=corr(x,y)
N = length(x);
sumaproducto = 0;
sumax2 = 0;
sumay2 = 0;

for j=1:N
    sumaproducto = sumaproducto+x(j)*y(j);
    sumax2 = sumax2+x(j)^2;
    sumay2 = sumay2+y(j)^2;
end

divisor1 = sqrt(sumax2-N*mean(x)^2);
divisor2 = sqrt(sumay2-N*mean(y)^2);
r = (sumaproducto-N*mean(x)*mean(y))/(divisor1*divisor2);
end