%skrypt z danymi do model3.mdl - symulacja pomieszczenia z kaloryferem
model='model4';
clear all;

%************Wielko��i wej�ciowe*****************************************

TzewN=-20;  % Temperatura zewnetrzna nominalna
TwewN=20;   % Temperatura wewnetrzna nominalna
TgzN=90;    % Temperatura grzejnika zasilania nominalna
TgpN=70;    % Temperatura grzejnika na powrocie nominalna
QN=5000;    % Moc elementu grzewczego
QdN=0;      % Moc dodatkowa
Cpp=1000;   % Cieplo wlasciwe powietrza
Cpw=4175;   % Cieplo wlasciwe wody
ROp=1.2;    % Gestosc powietrza
ROw=1000;   % Gestosc wody
Vw=32.4;    % objetosc pomieszczenia w m3
Vg=0.5;		% objetosc wody w m3
Step_time= 0000; % Czas opoznienia
max_step_size=10;
stop_time=20000;

%***********Wsp�czynniki i warunki poczatkowe**************************

Cvw=Cpp*ROp*Vw % 
Cvg=Cpw*ROw*Vg % 
fN=QN/(Cpw*ROw*(TgzN-TgpN));    % Nominalna funkcja przeplywu wody przez grzejnik 
kg=QN/(TgpN-TwewN);     % Wspolczynnik grzania
kstr=QN/(TwewN-TzewN);  % Wspolczynnik strat
Tzew0=TzewN;    % Temperatura zewnetrzna poczatkowa
Tgz0=TgzN;      % Temperatura grzejnika zasilania poczatkowa
Qd0=QdN;        % Dodatkowe zrodla ciepla nominalne 
f0=fN;          % Funkcjia przeplywu wody przez grzejnik

                        % Temperatura wewnetrzna pomieszczenia
Tw0= ((kg*kstr*Tzew0 + kg*Qd0 + Cpw*ROw*f0*Tgz0*(kg + kstr)) / (Cpw*ROw*f0*(kg + kstr) + kg*kstr)) * ((Cpw*ROw*f0)/kg + 1) - (Cpw*ROw*f0*Tgz0)/kg;
                        % Srednia temperatura w pomieszczeniu
Tgsr0= (kg*kstr*Tzew0 + kg*Qd0 + Cpw*ROw*f0*Tgz0*(kg + kstr)) / (Cpw*ROw*f0*(kg + kstr) + kg*kstr);


    
    dTzew0=5*500;       % Zmiana temperatury zewnetrznej poczatkowej
    dTgz0=0;        % Zmiana temperatury grzejnika zasilania             
    dQd0=0;    % Zmiany dodatkowych zrodel ciepla nominalnych
    df0=0;          % Zmiana funkcji przeplywu wody przez grzejnik
    
    sim('model4.mdl')
    %     sim('Kupf.mdl')
%     sim('moment.mdl')




% 
clf;
plot(t,Twew,'b');
    title('Twew(t)');
    xlabel('t [s]');
    ylabel('Twew [^oC]');
    grid on;
%     legend('Obiekt','Model');
%     
    
%     
%     clf;
% plot(t,Tgsr,'b',t_g,Twew_g,'g --');
%     title('Tgsr(t)');
%     xlabel('t [s]');
%     ylabel('Tgsr [^oC]');
%     grid on;
%     legend('Obiekt','Model');

%  plot(t,Tgsr);
%  grid on;

%  plot(t,Twew);
%  grid on;


%     hold on;
% plot(t_g,Twew_g,'g');    
% end
% drawnow; 

% x=size(Twew,1);
% Twu=Twew(x)'
% M0=0;
% M1=0;
% 
% for j=1:(x)
%     M0=M0+(Twu-Twew(j));
%     M1=M1+j*(Twu-Twew(j));
% end
% 
% M0
% M1
% Tw0=Tw0;
% 
% T11=(2*M1*(Twu-Tw0)-2*(M0)^2)/((Twu-Tw0)*(-M0-sqrt(-3*(M0)^2+4*M1*(Twu-Tw0))))
% T22=(2*M1*(Twu-Tw0)-2*(M0)^2)/((Twu-Tw0)*(-M0+sqrt(-3*(M0)^2+4*M1*(Twu-Tw0))))
% k=Twu-Tw0
% 
% aa=M0-M1/Twu;
% bb=M0-2*M0*Twu;
% cc=Twu*Twu-Twu+1;
% 
% delta=(-4*aa*cc)+bb*bb;
% 
% if delta~=0
% b1=(-bb-sqrt(delta))/(2*aa);
% b2=(-bb+sqrt(delta))/(2*aa);
%     if(b1>0)
%         b=b1;
%     else
%         b=b2
%     end
% 
% else
% b=bb/(2*aa);    
% end
% 
% 
% a=b/(M0*b-Twu);
% 
% T1=1/a
% T2=1/b
% a;
% b;
% 
% y=Twu;
% a1=(y - (M0*y)/(2*(M0 - y)) + y^2/(2*(M0 - y)) + sqrt(-3*M0^2 * y^2 + 4 * M0 * y^3)/(2*(M0 - y)) - (y*sqrt(-3* M0^2* y^2 + 4* M0* y^3))/(2* M0* (M0 - y)))/(M0 - y)
% b1=(M0* y - sqrt(-3* M0^2* y^2 + 4* M0* y^3))/(2* M0* (M0 - y))
% a2=(y - (M0* y)/(2 *(M0 - y)) + y^2/(2 *(M0 - y)) - sqrt(-3* M0^2* y^2 + 4* M0* y^3)/(2* (M0 - y)) + (y* sqrt(-3* M0^2* y^2 + 4* M0* y^3))/(2* M0* (M0 - y)))/(M0 - y)
% b2=(M0* y + sqrt(-3* M0^2* y^2 + 4* M0 *y^3))/(2* M0* (M0 - y))

