%% Import data from text file
% Script for importing data from the following text file:
%   filename: C:\Users\LAB\Desktop\combined_test\potentiodynamic_test_E_vs_I.cor
% Auto-generated by MATLAB on 03-Sep-2024 10:42:39
%% Set up the Import Options and import the data

opts = delimitedTextImportOptions("NumVariables", 3);

% Specify range and delimiter
opts.DataLines = [26, Inf];
opts.Delimiter = "\t";

% Specify column names and types
opts.VariableNames = ["E", "i", "t"];
opts.VariableTypes = ["double", "double", "double"];

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Import the data
potentiodynamictestEvsI = readtable("C:\Users\LAB\Desktop\combined_test\potentiodynamic_test_E_vs_I.cor", opts);
clear opts;

%
%% Clear temporary variables

potentiodynamictestEvsI(:,'logi') = log(abs(potentiodynamictestEvsI(:,'i')))



% anodic and current indices
cathodic_index = table2array(potentiodynamictestEvsI(:,'i') < 0);
anodic_index = table2array(potentiodynamictestEvsI(:,'i') > 0);



% Create plot of potentiodynamictestEvsI.logi and potentiodynamictestEvsI.E


h2 = plot( ...
    potentiodynamictestEvsI(cathodic_index ,:).logi, ...
    potentiodynamictestEvsI(cathodic_index ,:).E, ...
    "DisplayName","cathodic");
hold on;
h1 = plot( ...
    potentiodynamictestEvsI(anodic_index, :).logi, ...
    potentiodynamictestEvsI(anodic_index, :).E, ...
    "DisplayName","anodic");

anodic_index = table2array(potentiodynamictestEvsI(:,'i') < 0);


% Add xlabel, ylabel, title, and legend
xlabel("logi")
ylabel("E")
title("E vs. logi")
legend

hold off;

%

E = table2array(potentiodynamictestEvsI(:,"E"));
i = table2array(smoothdata(potentiodynamictestEvsI(:,"i"),"loess",100));

% anodic and cathodic indices
cathodic_index = table2array(potentiodynamictestEvsI(:,'i') < 0);
anodic_index = table2array(potentiodynamictestEvsI(:,'i') > 0);

cathodic_E = table2array(potentiodynamictestEvsI(cathodic_index,"E"));
cathodic_i = table2array(smoothdata(potentiodynamictestEvsI(cathodic_index,"i"),"loess",100));

anodic_E = table2array(potentiodynamictestEvsI(anodic_index,"E"));
anodic_i = table2array(smoothdata(potentiodynamictestEvsI(anodic_index,"i"),"loess",100));

E_OCP = -0.4615

[params, fittedCurve] = fitButlerVolmer(...
    E - E_OCP, ...
    i ...
);

plot( ...
    log(abs(fittedCurve(:,2))), ...
    fittedCurve(:,1) + E_OCP, ...
    "DisplayName","fittedcurve");
hold on;
plot( ...
    table2array(potentiodynamictestEvsI(:,"logi")), ... 
    table2array(potentiodynamictestEvsI(:,"E")), ...    
    "DisplayName","original");
hold off;

legend("Position", [0.14167,0.12103,0.19107,0.082143])

%
% Fit Butler-Volmer curve
%

function [params, fittedCurve] = fitButlerVolmer(overPotential, currentDensity)
% FITBUTLERVOLMER Fits the Butler-Volmer equation to Tafel data
%   [PARAMS, FITTEDCURVE] = FITBUTLERVOLMER(TAFELDATA) takes Tafel 
%   data in the form of a two-column matrix (overpotential in the first 
%   column, current density in the second) and fits the Butler-Volmer 
%   equation to it. The function returns the fitted parameters in the 
%   structure PARAMS and the fitted curve as a two-column matrix 
%   FITTEDCURVE. 

% Define the Butler-Volmer equation as a function handle
bvEqn = @(params, overpotential) params(1) * (exp(params(2) * overpotential) - exp(-params(3) * overpotential));

% Set initial parameter estimates (adjust if necessary)
initialParams = [1, 0.5, 0.5]; 

% Perform the nonlinear least squares fitting
params = lsqcurvefit(bvEqn, initialParams, overPotential, currentDensity);

% Generate the fitted curve
fittedCurve = [overPotential, bvEqn(params, overPotential)];

% Display the fitted parameters (optional)
fprintf('Fitted Parameters:\n');
fprintf('  Exchange Current Density (i0): %f\n', params(1));
fprintf('  Anodic Transfer Coefficient (alpha_a): %f\n', params(2));
fprintf('  Cathodic Transfer Coefficient (alpha_c): %f\n', params(3));
end