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

    cpp     =	1000;                   %cieplo wlasciwe powietrza
    ep		= 	1.2;                    %gestosc pow
    Vp		=	3 * 2.5 * 3;            %objetosc pomieszczenia
    Cvg     =	cpp * ep * Vp;          %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

    Qd0     =   0;                      %ciepła dodatkowe
    f0      =   fN;                     %przepływ wody w grzejniku
    Tgz0    =   0;                      %temp grzejnika
    Tzew0	=	0;                      %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łócenia
    Qd1     =   Qd0;
    t1      =   2000;                   %czas zakłócenia
    Tgz1    =   Tgz0;                   %temp grzejnika
    Tzew1	=	Tzew0+10;               %temp zewnętrzna
    f1      =   f0;                     %przepływ wody w grzejniku

    % Zmienne pomocnicze
    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];

    [A1,B1,C1,D1] = tf2ss(num1,den);
    [A2,B2,C2,D2] = tf2ss(num2,den);
    [A3,B3,C3,D3] = tf2ss(num3,den);

    % Symulacja
    simtime =   11000;

    model	=	'sim1_model';

    colors = ['r', 'g', 'b'];

    clf;

    hold on;
    sim(model);

    plot(t,Tw, colors(1));
    plot(t,TwTF, colors(2));
    plot(t,TwSS, colors(3));

    hold off;

    xlabel('czas [s]');
    ylabel('temperatura [C]');
    title('Wplyw Temperatury zewnętrznej na temperaturę pomieszczenia');
    legend('Model Rzeczywisty','Model Operatorowy','Model Przestrzeni Stanów');