%I (obiekt)
TwN = 20;                 %temperatura wewnetrzna
TzewN = -20;              %temperatura zewnetrzna
TgzN = 90;                %TgzN/TgpN = 90/70
TgpN = 70;      
Cpp = 1000;               %cieplo wlasciwe powietrza
Cpw = 4175;               %cieplo wlasciwe wody
rop = 1.2;                %gestosc powietrza
row = 1000;               %gestosc wody
Vp = 3 * 3 * 3;           %objetosc pomieszczenia(powietrza)
Vw = 1 * 1 * 0.1;         %objetosc grzejnika(wody)
QN = 5000;                %cieplo oddawane przez piec
QdN = 0;

T0=100; %opóźnienie

cvg=Cpw*row*Vw;
cvw=Cpp*rop*Vp;

TgsrN = (TgzN + TgpN) / 2;
TgpN = 2 * TgsrN - TgzN;
FN = QN / (2 * Cpw * row * (TgzN - TgsrN));
kg = QN / (TgsrN - TwN);
kstr = (QN + QdN) / (TwN - TzewN);

%II (warunki początkowe)
Tzew0=TzewN;
Tgz0=TgzN;
F0=FN;
Qd0=QdN;

%III (zmiany)
dTzew0=10;
dTgz0=0;
dF0=0;
dQd0=0;

Tgsr0 = ( 2 * Cpw*row*F0*Tgz0*(kg + kstr) + kg*(kstr*Tzew0 + Qd0) ) / ( kg * kstr + 2 * Cpw * row * F0 * (kg + kstr) );
Tw0 = (Qd0 + kg * Tgsr0 + kstr * Tzew0) / (kg + kstr);

%Transmitancja pokoju w postaci L/M  TwN/TzewN
L=[(cvw * Cpw * FN * row) (kg * kstr + Cpw * FN * row * (kg + kstr))];
M=[(cvw * cvg) (kg * (cvw + cvg) + Cpw * FN * row * cvw + cvg * kstr) (kg * kstr + Cpw * FN * row * (kg + kstr)) ];

%Równania stanu dla pokoju
a0=(kg*kstr + Cpw*FN*row*(kg+kstr)) / (cvw*cvg);
a1=(kg*(cvw+cvg) + Cpw*FN*row*cvw+cvg*kstr)/ (cvw*cvg);

b0=(kg*kstr + Cpw*FN*row*(kg+kstr))/ (cvw*cvg);
b1=(cvw*Cpw*FN*row)/ (cvw*cvg);

A=[0 1; -a0 -a1];
B=[0; 1];
C=[b0 b1];
D=[0];

sys1=ss(A,B,C,D);
sys2=tf(L,M);

%--------------------------------------------------------------------------

%Tgsr/TzewN dla grzejnika
L1=[Cpw*FN*kg*row];
M1=[(cvw*cvg) (kg*(cvw+cvg)+Cpw*FN*row*cvw+cvg*kstr) (kg*kstr+Cpw*FN*row*(kg+kstr)) ];

%Równania stanu dla grzejnika
a1=(kg*(cvw+cvg)+Cpw*FN*row*cvw+cvg*kstr) / (cvw*cvg);
a0=(kg*kstr+Cpw*FN*row*(kg+kstr)) / (cvw*cvg);

b1=0;
b0=(Cpw*FN*kg*row) / (cvw*cvg);

A1=[0 1; -a0 -a1];
B1=[0; 1];
C1=[b0 b1];
D1=[0];

sys3=ss(A1,B1,C1,D1);
sys4=tf(L1,M1);
