clear all;
clc;

%% Loading data to workspace
load TestData;
load RCData;
load OCVData;

%% Initialize simulation
Q0 = 6.05; %Cell initial capacity
Qnom = 7.6; %Cell nominal capacity

%% SOC breakpoints for R0, R1, R2, C1, C2
SOC_breakpoints = [0 10 20 30 40 50 60 70 80 90 100];

%% Thermal model constants
Cth = 273; %Thermal capacity
Rcond = 1.5; %Conduction resistor
Rconv = 3.6; %Convection resistor
Tamb = 25; %Ambient temperature

%% Run simulation
sim("Model_class.slx")

%% Visualization

% Electrical model visualization
figure(1);
subplot(2,2,1)
plot(SIM_time,SIM_Vbat_experimental,SIM_time,SIM_Vbat);
xlabel('time(s)')
ylabel(' Voltage [V]')
xlim([0 Time_test(364409)])
legend('Experimental voltage','Simulation voltage')
title('Simulation and experimental voltage comparison')
subplot(2,2,3)
plot(SIM_time,SIM_Ibat);
xlabel('time(s)')
ylabel(' Current [A]')
xlim([0 Time_test(364409)])
title('Battery current')
subplot(2,2,2)
plot(SIM_time,SIM_SOC);
xlabel('time(s)')
ylabel(' SOC [%]')
xlim([0 Time_test(364409)])
title('Estimated SOC')
subplot(2,2,4)
plot(SIM_time,SIM_Error_percentage);
xlabel('time(s)')
ylabel(' Error [%]')
xlim([0 Time_test(364409)])
title('Estimated error for voltage')
grid

% Thermal model visualization
SIM_Tamb = SIM_Tamb*ones(length(SIM_time),1);
figure(2);
subplot(2,2,1)
plot(SIM_time,SIM_Tin,SIM_time,SIM_Ts,SIM_time,SIM_Tamb);
xlabel('time(s)')
ylabel(' Temperature [ºC]')
xlim([0 Time_test(364409)])
legend('Tinside','Tsurface','Tambient')
title('Temperature comparison figure')
subplot(2,2,3)
plot(SIM_time,SIM_Ibat);
xlabel('time(s)')
ylabel(' Current [A]')
xlim([0 Time_test(364409)])
title('Battery current')
subplot(2,2,2)
plot(SIM_time,SIM_OCV,SIM_time,SIM_Vbat);
xlabel('time(s)')
ylabel(' Voltage [V]')
xlim([0 Time_test(364409)])
title('OCV and Vbat comparison')
subplot(2,2,4)
plot(SIM_time,SIM_qt);
xlabel('time(s)')
ylabel(' Heat source [W]')
xlim([0 Time_test(364409)])
title('Internal heat generation within the cell')
grid

figure(3);
plot(SIM_time,8*SIM_Vbat);
xlabel('time(s)')
ylabel(' Voltage [V]')
xlim([0 16000])

figure(4);
plot(SIM_time,SIM_Ibat);
xlabel('time(s)')
ylabel(' Current [A]')
xlim([0 16000])
ylim([-15 1])