%% 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  =   10;                 %wysokość impulsu
szerimp =   6;                        %szerokość impulsu
Tzew1	=	Tzew0+wysimp;           %temp zewnętrzna
f1      =   f0;                     %przepływ wody w grzejniku

%% Symulacja
simtime =   5500;
simstep =   1;
t1      =   50;                    %czas zakłócenia
t2      =   simtime;                %szerokość impulsu

model	=	'sim1_model';
trans	=	'sim2_model';

sim(model);

clf;
hold on;

Tw_max  = max(Tw);              %stan ustalony
t_max   = 500;                      %granica kreślenia wykresu

%plot(t,Inpu, 'b');                 %wyświetla pobudzenie
grid on;
TW_max = ones(t_max,1) * Tw_max;
plot(t(1:t_max),Tw(1:t_max), 'g');     %wyświetla odp. imp.
plot(t(1:t_max),TW_max, 'b');      %poziom stanu ustalonego

a_min = t1;

% oblicza punkt przecięcia się stycznej i poziomu stanu ustalonego
%for a=1:150;
%   if ( (Tw_max - t(a+t1)*((Tw(t1+2) - Tw(t1+1)))) < Tw_max - t(a_min)*(Tw(t1+2) - Tw(t1+1)) )
%       k = (Tw(t1+2) - Tw(t1+1)
%       z = Tw_max - t(a_min)*(Tw(t1+2) - Tw(t1+1))
%       a_min = a+t1
%   end
%end

%plot(t(a_min),t(a_min-t1)*(Tw(t1+2) - Tw(t1+1))+Tw0,'or');

tt = t(1:150);
plot(tt+t1,tt*(Tw(t1+2) - Tw(t1+1))+Tw0,'-k');

numm = (Tw_max - Tw0)/(Tzew1 - Tzew0)
denm = [180 1]
sim(trans);
plot(t(1:t_max),stepik(1:t_max)/42+Tw0, 'r');
hold off;
