

%clear all
clc
clf

plotte = input('Plot figure Toutside = 1\n states = 2\n c1 = 3\n c2 = 4 \n sum = 5\n c3 = 6\n c4 = 7');  

%% Parameters in use in the differential equations

a.vent = 0.05;             %[m^3/s] 75
a.room = 100;            %[m^3], 806
a.U = 0.03;               %[kW/m^2*K]
a.A = 1.8;                %[m^2]
a.Tw = 273 ;             %[K]
a.massys = 0.07;         %[kg/s] 19.7
a.radiator = 0.060;       %[m^3]
a.fus= 334;             %[kJ/kg]    \cite{11}
a.cpair = 1.005;       %[kJ/kg K], \cite{11}, 273.15 K
a.densair = 1.2929;     %[kg/m^3], \cite{11}, 0 [C]
a.cpwater = 4.220;      %[kJ/kg K], \cite{11}, 273.15 K
a.denswater = 999.87 ;     %[[kg/m^3]\cite{11}
 
%% Outside temperature

var.time = [0 7 12 15 17 20 24]*3600'; %[h]
var.signals.values = [20 23 30 30 25 23 20]'+273;
var.signals.dimensions = [1]';

if plotte == 1 

figure(1);

plot(var.time/3600,var.signals.values-273)
title('Outdoor Temperature')
xlabel('Time [h]')
ylabel('Outdoor temperature [deg C]')

end
                      
%% Case switch

%normal values
alpha1 = 0;
alpha2 = 0;
Qchiller = 0;
Qchiller2 = 0;
Qchiller3 = 0;
tswitch1 = 36000;       %10h
tswitch2 = 54000;       %15h
tswitch3 = 72000;       %20h
tswitch4 = 0;           %0h
tswitch5 = 0;

switch 4
%case one
    case 1
%case two    
    case 2
        alpha1 = 1.0;
%case three
    case 3
        alpha1 = 0.5;
        alpha2 = -0.5;
%case four
    case 4
        alpha1 = .5;
        alpha2 = -.5;
        Qchiller = 0.5;
        Qchiller2 = -0.5;
        %Qchiller3 = -0.5;
        %tswitch1 = 18000;   %5h
        %tswitch2 = 36000;   %10h
        tswitch3 = 61200;    %17h
        tswitch4 = 79200;    %22h
        
        
end


%% Simulink simulation

tsim = 3600*24;     %24 hour simulation over one day

sim('master190413')


%% simulation values

t = simout.time/3600;
Troom = simout.signals.values(:,1)-273;
Tradiator = simout.signals.values(:,2)-273;
Mass_ice = simout.signals.values(:,3);
Mass_water = simout.signals.values(:,4);
T_outdoor = simin.signals.values(:,1)-273;
split = simin.signals.values(:,2);
Q_chiller = simin.signals.values(:,3);

%% Plots of the different states

set(0,'DefaultAxesFontSize',14)

if plotte == 2

figure(2);

subplot(4,1,1);
plot(t,Troom)
axis([0 24 10 32]);
title('Temp. room')
%xlabel('Time [h]')
ylabel('Temp. room [\circC]')

subplot(4,1,2);
plot(t,Tradiator)
axis([0 24 0 32]);
title('Temp. radiator')
%xlabel('Time [h]')
ylabel('Temp. radiator [\circC]')

subplot(4,1,3);
plot(t,Mass_ice)
axis([0 24 15 60]);
title('Mass ice')
%xlabel('Time [h]')
ylabel('Mass ice [kg]')

subplot(4,1,4);
plot(t,Mass_water)
axis([0 24 70 125]);
title('Mass water')
xlabel('Time [h]')
ylabel('Mass water [kg]')

end

%% Case one

% Change in room and radiator temperature with split = 0.

if plotte == 3 
    
figure(3);

subplot(3,1,1)
plot(t,T_outdoor)
axis([0 24 18 32]);
title('Outdoor temp.')
xlabel('Time [h]')
ylabel('Outdoor temp. [\circC]')

subplot(3,1,2);
plot(t,Troom)
axis([0 24 15 32]);
set(gca,'YTick',[15 20 25 30])
title('Temp. room')
xlabel('Time [h]')
ylabel('Temp. room [\circC]')

subplot(3,1,3);
plot(t,Tradiator)
axis([0 24 14 30]);
title('Temp. radiator')
xlabel('Time [h]')
ylabel('Temp. radiator [\circC]')
    

end


%% case two

% change in room and raditor temperature when step in split at t=10h
% step in split from 0 to 1.

if plotte == 4
    
figure(4);

subplot(3,1,1);
plot(t,Troom)
axis([0 24 10 25]);
title('Temp. room')
xlabel('Time [h]')
ylabel('Temp. room [\circC]')


subplot(3,1,2)
plot(t,Tradiator)
axis([0 24 0 25]);
title('Temp. radiator')
xlabel('Time [h]')
ylabel('Temp. radiator [\circC]')


subplot(3,1,3)
plot(t,split)
axis([0 24 0 1]);
title('Split fraction')
xlabel('Time [h]')
ylabel('\alpha')

end


%% Case three

% Changes in mass_water and mass_ice during cooling.
% Parameters: split1 = 0.5 at time = 10 h, and split2 = -0.5 at time = 15 h
% Q_chiller = 0


if plotte == 5
 
%The sum of the mass_ice and mass_water.

Tank = sum(simout.signals.values(:,3:4),2) 

figure(5)

plot(t,Tank)
axis([0 24 120 140]);
title('Total mass');
xlabel('Time [h]');
ylabel('Total mass [kg]');

end 

if plotte == 6
    
figure(6);

subplot(3,1,1);
plot(t,Mass_ice)
axis([0 24 1 60]);
set(gca,'YTick',[10 20 30 40 50 60])
title('Mass ice')
xlabel('Time [h]')
ylabel('Mass ice [kg]')

subplot(3,1,2);
plot(t,Mass_water)
axis([0 24 70 150]);
title('Mass water')
xlabel('Time [h]')
ylabel('Mass water [kg]')

subplot(3,1,3)
plot(t,split)
axis([0 24 0 1]);
title('Split fraction')
xlabel('Time [h]')
ylabel('\alpha')

end

%% Case four

% change in mass_water and mass_ice when split and Q_chiller is turned on.
% split = 0.5 at time = 10 h
% split = -0.5 at time = 15 h
% Q_chiller = 1.5 at time = 20 h

if plotte == 7 
    
figure(7);

subplot(4,1,1);
plot(t,Mass_ice)
axis([0 24 0 60]);
set(gca,'YTick',[10 20 30 40 50 60])
title('Mass ice')
xlabel('Time [h]')
ylabel('Mass ice [kg]')

subplot(4,1,2);
plot(t,Mass_water)
axis([0 24 70 150]);
set(gca,'YTick',[80 100 120 140])
title('Mass water')
xlabel('Time [h]')
ylabel('Mass water [kg]')

subplot(4,1,3)
plot(t,split)
axis([0 24 0 1]);
title('Split fraction')
xlabel('Time [h]')
ylabel('\alpha')

subplot(4,1,4);
plot(t,Q_chiller)
axis([0 24 -1.5 3]);
title('Q chiller')
xlabel('Time [h]')
ylabel('Q chiller [kW]')

end

% %% lagre
% 
% y = 1;
% res_Troom{y} = Troom;
% res_Tradiator{y} = Tradiator;
% res_split{y} = split;
% res_Toutdoor{y} = T_outdoor;
% res_time{y}= t;
% res_ice{y} = Mass_ice;
% res_water{y} = Mass_water;
% 
%   
% for y = 1:3
% 
% subplot(4,1,1);
% plot(res_time{y},res_Troom{y})
% axis([0 24 10 32]);
% title('Temp. room')
% %xlabel('Time [h]')
% ylabel('Temp. room [\circC]')
% 
% hold on
% 
% subplot(4,1,2);
% plot(res_time{y},res_Tradiator{y})
% axis([0 24 0 32]);
% title('Temp. radiator')
% %xlabel('Time [h]')
% ylabel('Temp. radiator [\circC]')
% 
% hold on
% 
% subplot(4,1,3);
% plot(res_time{y},res_ice{y})
% axis([0 24 15 60]);
% title('Mass ice')
% %xlabel('Time [h]')
% ylabel('Mass ice [kg]')
% 
% hold on
% 
% subplot(4,1,4);
% plot(res_time{y},res_water{y})
% axis([0 24 70 125]);
% title('Mass water')
% xlabel('Time [h]')
% ylabel('Mass water [kg]')
% 
% hold on
% 
% end
 