%% MSE's script - Treating Chronoamperometry-EPR data
%% 2021 - 2022 - Roessler Group - Imperial College London
%%
% Chornoamp-EPR : double integral of EPR spectra

% With this script you can run data in folders "Fig6_S15" and "S14_S16"
%%%%%%

clear all
close all

%%---------EPR 

% Change path or comment it if you save this script in the same folder
% where your data are
path(path,'/Users/j/Documents/2021_Imperial/Projects/FE-EPR/Paper/Data/Fig6_S15/EPR');

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

% 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);
xlabel('Time (s)');
ylabel('Profil of EPR signals');

%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);
xlabel('Magnetic field (mT)');
ylabel('EPR signal (a.u.)');


% Choose the region to integrate on the EPR spectra directly on figure 2
% Choose the lowest field position then click shift to hold the Datatip 
% then choose the highest field position
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);
xlabel('Magnetic field (mT)');
ylabel('Double integration of EPR signal (a.u.)');

% 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)*10^6; %uM

figure(4)
scatter(x,Conc);
xlabel('Time (s)');
ylabel('Concentration of STEMPO (uM)');

rawdata_03 = [x Conc.'];
path='/Users/j/Documents/2021_Imperial/Projects/FE-EPR/Paper/Data/Fig6_S15/EPR/';

filename =[path 'rawdata_ChronoampEPR_03.xlsx'];
writematrix(rawdata_03,filename,'Sheet',1)
