%% Maryam Seif-Eddine's FE-EPR script for analysing CV-EPR data
%% 2021-2022 - Roessler Research Group - Imperial College London
%% Please cite publication (Seif-Eddine et al. Nature Chemistry 2024) if using this script

clc
clear all
close all

%%

%%---------EPR - data---------- 
 
%%
 
%Change path
%path(path,'/Users/j/Documents/2021_Imperial/Projects/FE-EPR/STEMPO/STEMPO/EPR/July22/10_07'); 

%Change file name
[B,spec,Params] = eprload('CVEPR_STEMPO_02');

%Components
TimeEPR = B{1,2}; % is the total time of programmed EPR increments
Y = B{1,1}*1e-1; % magnetic field  
Z = spec; % data

%Clean the EPR matrix from the zeros
zerocol = find(sum(abs(Z)) == 0,1);% will find the 0s in the matrix and return the index of the fst 0  
X = reshape(1:zerocol,1,[])'; % nb of increments as X axis
Xaxis = X(1:zerocol-1); % to limit the time to the real time of acquisition and not the equivalent time of the number of increments added
EPRdata = spec(1:size(Z,1), 1:zerocol-1); % to limit the data to the real time of aquisition
TimeEPR = TimeEPR(1:zerocol+1);  % total time of measurement only

%Baseline Correction
EPRdata = basecorr(EPRdata,1,2); % 2d order polynomial

%Smooth data - if neede -
EPRdata =  datasmooth(EPRdata,20,'savgol');

%EPR part of the first figure
figure(1)
subplot(2,1,1)
z=peaks(TimeEPR,Y);
ribbon(Y,EPRdata,0.5,'Color','black');
%Format of the figure
xlabel('EPR increments');
hxLabel = get(gca,'xlabel');
set(hxLabel,'rotation',7,'VerticalAlignment','middle');
axis tight
ylabel('Magnetic field (mT)');
hyLabel = get(gca,'ylabel');
set(hyLabel,'rotation',-10,'VerticalAlignment','middle');
zlabel('Intensity (a.u.)');
set (gca, 'YDir','reverse'); 
set(gca,'FontSize',12)

shading interp;

%%

 %%---------Electrochemistry - data---------

 
%%

%Change path
%path(path,'/Users/j/Documents/2021_Imperial/Projects/FE-EPR/STEMPO/STEMPO/Electrochem/10_07_22_20mMBA'); 

%Change file name
filename = 'CVEPR_STEMPO_02 nosub O2 flow4 5mVps 0p2to1V 1sc'; 

% Read out raw data
raw = readtable(filename,'HeaderLines',1); % read the complete .txt file without the header
we_potential = raw{:,1}+0.211; % potential of working electrode in V
current = raw{:,3}; % measured current in A
TimeCV = raw{:,2};
scan_nb = raw{:,4};

totscan = max(raw{:,4}); % number of performed scans
points = length(we_potential); % length of the raw table
length_matrix = points/totscan; % length of each CV scan

%Re-classify data: each scan in a separate cell
for n = 1:totscan 
    scan{n} = [we_potential(1+(n-1)*length_matrix:n*length_matrix) current(1+(n-1)*length_matrix:n*length_matrix) TimeCV(1+(n-1)*length_matrix:n*length_matrix)];
end

%CV data in the first figure: plotting all the scans
figure(1)
subplot(2,1,2)
box on
hold on
for n = 1:totscan
   plot(scan{n}(:,1),scan{n}(:,2)*1e6,'linewidth',2.0,'Color','black')
end
%Format of the figure
axis tight
xlabel('Potential (V vs SHE)')
set(gca,'XMinorTick','on')
ylabel('Current (µA)')
set(gca,'FontSize',12)

%%
    % Preparing for combining EPR and Echem data
    
%%

%Deduce the echem waiting time from the echem time vector
TimeCV = TimeCV(:,1)-TimeCV(1,1); 

%Define the time at the end of the CV experiment (including all scans)
endCV = round(TimeCV(end));

%Find the equivalent time in the EPR time vector
reallengthEPR = interp1(TimeEPR, 1:length(TimeEPR), endCV, 'nearest');
TimeEPR = TimeEPR(1:reallengthEPR+1);

%Asking the user to enter the number of CV scan to be analysed
x = input('Enter the CV scan number: '); 

%Narrowing down the data sets to the scan number provided by the user
we_potentialx = scan{x}(:,1); 
currentx = scan{x}(:,2);
TimeCVscanx = scan{x}(:,3);
TimeCVscanx = TimeCVscanx-TimeCVscanx(1,1);

%Finding the equivalent data set in the EPR matrix
lengthonescanEPR = length(TimeEPR)/totscan;
lengthonescanEPR = round(lengthonescanEPR);

if x == 1
TimeaxisEPR = TimeEPR(1 : x*lengthonescanEPR);
aa = 1;
else
aa = (x-1)*lengthonescanEPR;
bb = x*lengthonescanEPR-2*x;
TimeaxisEPR = TimeEPR(aa : bb);
end

bb = x*lengthonescanEPR-2*x;
EPRdatax = EPRdata(:,aa:bb); % to limit the data to the wanted CV scan
%%
%--------------------------------------------------------------------------
                               % Definitions
%--------------------------------------------------------------------------

%Maximum and miniumu of the delimited data sets from above
ii = length(we_potentialx)/2;
we_potentialx2 = we_potentialx(ii:length(we_potentialx));
[Min,Imin] = min(we_potentialx);
[Min2,Imin2] = min(we_potentialx2);
Imin2 = 2*Imin2-2;
[Max,Imax] = max(we_potentialx);

%Start potential and its time
start_pot = we_potential(1,1);
start_pot = round(start_pot);
tCVi = TimeCV(1,1);

%Lowest potential and its time
lowest_pot = Min;
lowest_pot = round(lowest_pot);
tCV_lowest_pot_1 = TimeCVscanx(Imin);

%Highest potential and its time
highest_pot = Max;
tCV_highest_pot_1 = TimeCVscanx(Imax)
%%
%--------------------------------------------------------------------------
              % Condition: IF start potential = lowest potential
%--------------------------------------------------------------------------

% start potential = lowest potential means we start from and stop at the
% same potential
if start_pot == lowest_pot

% l1: 1st lowest potential
tCV_l1 = TimeCVscanx(1,1);
indx_pot_CV_l1 = find(TimeCVscanx == tCV_l1);

% h1: 1st highest potential
tCV_h1 = tCV_highest_pot_1;
indx_pot_CV_h1 = find(TimeCVscanx == tCV_h1);

% l2: 2d lowest potential
indx_scnd_lowest_pot = Imin2;
tCV_l2 = TimeCVscanx(indx_scnd_lowest_pot);

%Forward CV scan data: from the lowest to the highest potential
tCV_l1h1 = TimeCVscanx(indx_pot_CV_l1:indx_pot_CV_h1);
Pot_l1h1 = we_potentialx(indx_pot_CV_l1:indx_pot_CV_h1);
Current_l1h1 = currentx(indx_pot_CV_l1:indx_pot_CV_h1);
forward = [tCV_l1h1;Pot_l1h1;Current_l1h1];

indx_l1 = interp1(TimeaxisEPR, 1:length(TimeaxisEPR), tCV_l1, 'nearest');
tEPR_l1 = TimeaxisEPR(indx_l1);
indx_h1 = interp1(TimeaxisEPR, 1:length(TimeaxisEPR), tCV_h1, 'nearest');
tEPR_h1 = TimeaxisEPR(indx_h1);

tEPR_l1h1 = TimeaxisEPR(indx_l1:indx_h1);
EPRdatax_l1h1 = EPRdatax(indx_l1:indx_h1); 

%Forward CV scan - Corresponding EPR data
EPRdatax_l1h1 = EPRdatax(:,indx_l1:indx_h1); 
EPRdatax_l1h1_bkg = basecorr(EPRdatax_l1h1,1,2);
field= Y.';

%Smoothing EPR data - if needed -
EPRdatax_l1h1_bkg =  datasmooth(EPRdatax_l1h1_bkg,20,'savgol');

%Figure of the EPR spectra corresponding to the forward CV scan only
figure('DefaultAxesFontSize',12)
[m,n] = size(EPRdatax_l1h1_bkg);
figure(2)
cc = turbo(n);
for i=1:n
    plot(field,EPRdatax_l1h1_bkg(:,i),'Color',cc(i,:),'LineWidth',1);
    hold on
    i=i+1;
    axis tight
xlabel('Magnetic Field (mT)');
ylabel('EPR intensity (a.u.)');
end
title('Forward CV scan');
figure('DefaultAxesFontSize',18);

% Choose the region to integrate on the EPR spectra directly on figure 2
disp('Choose the lowest field position, then click & hold shift before choosing the highest field position')
disp('! final outcome dependent on chosen posisions !')
dcm_obj = datacursormode(figure(2));
set(dcm_obj,'DisplayStyle','datatip',...
'SnapToDataVertex','off','Enable','on')
pause();
Datatip = getCursorInfo(dcm_obj);
p1 = extractfield(Datatip,'Position');

pmax=(p1(1,1));
pmin=(p1(1,3));

field= Y.';
dfield = mean(diff(field));

% Find the indexes of the chosen field positions
indx_pmax = interp1(field,1:length(field),pmax,'nearest');
indx_pmin = interp1(field,1:length(field),pmin,'nearest');

% Limiting the EPR data to the chosen field positions
signal_selc = field(indx_pmin:indx_pmax); %region of selected signal
EPRdata_for_bkg1 = EPRdatax_l1h1_bkg(indx_pmin:indx_pmax,:);

[f,g]=size(EPRdata_for_bkg1);

% Integrate then double integrate the corresponding data
for i = 1:g
    
EPRdata_for_bkg_int(:,i) = cumtrapz(EPRdata_for_bkg1(:,i))*dfield; %int for integral

EPRdata_for_bkg_dint(:,i) = sum(EPRdata_for_bkg_int(:,i))*dfield; %dint for double integral

end

% Represent the integrated data in a figure
figure(3)
cc = turbo(n);
for i=1:n
    plot(signal_selc,EPRdata_for_bkg_int(:,i),'Color',cc(i,:),'LineWidth',1);
    hold on
    i=i+1;
    axis tight
xlabel('Magnetic Field (mT)');
ylabel('Integral EPR intensity (a.u.)');

end
title('Forward CV scan');
figure('DefaultAxesFontSize',18);

intensities_for = EPRdata_for_bkg_dint.'; % Double integrals of EPR on th forward scan are now intensities_for

%Finding the potential of each EPR data point (EPR spectrum)
ratio_for = length (Pot_l1h1)/g; %pot for potential. count steps that correspond to EPR on Potential axis
indx_x_for = 1:ratio_for:length(Pot_l1h1); %find indexes of the steps above on the potential axis
indx_x_for=round(indx_x_for);
x_for =(Pot_l1h1(indx_x_for));%find the corresponding potential at each index from above. x_for: potential points from forward CV that correspond to each EPR on the forward CV


%%
%Reverse CV scan data: from the highest to the 2d lowest (which is
%equivalent to the 1st lowest under the condition where star potential =
%stop potential

indx_l2  = length(TimeaxisEPR);
tEPR_l2 = TimeaxisEPR(indx_l2);

tCV_h1l2 = TimeCVscanx(indx_pot_CV_h1:indx_scnd_lowest_pot);
Pot_h1l2 = we_potentialx(indx_pot_CV_h1:indx_scnd_lowest_pot);
Current_h1l2 = currentx(indx_pot_CV_h1:indx_scnd_lowest_pot);

tEPR_h1l2 = TimeaxisEPR(indx_h1:indx_l2-1);
EPRdatax_h1l2 = EPRdatax(indx_h1:indx_l2-1);


%Reverse CV scan - Corresponding EPR data
EPRdatax_h1l2 = EPRdatax(:,indx_h1:indx_l2-2);
EPRdatax_h1l2_bkg = basecorr(EPRdatax_h1l2,1,2);

%Smoothing EPR data - if needed -
EPRdatax_h1l2_bkg =  datasmooth(EPRdatax_h1l2_bkg,20,'savgol');

%Figure of the EPR spectra corresponding to the reverse CV scan only
figure('DefaultAxesFontSize',12)
figure(4)
[m,n] = size(EPRdatax_h1l2_bkg);
cc = turbo(n);
for i=1:n
    plot(field,EPRdatax_h1l2_bkg(:,i),'Color',cc(i,:),'LineWidth',1);
    hold on
    i=i+1;
    axis tight
xlabel('Magnetic Field (mT)');
ylabel('EPR intensity (a.u.)');

end
title('Reverse CV scan');
figure('DefaultAxesFontSize',12);

% Choose the region to integrate on the EPR spectra directly on figure 4
disp('REPEAT: Choose the lowest field position, then click & hold shift before choosing the highest field position')
dcm_obj = datacursormode(figure(4));
set(dcm_obj,'DisplayStyle','datatip',...
'SnapToDataVertex','off','Enable','on')
pause();
Datatip = getCursorInfo(dcm_obj);
p1 = extractfield(Datatip,'Position');

pmax=(p1(1,1));
pmin=(p1(1,3));

field= Y.';
dfield = mean(diff(field));

% Find the indexes of the chosen field positions
indx_pmax = interp1(field,1:length(field),pmax,'nearest');
indx_pmin = interp1(field,1:length(field),pmin,'nearest');

% Limiting the EPR data to the chosen field positions
signal_selc = field(indx_pmin:indx_pmax);
EPRdata_rev_bkg1 = EPRdatax_h1l2_bkg(indx_pmin:indx_pmax,:);

[f,g]=size(EPRdata_rev_bkg1);

% Integrate then double integrate the corresponding data
for i = 1:g
    
EPRdata_rev_bkg_int(:,i) = cumtrapz(EPRdata_rev_bkg1(:,i))*dfield;

EPRdata_rev_bkg_dint(:,i) = sum(EPRdata_rev_bkg_int(:,i))*dfield;

end

% Represent the integrated data in a figure
figure(5)
cc = turbo(n);
for i=1:n
    plot(signal_selc,EPRdata_rev_bkg_int(:,i),'Color',cc(i,:),'LineWidth',1);
    hold on
    i=i+1;
    axis tight
xlabel('Magnetic Field (mT)');
ylabel('Integral EPR intensity (a.u.)');

end
title('Reverse CV scan');
figure('DefaultAxesFontSize',12);

intensities_rev = EPRdata_rev_bkg_dint.'; % Double integrals of EPR on th reverse scan are now intensities_rev

%Finding the potential of each EPR data point (EPR spectrum)
[Min3,Imin3] = min(Pot_h1l2);
length_Poth1l2 = length(Pot_h1l2);
ratio_rev =(length (Pot_h1l2)/g);
indx_x_rev = 1:ratio_rev:(length_Poth1l2)-1;
x_rev = Pot_h1l2(round(indx_x_rev));

end
%%
%-----------------------------------
        %Nernst fitting
%-----------------------------------

% forward
% Nernst equation

Ex_for = x_for;  
time_for = tEPR_l1h1(1:length(intensities_for));
Potential_for = Ex_for;

scanrate = 0.005;% scan rate V/s

% Error on the potential of the EPR data
for i = 1:length(time_for)-1
dtfor(i,1) =  (time_for(i+1,1) - time_for(i,1));

dPotential_for(i,1) = scanrate * dtfor(i,1);
end
dPotential_for(length(time_for),1) = dPotential_for(length(time_for)-1,1);

figure('DefaultAxesFontSize',12)

figure(6)
subplot (1,2,1);
scatter(Potential_for,intensities_for);
% errorbar(Potential_for,intensities_for,dPotential_for,'horizontal','LineStyle','none');


xlim([0.4 1.2]);
xlabel('Applied potential (V vs SHE)');
title ('Forward');
ylim([-0.02 0.45]);
ylabel('Double integral EPR - forward CV scan');

hold all

rng default % for reproducibility
x = Potential_for;
xlow = x(x<0.5);
xhigh = x(x>0.5);

y = intensities_for;
ylow = y(x<0.5);
yhigh = y(x>0.5);

hhigh=plot(xhigh,yhigh);

hhigh.LineWidth = 2;


start = [0.3; 0.951];
options = optimset('TolX',1e-5);
[CfornoxEox,residualfor,rrsfor] = fminsearch (@(CfornoxEox) fitfor(CfornoxEox,xhigh,yhigh,hhigh),start,options);

%%
% Reverse
% Nernst equation

% 
Ex_rev = x_rev;  
Potential_rev = Ex_rev;
time_rev = tEPR_h1l2(1:length(intensities_rev));

scanrate = 0.005;% scan rate V/s

% Error on the potential of each EPR spectrum
for i = 1:length(time_rev)-1
dtrev(i,1) =  (time_rev(i+1,1) - time_rev(i,1));

dPotential_rev(i,1) = scanrate * dtrev(i,1);
end
dPotential_rev(length(time_rev),1) = dPotential_rev(length(time_rev)-1,1);

figure('DefaultAxesFontSize',12)

figure(6)
subplot (1,2,2);
Potential_Rev = Ex_rev;

scatter(Potential_Rev,intensities_rev);
% errorbar(Potential_Rev,intensities_rev,dPotential_rev,'horizontal','LineStyle','none');

hold on

Potential_Rev = Ex_rev;

xlim([0.4 1.2]);
xlabel('Applied potential (V vs SHE)');
title ('Reverse');
ylim([-0.02 0.45]);
ylabel('Double integral EPR - reverse CV scan');

hold all

% Simplex - Derivative-free method using fminsearch
rng default % for reproducibility
x = Potential_rev;
y = intensities_rev;
xhigh = x(x>0);
yhigh = y(x>0);

hhigh=plot(xhigh,yhigh);

hhigh.LineWidth = 2;

start = [0.3; 0.951];
options = optimset('TolX',1e-4);

[CrevnoxEoxrd,residualrevrd,rrsrevrd] = fminsearch (@(CrevnoxErd) fitrevrd(CrevnoxErd,xhigh,yhigh,hhigh),start,options);

set (gca, 'XDir','reverse');

close(figure(7))
close(figure(8))
close(figure(9))