% Script to load in several CV measurements (Ivium potentiostat)
% Determine peak heights and plot against scan rate (diffusion vs. surface
% bound)

% from low potential to high potential

% Jana Eisermann and Yunfei Dang

clear, clc, close all

%% Load Excel Sheet

path(path,'/Users/j/Documents/2021_Imperial/Projects/FE-EPR/STEMPO/STEMPO/Electrochem/210824_Sputtered-elec/');
          
filename = {'CV_eleccell_24 5mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_23 10mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_22 20mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_21 30mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_20 40mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_19 50mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_18 60mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_17 70mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_16 80mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_15 90mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_01 100Vps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_02 125mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_03 150mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_04 175mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_05 200mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_06 225mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_07 250mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_08 275mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_09 300mVps pH8 500mMNa2CO3 AgAgCl.xlsx',...
             'CV_eleccell_10 325mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_11 350mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_12 375mVps pH8 500mMNa2CO3 AgAgCl.xlsx',... 
             'CV_eleccell_13 400mVps pH8 500mMNa2CO3 AgAgCl.xlsx'};

             
sample = length(filename);
scan_rate = [5,10,20,30,40,50,60,70,80,90,100,125,150,175,200,225,250,275,300,325,350,375,400]; % write down used scan rate in order in mV/s
scan_rate_2 = sqrt(scan_rate); % calculate square root of scan rate
log_scan_rate = log10(scan_rate/1000); % calculate log of scan rate in V/s

%% Read out raw data

for n = 1:sample
    [Voltage{n},Current{n}] = CVload_mA_SHE(filename{1,n},2); % read file,import potential V and current A, choose scan No.2
    CV_y{n} = Current{n};
    CV_x{n} = Voltage{n};
end


% Plot CV measurements

figure(1)
box on
hold on
for n = 1:sample
    plot(CV_x{n},CV_y{n},'linewidth',1.0);
end
axis tight
xlabel('Potential (V) vs SHE')
ylabel('Current (mA)')
% xlim([-0.3 1.7])
% ylim([-1.1 1.2])

set(gca,'FontSize',18)
title('CV of immobilised STEMPO on ITO','Fontsize',18);

% lgd = cell(1,sample); % predefine size 
for n = 1:sample
    lgd{n} = sprintf('%d mV/s', scan_rate(1,n)); % adjusts legend handles based on data set
end
legend(lgd)

%% Potential positioning

% separate CV into forward and back scan

for n = 1:sample
    data_1{n}(:,1) = CV_x{n}(1:round(length(CV_x{n})/2)); % forward scan x-values (voltage)
    data_1{n}(:,2) = CV_y{n}(1:round(length(CV_y{n})/2)); % forward scan y-values (current)
    
    data_2{n}(:,1) = CV_x{n}(round(length(CV_x{n})/2+1):length(CV_x{n})); % back scan x-values
    data_2{n}(:,2) = CV_y{n}(round(length(CV_y{n})/2+1):length(CV_y{n})); % back scan y-values
    
    data_length(n) = length(data_1{n}); % read out number of points for forward/back scan
end


% maximum current
for n = 1:sample
    cut_max(n) = round(0.2*data_length(n)); % data point to 'cut' search for local minimum (%-value of total data length)
    [y(n),x(n)] = max(data_1{n}(1:data_length(n)-cut_max(n),2)); % determine maximum current (y) and data point (x)
    max_volt(n) = data_1{n}(x(n),1); % determine voltage value for data point x
end


% minimum current
for n = 1:sample
    cut(n) = round(0.2*data_length(n)); % data point to 'cut' search for local minimum (%-value of total data length)
    [yy(n),xx(n)] = min(data_2{n}(1:data_length(n)-cut(n),2)); % determine minimum current (yy) and data point(xx); 1:data_length(n)-cut(n) defines range for search
    min_volt(n) = data_2{n}(xx(n),1); % determine voltage value for data point xx
end

% E_mid calculation

for n = 1:sample

    E_mid(n) = 0.5*(max_volt(n)+min_volt(n)); % midpoint redox potential
    
end

E0 = mean(E_mid); % the average E

%% Trumpet plot

figure(2)
for n = 1:sample
delta_E_A(n) = max_volt(n) - E_mid(n);
delta_E_C(n) = min_volt(n) - E_mid(n);
end

% plot trumpet plot
hold on
Sca_E_A = scatter(log_scan_rate,delta_E_A,80,'r','filled');
Sca_E_C = scatter (log_scan_rate,delta_E_C,80,'k','filled');

% looking for delta_E larger than 200 mV (0.2 V)
delta_E = delta_E_A - delta_E_C; % calculate peak separation (delta_E)
E_length = length(delta_E);
k = find(delta_E > 0.2,1); % find out the minimun delta_E > 200mV

% a loop to browse through all delta_E, to check if it > 200 mV
% i = 1;
% while(i<(E_length + 1)) 
%     if delta_E(i) > 0.2
%         delta_E_lowlimit = delta_E(i);
%         break  
%     end
%     i=i+1;
% end


% anodic linear fitting
P = polyfit(log_scan_rate(k:E_length),delta_E_A((k:E_length)),1);
x1 = ((log_scan_rate(k)-0.5):0.001:(log_scan_rate(E_length)));
delta_E_A_fit = P(1)*x1 + P(2);
plot (x1,delta_E_A_fit,'r-.', 'linewidth',1.5);

% cathodic linear fitting
Q = polyfit(log_scan_rate(k:E_length),delta_E_C((k:E_length)),1);
delta_E_C_fit = Q(1)*x1 + Q(2);
plot (x1,delta_E_C_fit,'k-.', 'linewidth',1.5);

% draw the zero base line
x2 = (-3:0.001:-0.50);
y2 = x2 .* 0;
plot (x2,y2,'b--', 'linewidth',1.0);
set(gca,'FontSize',18)
box on
% xlim([-3 -0.48])
% ylim([-0.5 0.5])
xlabel('log(\nu) (V/s)')
ylabel('\DeltaE_{p} (V)')

legend([Sca_E_A Sca_E_C],{'Anodic peaks','Cathodic peaks'},'Location','northwest','Fontsize',18) 

% critical scan rate
length_x1 = length(x1);

% calculate average critical scan rate
for n = 1:length_x1
    if delta_E_A_fit(n) > delta_E_C_fit(n)   % find the data point E_A > E_C
        log_v_c = x1(n);
        y_log_v_c = delta_E_A_fit(n);
        break
    end
end
v_c = 10^log_v_c*1000;
% text(log_v_c-0.2,y_log_v_c * 8,['v_{c}=' num2str(v_c) ' V/s'],'FontSize',14,'color','b')
text(log_v_c-0.2,y_log_v_c * 8,sprintf('v_{c}= %0.3f V/s',v_c),'FontSize',14,'color','b');

% Laviron analysis
% Anodic
[v_c_A,one_alpha,kapp_A] = Laviron(P(1),P(2));
% text(log_v_c-0.2 * (5),y_log_v_c * (-8),['v_{c,A}=' num2str(v_c_A) ' V/s; 1-a=' num2str(one_alpha) '; k_{app,A}=' num2str(kapp_A) ' s-1'],'FontSize',14,'color','r')
text(log_v_c-0.2 * (5),y_log_v_c * (-8),sprintf('v_{c,A}= %0.3f V/s; 1-a= %0.3f; k_{app,A}= %0.3f  s^{-1}',v_c_A,one_alpha,kapp_A),'FontSize',14,'color','r')

% Cathodic
[v_c_C,alpha,kapp_C] = Laviron(Q(1),Q(2));
% text(log_v_c-0.2 * (5),y_log_v_c * (-11),['v_{c,C}=' num2str(v_c_C) ' V/s; a=' num2str(alpha) '; k_{app,C}=' num2str(kapp_C) ' s-1'],'FontSize',14,'color','k')
text(log_v_c-0.2 * (5),y_log_v_c * (-11),sprintf('v_{c,C}= %0.3f V/s; a= %0.3f; k_{app,C}= %0.3f  s^{-1}',v_c_C,alpha,kapp_C),'FontSize',14,'color','k')


%% Remove data with scan rate > critical scan rate

number = find(scan_rate > v_c,1);
number = number -1;

scan_rate = scan_rate(1:number);
scan_rate_2 = scan_rate_2(1:number);
data_1 = data_1(1:number);
data_2 = data_2(1:number);
sample = number;

%% Peak height determination
% data 1 (forward scan)

for n = 1:sample
    
    edge_1_left(n) = round(0.25*data_length(n)); % define step width from the maximum data point x (%-value of total data length)
    edge_1_right(n) = round(0.25*data_length(n)); % define step width from the maximum data point x (%-value of total data length)
    % pick points for 'baseline'
    off_volt_1(n) = data_1{n}(x(n)-edge_1_left(n),1); % voltage value of left point; position defined via x(n)-edge_1_left(n)
    off_cur_1(n) = data_1{n}(x(n)-edge_1_left(n),2); % fitting current to off_volt_1
    ende_volt_1(n) = data_1{n}(x(n)+edge_1_right(n),1); % voltage value of right point; position defined via x(n)+edge_1_right(n)
    ende_cur_1(n) = data_1{n}(x(n)+edge_1_right(n),2); % fitting current to ende_volt_1
    % combine to calculate parameters for baseline
    A(:,n) = [off_volt_1(n) ende_volt_1(n)]; % x-values for both offset and ende point
    B(:,n) = [off_cur_1(n) ende_cur_1(n)]; % y-values for both offset and ende point
    C(n,:) = [[1; 1]  A(:,n)]\B(:,n); % calculate baseline
    slope_1(n) = C(n,2); % read out slope of baseline
    intercept_1(n) = C(n,1); % read out intercept of baseline
    % apply baseline parameters
    x_base_1{n} = data_1{n}(x(n)-edge_1_left(n):x(n)+edge_1_right(n),1); % define range for x-values (by using 'x(n)-edge_1_left(n):x(n)+edge_1_right(n)')
    y_base_1{n} = C(n,2)*x_base_1{n}+C(n,1); % calculate y-values in the given x_base range
    y_cross_1(n) = C(n,2)*max_volt(n)+C(n,1); % determine crossing point from maximum voltage point and baseline
    cur_diff_1(n) = y(n) - y_cross_1(n); % calculate Peak height

end

% data 2 - apply similar routine to data 1 set
% keep in mind: data point 1 for back scan is at the highest measured
% potential

for n = 1:sample
    
    edge_2_left(n) = round(0.25*data_length(n));    
    edge_2_right(n) = round(0.25*data_length(n));
    % pick points for 'baseline'
    off_volt_2(n) = data_2{n}(xx(n)+edge_2_left(n),1); % xx(n)+edge_2_left(n): left data point for baseline
    off_cur_2(n) = data_2{n}(xx(n)+edge_2_left(n),2);
    ende_volt_2(n) = data_2{n}(xx(n)-edge_2_right(n),1); % xx(n)-edge_2_right(n): left data point for baseline
    ende_cur_2(n) = data_2{n}(xx(n)-edge_2_right(n),2);
    % combine to calculate parameters for baseline    
    AA(:,n) = [off_volt_2(n) ende_volt_2(n)];
    BB(:,n) = [off_cur_2(n) ende_cur_2(n)];
    CC(n,:) = [[1; 1]  AA(:,n)]\BB(:,n);
    slope_2(n) = CC(n,2);
    intercept_2(n) = CC(n,1);    
    % apply baseline parameters
    x_base_2{n} = data_2{n}(xx(n)-edge_2_right(n):xx(n)+edge_2_left(n),1);
    y_base_2{n} = CC(n,2)*x_base_2{n}+CC(n,1);
    y_cross_2(n) = CC(n,2)*min_volt(n)+CC(n,1);
    cur_diff_2(n) = yy(n) - y_cross_2(n);

end

% Plot to check automatic picking 
% choose number of subplot rows and column according to n
nrows = floor(sqrt(n));
ncols = ceil(n/nrows);

figure(3)
set(gca,'FontSize',18)
for subp = 1:n
    subplot(nrows, ncols, subp); % define the layout of the subplots based on the number of measurements
    
    box on
    hold on
    plot(data_1{subp}(:,1),data_1{subp}(:,2),'r') % plot forward scan
    plot(data_2{subp}(:,1),data_2{subp}(:,2),'k') % plot back scan
    
    plot(x_base_1{subp},y_base_1{subp},'r') % plot baseline for forward scan
%     xline(max_volt(1,subp),'r') % vertical line crossing maximum voltage
    xline(max_volt(subp),'r') % vertical line crossing maximum voltage

    plot(x_base_2{subp},y_base_2{subp},'k') % plot baseline for back scan
%     xline(min_volt(1,subp),'k') % vertical line crossing minimum voltage
    xline(min_volt(subp),'k') % vertical line crossing minimum voltage
    
    axis tight
    title(sprintf('scan rate %d mV/s', scan_rate(1,subp))); % adjust title automatically based on measurement
    xlabel('Voltage in V')
    ylabel('Current in mA')

end




% %% Linear Fits for Scan rate plots
% %mdl-command used because it calculates Rsq-value 
% 
mdl_1 = fitlm(scan_rate,cur_diff_1); % fitting routine forward scan (vs. v)
Rsq_1 = mdl_1.Rsquared.Ordinary; % quality of fit

mdl_11 = fitlm(scan_rate,cur_diff_2); % fitting routine back scan (vs. v)
Rsq_11 = mdl_11.Rsquared.Ordinary; % quality of fit

mdl_2 = fitlm(scan_rate_2,cur_diff_1); % fitting routine forward scan (vs. v^(1/2))
Rsq_2 = mdl_2.Rsquared.Ordinary; % quality of fit

mdl_22 = fitlm(scan_rate_2,cur_diff_2); % fitting routine forward scan (vs. v^(1/2))
Rsq_22 = mdl_22.Rsquared.Ordinary; % quality of fit
% 
% Scan Rate Plots
figure(4)
set(gca,'FontSize',18)
subplot(1,2,1)
box on
hold on 
plot(scan_rate,cur_diff_1,'kx','markersize',16)
plot(scan_rate,cur_diff_2,'rx','markersize',16)
% plot(scan_rate,mdl_1.Fitted,'k')
% plot(scan_rate,mdl_11.Fitted,'r')

mdl_x = (0:0.01:max(scan_rate));
mdl_1_a = table2array(mdl_1.Coefficients(2,1));
mdl_1_b = table2array(mdl_1.Coefficients(1,1));
mdl_1_y = mdl_x * mdl_1_a + mdl_1_b;

mdl_11_a = table2array(mdl_11.Coefficients(2,1));
mdl_11_b = table2array(mdl_11.Coefficients(1,1));
mdl_11_y = mdl_x * mdl_11_a + mdl_11_b;

plot(mdl_x,mdl_1_y,'k','linewidth',1)
plot(mdl_x,mdl_11_y,'r','linewidth',1)
set(gca,'FontSize',18)
axis tight
xlabel('\nu (mV/s)')
ylabel('Peak current (mA)')
tex = text(scan_rate(1,2),0.9*cur_diff_1(1,sample),['R^2 = ' num2str(Rsq_1)]);
tex = text(scan_rate(1,2),0.1*cur_diff_1(1,sample),['R^2 = ' num2str(Rsq_11)],'color','r');

subplot(1,2,2)
box on
hold on 
plot(scan_rate_2,cur_diff_1,'kx','markersize',16)
plot(scan_rate_2,cur_diff_2,'rx','markersize',16)
plot(scan_rate_2,mdl_2.Fitted,'k')
plot(scan_rate_2,mdl_22.Fitted,'r')
set(gca,'FontSize',18)
axis tight
xlabel('\nu^{1/2} (mV/s)^{1/2}')
ylabel('Peak current (A)')
tex = text(scan_rate_2(1,2),0.9*cur_diff_1(1,sample),['R^2 = ' num2str(Rsq_2)]);
tex = text(scan_rate_2(1,2),0.1*cur_diff_1(1,sample),['R^2 = ' num2str(Rsq_22)],'color','r');