% Based on Equation (4) from Forgez et al. 2010

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Thermal inputs
Rout = 6.651562529;
cp = 978.2969842;
Rin = 1.408390014;
m = 0.095;

CpRinRout = m*cp*(Rin+Rout);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find the length of time of the experiment
tfin = length(texp);

% Calculate the open circuit voltage at each time, this is based on the
% state of charge of the cell
% Initially predefine a vector
OCVatt = zeros(tfin,1);
% Go through each time
for i = 1:tfin
    % Find the location of the state of charge
    [~, idx] = min(abs(SoC-(100*SoCexp(i))));
    % Use this location to find the open circuit voltage
    OCVatt(i) = OCV(idx);
end

% Calculate the average value for the ambient temperature
TaAverage = mean(Taexp);

% Calculate the heat generation at each time
%Qdot = abs(Iexp.*(Vexp-OCVatt));

Qdot = abs(Iexp.*(Vexp-OCVatt));;


% Predefine a vector for the surface temperature
Ts = zeros(tfin,1);

% The initial temperature is the same as the initial temperature from the 
% experiment
Ts(1) = Tsexp(1);

E_heat = 0;
E_heat_exp = 0;

% Go through each time
for i = 2:tfin
    
    % Calculate the surface temperature
    Ts(i) = (((TaAverage-Ts(i-1))/CpRinRout + (Qdot(i)*Rout)/CpRinRout) * (texp(i) - texp(i-1))) + Ts(i-1);
    
    E_heat = E_heat + (Ts(i) - TaAverage) / Rout;
    E_heat_exp = E_heat_exp + (Tsexp(i) - TaAverage) / Rout;
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE PLOTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 12)

figure('units','centimeters','position',[10 10 8.8 5.0])

hold on

subplot(2,1,1);

p1 = plot(texp/60,Ts,'b',texp/60,Tsexp,'g');

xlabel('Time (mins)')
ylabel(sprintf('Temperature (%cC)', char(176)))
legend('Calculated Temperature','Measured Temperature','Location','northwest')

subplot(2,1,2);
p3 = plot(texp/60,abs(Tsexp-Ts));
%p3 = plot(texp/60,100*abs((Tsexp-Ts)./(Tsexp+273)));
p3.Color = [0 0 1.0000];
p3.LineStyle = '-';
p3.LineWidth = 1;
p3.Marker = 'none';

xlabel('Time (mins)')
ylabel(sprintf('Absolute Error (%cC)', char(176)))
%ylabel(sprintf('Percentage Error (%)', char(176)))


%legend('Calculated Temperature','Measured Temperature','northwest')