clear;

%% Stale

PgN = 15000;                % moc nominalna grzejnika
TzewN = -20;                % nominalna temp zewn
TwewN = 20;                 % temp wewnetrzna
TgzN = 40;                  % nominalna temp grzejnika
TsN = 15;                   % temp scian

Vw = 2 * 3 * 4;             % kubatura pomieszczenia
cp = 1000;                  % c. walsc. powietrza
ep = 1.2;                   % gestosc powietrza
Cv = cp * ep * Vw;          % poj. cieplna pomieszenia

Vs = 0.1 * Vw;              % objetosc scian
cps = 800;                  % c. wlasc. betonu
eps = 2300;                 % gestosc betonu
Cs = cps * eps * Vs;        % poj. cieplna scian

ks = PgN / (TsN - TzewN);
kw = PgN / (TwewN - TsN);
fN = PgN / (cp * ep * (TgzN - TwewN));

%% Warunki poczatkowe
f0    = fN*0.5;
Tgz0  = TgzN+1;
Tzew0 = TzewN;

M1 = (cp * ep * f0 + kw);
M2 = (ks + kw);

Twew0 = (cp * ep * f0 * (ks + kw) * Tgz0  + kw * ks * Tzew0)/ (M1 * M2 - kw^2);
%Twew0 = (M1 * cp * ep * f0 * Tgz0 + 
Ts0 = (kw * Twew0 + ks * Tzew0) / (ks + kw) ;

%% Zaklocenia
t1    = 30000;
f1    = f0;
Tgz1  = Tgz0+1;
Tzew1 = Tzew0;
Twew1 = Twew0;
Ts1   = Ts0;

%% Symulacja

simtime = 100000;
simmodel = 'modello';
colors = ['r', 'g', 'b', 'c', 'm', 'y', 'k'];

sim(simmodel);

cla;
hold on;
axis auto;
grid on;
plot(t,Twew,'g');
plot(t,Ts,'r');

legend('T_{wew}(t)','T_s(t)');

