% Representation of the experimental dataset for 
% "Concentration-Polarization Electroosmosis Near Insulating Constrictions
%   within Microfluidic Channels"

close all;clear;clc

%% General parameters

% Measured frequencies
f = [30 44 65 94 137 200]; % Hz
% Microscope pixel-micron conversion
conv = 4.025; % px/fr

%% Measurements 1.7 mS/m

% 30 Hz. Pixels per frame
vM_30Hz_1 = mean([3.0376 2.7751 2.7774 2.2570]);
std_30Hz_1 = std([3.0376 2.7751 2.7774 2.2570]);

% 44 Hz. Pixels per frame
vM_44Hz_1 = mean([2.5477 2.0106 2.4265 2.2131]);
std_44Hz_1 = std([2.5477 2.0106 2.4265 2.2131]);

% 65 Hz. Pixels per frame
vM_65Hz_1 = mean([2.3234 2.0543 2.8374 2.4402]);
std_65Hz_1 = std([2.3234 2.0543 2.8374 2.4402]);

% 94 Hz. Pixels per frame
vM_94Hz_1 = mean([2.2319 2.2583 1.9990 1.6754]);
std_94Hz_1 = std([2.2319 2.2583 1.9990 1.6754]);

% 137 Hz. Pixels per frame
vM_137Hz_1 = mean([1.9655 1.8381 1.7088 1.7558]);
std_137Hz_1 = std([1.9655 1.8381 1.7088 1.7558]);

% 200 Hz. Pixels per frame
vM_200Hz_1 = mean([1.6361 1.9362 1.4231 1.3611]);
std_200Hz_1 = std([1.6361 1.9362 1.4231 1.3611]);

% Gathering of data. Pixels per frame
pxfr_1 = [vM_30Hz_1 vM_44Hz_1 vM_65Hz_1 vM_94Hz_1...
    vM_137Hz_1 vM_200Hz_1];
std_1 = [std_30Hz_1 std_44Hz_1 std_65Hz_1 std_94Hz_1...
    std_137Hz_1 std_200Hz_1];

% Conversion of data to SI
v_1 = pxfr_1*conv*1e-6; % m/s
std_1 = std_1*conv*1e-6; % m/s

%% Third Set 6.1 mS/m

% 30 Hz. Pixels per frame
vM_30Hz_6 = mean([1.2897 0.89817 0.98986]);
std_30Hz_6 = std([1.2897 0.89817 0.98986]);

% 44 Hz. Pixels per frame
vM_44Hz_6 = mean([1.0789 1.1396 0.84538 0.78205]);
std_44Hz_6 = std([1.0789 1.1396 0.84538 0.78205]);

% 65 Hz. Pixels per frame
vM_65Hz_6 = mean([1.0193 0.96855 0.8220 0.90247]);
std_65Hz_6 = std([1.0193 0.96855 0.8220 0.90247]);

% 94 Hz. Pixels per frame
vM_94Hz_6 = mean([0.8073 0.77941 0.75343 0.66789]);
std_94Hz_6 = std([0.8073 0.77941 0.75343 0.66789]);

% 137 Hz. Pixels per frame
vM_137Hz_6 = mean([0.82968 0.68868 0.62682 0.50486]);
std_137Hz_6 = std([0.82968 0.68868 0.62682 0.50486]);

% 200 Hz. Pixels per frame
vM_200Hz_6 = mean([0.69229 0.58417]);
std_200Hz_6 = std([0.69229 0.58417]);

% Gathering of data. Pixels per frame
pxfr_6 = [vM_30Hz_6 vM_44Hz_6 vM_65Hz_6 vM_94Hz_6 vM_137Hz_6 vM_200Hz_6];
std_6 = [std_30Hz_6 std_44Hz_6 std_65Hz_6 std_94Hz_6 std_137Hz_6...
    std_200Hz_6];

% Conversion of data to SI
v_6 = pxfr_6*conv*1e-6; % m/s
std_6 = std_6*conv*1e-6; % m/s

%% Third Set 12.2 mS/m

% 30 Hz. Pixels per frame
vM_30Hz_12 = mean([0.64083 0.84399 0.78898]);
std_30Hz_12 = std([0.64083 0.84399 0.78898]);

% 44 Hz. Pixels per frame
vM_44Hz_12 = mean([0.73181 0.78349]);
std_44Hz_12 = std([0.73181 0.78349]);

% 65 Hz. Pixels per frame
vM_65Hz_12 = mean([0.85732 0.74582 0.76557]);
std_65Hz_12 = std([0.85732 0.74582 0.76557]);

% 94 Hz. Pixels per frame
vM_94Hz_12 = mean([0.83860 0.84542 0.47365 0.47888]);
std_94Hz_12 = std([0.83860 0.84542 0.47365 0.47888]);

% 137 Hz. Pixels per frame
vM_137Hz_12 = mean([0.53765 0.56817 0.3645 0.43093]);
std_137Hz_12 = std([0.53765 0.56817 0.3645 0.43093]);

% 200 Hz. Pixels per frame
vM_200Hz_12 = mean([0.45985 0.3585]);
std_200Hz_12 = std([0.45985 0.3585]);

% Gathering of data. Pixels per frame
pxfr_12 = [vM_30Hz_12 vM_44Hz_12 vM_65Hz_12 vM_94Hz_12...
    vM_137Hz_12 vM_200Hz_12];
std_12 = [std_30Hz_12 std_44Hz_12 std_65Hz_12 std_94Hz_12...
    std_137Hz_12 std_200Hz_12];

% Conversion of data to SI
v_12 = pxfr_12*conv*1e-6; % m/s
std_12 = std_12*conv*1e-6; % m/s

%% Figure 4: Experimental data
figure('Name','Triangular Constrictions. Experimental values')
hold on
box on
axis([25 240 1 20])
ax=gca;
ax.FontSize = 20;
ax.LineWidth = 1.5;
ax.FontName = 'Times New Roman';
set(gca, 'XScale','log', 'YScale','log')
xlabel('\fontname{Times New Roman} \fontsize{24} f (Hz)')
ylabel('\fontname{Times New Roman} \fontsize{24} Mean Velocity Magnitude (\mum/s)')

% Data conversion to microns
v_1 = v_1/1e-6; std_1=std_1/1e-6;
v_6 = v_6/1e-6; std_6=std_6/1e-6;
v_12 = v_12/1e-6; std_12=std_12/1e-6;

h1 = errorbar(f(1:end),v_1(1:end),std_1(1:end),'ok',...
    'LineWidth',1.5,'MarkerSize',10);

zf1 = linspace(f(1),f(end));
p = polyfit(log(f(1:end)),log(v_1(1:end)),1);
PF1 = p(1);
yE1 = exp(p(1)*log(zf1)+p(2));
colorR = [236,176,181]/255;
plot(zf1,yE1,'--k','LineWidth',2,'color',colorR);

h2 = errorbar(f(1:end),v_6(1:end),std_6(1:end),'sk',...
    'LineWidth',1.5,'MarkerSize',10);

zf1 = linspace(f(1),f(end));
p = polyfit(log(f(1:end)),log(v_6(1:end)),1);
PF6 = p(1);
yE1 = exp(p(1)*log(zf1)+p(2));
colorG = [195,245,199]/255;
plot(zf1,yE1,'--k','LineWidth',2,'color',colorG);

h3 = errorbar(f(1:end),v_12(1:end),std_12(1:end),'^k',...
    'LineWidth',1.5,'MarkerSize',10);

zf1 = linspace(f(1),f(end));
p = polyfit(log(f(1:end)),log(v_12(1:end)),1);
PF12 = p(1);
yE1 = exp(p(1)*log(zf1)+p(2));
colorB = [194,224,252]/255;
plot(zf1,yE1,'--k','LineWidth',2,'color',colorB);

legend([h1 h2 h3],'1.7 mS/m','6.1 mS/m','12.2 mS/m','Location','southwest')