%% 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


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

sim(model);

clf;
hold on;

%plot(t,Inpu, 'r');                 %wyświetla pobudzenie
plot(t,Tw-Tw0, 'g');                %wyświetla odp. imp.

%% normalizacja do całki z pobudzenia

m0 = m0./(wysimp*szerimp);
m1 = m1./(wysimp*szerimp);
m2 = m2./(wysimp*szerimp);
m3 = m3./(wysimp*szerimp);
m4 = m4./(wysimp*szerimp);
m5 = m5./(wysimp*szerimp);

%% pierwsza macierz rownania dla transmitancji: st.l. 2, st.m. 3
Ml = [            0                0                 0   1  0 0;
      m0(simtime+1)                0                 0   0 -1 0;
      m1(simtime+1)    -m0(simtime+1)                0   0  0 1; 
      m2(simtime+1)/2  -m1(simtime+1)    m0(simtime+1)   0  0 0;
      m3(simtime+1)/6  -m2(simtime+1)/2  m1(simtime+1)   0  0 0;
      m4(simtime+1)/24 -m3(simtime+1)/6  m2(simtime+1)/2 0  0 0];
  
%% druga macierz rownania dla transmitancji: st.l. 2, st.m. 3  
Mp = [ m0(simtime+1);
       m1(simtime+1);
       m2(simtime+1)/2;
       m3(simtime+1)/6;
       m4(simtime+1)/24;
       m5(simtime+1)/120];

%% rozwiązanie układu rownań Ml * A = Mp (pinv ponieważ Ml jest osobliwa)
A = pinv(Ml) * Mp;

numm = [A(6) A(5) A(4)];           % [b2 b1 b0]
denm = [A(3) A(2) A(1) 1];          % [a3 a2 a1 1]
TF = tf(numm, denm);

fprintf('%f %f %f\n', A(6), A(5), A(4));
fprintf('%f %f %f %f\n', A(3), A(2), A(1), 1);

% impulsowa z CST
%impulse(TF);

%% porownanie modeli
Tzew1	=	Tzew0+20;               %temp zewnętrzna

t1      =   100;                    %czas zakłócenia
t2      =   simtime;                %szerokość impulsu

sim(model);
sim(trans);


plot(t,Tw, 'r');                    %odpowiedź skokowa modelu rzeczywistego
plot(t,stepik+Tw0, 'b');            %odpowiedź skokowa modelu operatorowego

legend('impuls - momenty','skok - rzeczywisty','skok - transmitancja');

hold off;
