%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the average temperature of the last 1000s
TAvLast1000 = mean(Tsexp(end-1000:end));
tCpRinRout = round(find(Tsexp>TAvLast1000-(TAvLast1000-Taexp(1))*0.05,1)/10);
% Calculate the length of the vectors
tfin = find(texp==tCpRinRout);
% Calculate the average value for the ambient temperature
TaAverage = mean(Taexp(1:tfin));
% Calculate the heat generation at each time
Qdot = abs(Iexp.*(Vexp-OCV));
% Predefine a vector for the least squares error, the value for CpRinRout
% is varied between 1 and 3000
Err = zeros(3000,1);
% Go through each value of CpRinRout
for CpRinRout = 1:3000
% 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);
% Go through each timestep
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);
% Calculate the difference between the modelled surface temperature
% (Ts(i)) and the measured surface temperature (Tsexp(i)) and square
% the difference for the least squared solution. Add this to the
% error already obtained. Do this for every time (i) and for every
% value of CpRinRout
Err(CpRinRout) = Err(CpRinRout) + (Ts(i) - Tsexp(i))^2;
end
end
% The best value for CpRiRout is the value with the least error
[~,CpRinRoutLS] = min(Err)