%% Maryam Seif-Eddine's FE-EPR script for treating Chronoamperometry-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

%%
%%---------Chornoamp-EPR : double integral of EPR spectra 

% Change path
%path(path,'/Users/j/Documents/2021_Imperial/Projects/FE-EPR/STEMPO/STEMPO/EPR/Apr22/24_04');

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

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

%Clean the EPR matrix from the zeros
[A,limindex]=find((~Z),1); % will find the 0s in the matrix and return the index of the fst 0  
X = reshape(1:limindex,1,[])'; % nb of increments as X axis
Xaxis = X(1:limindex-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:limindex-1); % to limit the data to the real time of aquisition

TimeEPR = TimeEPR(1:limindex+1);  % total time of measurement only

%Baseline Correction
EPRdata_bkg = basecorr(EPRdata,1,2);

%Smooth data - if needed -
EPRdata_bkg =  datasmooth(EPRdata_bkg,50,'savgol');

%EPR figure; profile vs time
figure('DefaultAxesFontSize',14)
figure(1)
[m,n] = size(EPRdata);
ratio = length (TimeEPR)/n;
indx_x = 1:ratio:length(TimeEPR);
indx_x=round(indx_x);
x =(TimeEPR(indx_x));
plot(x, EPRdata_bkg);

%Limiting figure by colors
figure(2)

indx_x1 = interp1(x,1:length(x),100,'nearest');
indx_x2 = interp1(x,1:length(x),200,'nearest');

EPRdata_bkg_x1 = EPRdata_bkg(:,1:indx_x1);
EPRdata_bkg_x2 = EPRdata_bkg(:,indx_x1:indx_x2);
EPRdata_bkg_x3 = EPRdata_bkg(:,indx_x2:n);

plot(Y,EPRdata_bkg_x1,'Color','r','LineWidth',1);
hold on
plot(Y,EPRdata_bkg_x2,'Color','b','LineWidth',1);
hold on
plot(Y,EPRdata_bkg_x3,'Color','r','LineWidth',1);

% 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.';
dB = mean(diff((Y)));

% 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_bkg1 = EPRdata_bkg(indx_pmin:indx_pmax,:);

% Integrating EPR spectra
EPRdata_bkg_int = cumtrapz(EPRdata_bkg1)*dB;

figure(3)

EPRdata_bkg_int1 = EPRdata_bkg_int(:,1:indx_x1);
EPRdata_bkg_int2 = EPRdata_bkg_int(:,indx_x1:indx_x2);
EPRdata_bkg_int3 = EPRdata_bkg_int(:,indx_x2:n);

plot(signal_selc,EPRdata_bkg_int1,'Color','r','LineWidth',1);
hold on
plot(signal_selc,EPRdata_bkg_int2,'Color','b','LineWidth',1);
hold on
plot(signal_selc,EPRdata_bkg_int3,'Color','r','LineWidth',1);

% Double integrating EPR spectra
EPRdata_bkg_dint = cumtrapz(EPRdata_bkg_int)*dB;

% Findingthe maximum of each double integration for quantification
EPRdata_dint_max = EPRdata_bkg_dint(end,:);

% Equaion of the calibration curve realised with different concentrations of 4-aminoTEMPO
% y = 8.5 * x (mM)

Conc  = ((EPRdata_dint_max)/8525); %M

figure(4)
scatter(x,Conc);