%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)