% Parametry obiektu
PgN     =  5000;                    %moc nominalna pieca
TwN     =    20;                    %temp wewn
TgzN	=    90;                    % TgzN/TgpN = 90/70
TgpN	=    70;
TzewN   =   -20;                    %temp zewn

cpw     =	4175;                   %cieplo wlasciwe wody
ew		= 	1000;                   %gestosc wody
Vw		=	0.1 * 1 * 1;            %objetosc grzejnika
Cvw     =	cpw * ew * Vw;          %pojemnosc cieplna

cpp     =	1000;                   %cieplo wlasciwe powietrza
ep		= 	1.2;                    %gestosc pow
Vp		=	3 * 2.5 * 3;            %objetosc pomieszczenia
Cvg     =	cpp * ep * Vp;          %pojemnosc cieplna

TgsrN	=	(TgzN + TgpN)/2;        %śrwdnia nominalna temp grzejnika
kg		=	PgN / (TgsrN - TwN);	%wsp wymiany - grzejnik pow
kstr	=	PgN / (TwN - TzewN);	%wsp strat
fN      =   PgN / (2 * cpw * ew * (TgzN - TgsrN));

% Warunki Początkowe
Qd0     =   0;                      %ciepła dodatkowe
f0      =   fN;                     %przepływ wody w grzejniku
Tgz0    =   TgzN;                   %temp grzejnika            
Tzew0	=	TzewN;                  %temp zewnętrzna

Qd0     =   0;                      %ciepła dodatkowe
f0      =   fN;                     %przepływ wody w grzejniku
Tgz0    =   TgzN;                   %temp grzejnika            
Tzew0	=	-30;                  %temp zewnętrzna

%średnia temp grzejnika
Tgsr0   = (2 * cpw * ew * fN * kstr * Tgz0 + kg * (Qd0 + 2 * cpw * ew * fN * Tgz0 + kstr * Tzew0))/(kg * kstr + 2 * cpw * ew * fN * (kg + kstr));
%średnia temp wnętrza
Tw0     = (Qd0 + kg * Tgsr0 + kstr * Tzew0) / (kg + kstr);

% Zakłócenia
Qd1     =   Qd0;  
t1      =   2000;                   %czas zakłucenia
Tgz1    =   Tgz0;                   %temp grzejnika
Tzew1	=	Tzew0;                  %temp zewnętrzna
f1      =   f0;                     %przepływ wody w grzejniku

% Symulacja
model	=	'sim1_model';

params  = [0, 1000, 2000, 3000, 4000, 5000, 6000];
colors = ['r', 'g', 'b', 'c', 'm', 'y', 'k'];

clf;

subplot(2,1,1);

hold on;
grid on;
for a = 1:size(params,2);
    Qd1 = params(a);
    sim(model);
    plot(t,Tw, colors(a));    
end
hold off;

xlabel('czas [s]');
ylabel('temperatura [C]');
title('Wpływ Q_d na temperaturę pomieszczenia');
legend('Q_d=0','Q_d=1000','Q_d=2000','Q_d=3000','Q_d=4000','Q_d=5000','Q_d=6000');

subplot(2,1,2);

hold on;
grid on;
for a = 1:size(params,2);
    Qd1 = params(a);
    sim(model);
    plot(t,Tgsr, colors(a));    
end
hold off;

xlabel('czas [s]');
ylabel('temperatura [C]');
title('Wpływ Q_d na temperaturę średnią grzejnika');
legend('Q_d=0','Q_d=1000','Q_d=2000','Q_d=3000','Q_d=4000','Q_d=5000','Q_d=6000');
