%% Czyszczenie

clear

%% 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
Vg		=	0.1 * 1 * 1;            %objetosc grzejnika
Cvg     =	cpw * ew * Vg;          %pojemnosc cieplna

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

TgsrN	=	(TgzN + TgpN)/2;        %średnia 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

%ś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łocenia
Qd1     =   Qd0;  
Tgz1    =   Tgz0;                   %temp grzejnika
% impuls - 100-1000 * wart ustal
wysimp  =   100*Tw0;                %wysokość impulsu
szerimp =   6;                      %szerokość impulsu
Tzew1	=	Tzew0+wysimp;           %temp zewnętrzna
f1      =   f0;                     %przepływ wody w grzejniku

%% Symulacja
simtime =   1500;
simstep =   1;
t1      =   000;                    %czas zakłócenia
t2      =   t1 + szerimp * simstep; %szerokość impulsu


%% Transmitancja

num1 = 2*kg*cpw*ew*f0;
num2 = [(Cvg*kstr) 2*cpw*ew*f0*kstr+kg*kstr];
num3 = [Cvg 2*cpw*ew*f0+kg];
den = [(Cvw*Cvg) (2*Cvw*cpw*ew*f0 + kg*Cvw + kg*Cvg + kstr*Cvg) (2*cpw*ew*f0)*(kg+kstr) + kg*kstr];

TF1 = tf(num1, den)
TF2 = tf(num2, den);
TF3 = tf(num3, den);

cla;
[mag,phase,w] = bode(TF1);
bode(TF1);
grid on;
hold on;

[w(1) phase(:,:,1)]                 % dla przes fazowego 0st
[w(19) phase(:,:,19)]               % dla przes fazowego -45st

num_a = mag(:,:,1);                 % wzmoznienie
den_a = [1/w(1) 1/w(19) 1]

TF = tf(num_a, den_a);
bode(TF);

legend('Obiekt rzeczywisty','Model');