close all;
clear all;

%Przeplyw
cpp = 1008; %J/kg K, powietrze
rop = 1.185; %kg/m3, powietrze
cpw = 4175.4; %J/kg K, woda
row = 997.04; %kg/m3,  woda

%Wartosci nominalne
Twewn=20;
Tzewn=-20;
Tgzn1=90;
Tgzn2=90;
Tgpn=70;
Tkzn=90;
Tkpn=70;
qgn=5000;
qgn1=5000;
qkn=5000;
Top=7000;
fmgn1=(qgn/(cpw*(Tgzn1-Tgpn)));
fmgn2=(qgn1/(cpw*(Tgzn2-Tgpn)));
fmkn=(qkn/(cpw*(Tkzn-Tkpn)));
simTime=1000000;

%Wartosci poczatkowe
%t=2000;
Tzew0=Tzewn;
Tgz1=Tgzn1;
Tgz2=Tgzn2;
fmg1=fmgn1;
fmg2=fmgn2;
fmk0=fmkn;
Tkz0=Tkzn;
Tkp0=Tkpn;
qg=qgn;
qk0=qkn;
qt1=0;
qt2=0;
Top_wek=[];

%Skoki
t0=2000;
dTzew=0;
dTgz=0;
dfmg1=0;
dfmg2=0;
dqg=0;
dqt=0;
dfmk=0;
dqk=0;
dTkp=0;

%Grzejnik
hg=0.4;
wg=0.2;
lg=1.5;

%Pomieszczenie
h1=2.5;
w1=5;
l1=4;

h2=2.5;
w2=3;
l2=4;

%Kocioł
hk=0.4;
wk=0.2;
lk=1.5;

%Pojemności
Vp1=h1*w1*l1;
Vp2=h2*w2*l2;
Vg1=hg*wg*lg;
Vgw=hk*wk*lk;

%Pojemnosci cieplne
Cvg = cpw*row*Vg1;
Cvw1 = cpp*rop*Vp1;
Cvw2 = cpp*rop*Vp2;
Cvk = cpw*row*Vgw;

%Identyfikacja parametrow
Kcw=qgn/(Twewn-Tzewn);
Kcg=qgn/(Tgpn-Twewn);

%Stany równowagi
Twew1=(cpw*fmg1*Kcg*Tgz1+Kcw*(Kcg+cpw*fmg1)*Tzew0)/(Kcg*Kcw+cpw*fmg1*(Kcg+Kcw));
Tgp1=((Kcg+Kcw)*Twew1-Kcw*Tzew0)/Kcg;
Twew2=(cpw*fmg2*Kcg*Tgz2+Kcw*(Kcg+cpw*fmg2)*Tzew0)/(Kcg*Kcw+cpw*fmg2*(Kcg+Kcw));
Tgp2=((Kcg+Kcw)*Twew2-Kcw*Tzew0)/Kcg;

%Symulacje
%{
sim('pomieszczenie');
subplot(5,1,1);
hold on, grid on, title('Zmiana temperatury pomieszczenia');
xlabel('czas');
ylabel('temperatura');
plot(ans.t,ans.Twew);

subplot(5,1,2);
hold on, grid on, title('Zmiana temperatury wody w grzejniku');
xlabel('czas');
ylabel('temperatura');
plot(ans.t,ans.Tgp);

subplot(5,1,3);
hold on, grid on, title('Zmiana temperatury w kotle');
xlabel('czas');
ylabel('temperatura');
plot(ans.t,ans.Tkz);

subplot(5,1,4);
hold on, grid on, title('Zmiana temperatury pomieszczenia');
xlabel('czas');
ylabel('temperatura');
plot(ans.t,ans.Twew1);

subplot(5,1,5);
hold on, grid on, title('Zmiana temperatury wody w grzejniku');
xlabel('czas');
ylabel('temperatura');
plot(ans.t,ans.Tgp1);
%}
%CHARAKTERYSTYKI STATYCZNE
Tzew_wek = [-20:40];
Twew1_wek = [];
Twew2_wek = [];
Tkp_wek = [];
fmg_wek1 = linspace(0, fmgn1, 61);
fmg_wek2 = linspace(0, fmgn2, 61);
for i=1:length(Tzew_wek) %przeplyw nominal
    Twew1_wek(i)=(cpw*fmg1*Kcg*Tgz1+Kcw*(Kcg+cpw*fmg1)*Tzew_wek(i))/(Kcg*Kcw+cpw*fmg1*(Kcg+Kcw));
    Twew2_wek(i)=(cpw*fmg2*Kcg*Tgz2+Kcw*(Kcg+cpw*fmg2)*Tzew_wek(i))/(Kcg*Kcw+cpw*fmg2*(Kcg+Kcw));
    Tgp1_temp=((Kcg+Kcw)*Twew1_wek(i)-Kcw*Tzew_wek(i))/Kcg;
    Tgp2_temp=((Kcg+Kcw)*Twew2_wek(i)-Kcw*Tzew_wek(i))/Kcg;
    Tkp_wek(i) = (Tgp1_temp * fmg1 + Tgp2_temp * fmg2) / (fmg1 + fmg2);
end

figure(2)
subplot(4, 1, 1)
grid on;
hold on;
plot(Tzew_wek, Twew1_wek,'-.r','LineWidth',2)
plot(Tzew_wek, Twew2_wek,'-b','LineWidth',1)
plot(Tzewn,Twew1,'*','LineWidth',5)
xlabel('Tzew[°C]');
ylabel('Twew[°C]');
legend('Twew1','Twew2','TwewN');

subplot(4, 1, 2)
grid on;
hold on;
plot(Tzew_wek, Tkp_wek,'LineWidth',2)
plot(Tzewn, Tkp0,'*','LineWidth',5)
xlabel('Tzew[°C]');
ylabel('Tkp[°C]');
legend('Tkp','TkpN');

%%%%%
Twew1_wek = [];
Twew2_wek = [];
Tkp_wek = [];

for i=1:length(fmg_wek1) %przeplyw nominal
    Twew1_wek(i) = (cpw*fmg_wek1(i)*Kcg*Tgz1+Kcw*(Kcg+cpw*fmg_wek1(i))*Tzew0)/(Kcg*Kcw+cpw*fmg_wek1(i)*(Kcg+Kcw));
    Twew2_wek(i) = (cpw*fmg_wek2(i)*Kcg*Tgz2+Kcw*(Kcg+cpw*fmg_wek2(i))*Tzew0)/(Kcg*Kcw+cpw*fmg_wek2(i)*(Kcg+Kcw));
    Tgp1_temp=((Kcg+Kcw)*Twew1_wek(i)-Kcw*Tzewn)/Kcg;
    Tgp2_temp=((Kcg+Kcw)*Twew2_wek(i)-Kcw*Tzewn)/Kcg;
    Tkp_wek(i) = (Tgp1_temp * fmg_wek1(i) + Tgp2_temp * fmg_wek2(i)) / (fmg_wek1(i)+fmg_wek2(i));
end

subplot(4, 1, 3)
grid on;
hold on;
plot(fmg_wek1, Twew1_wek,'-.r','LineWidth',2)
plot(fmg_wek2, Twew2_wek,'-b','LineWidth',1)
plot(fmgn1,Twew1,'*','LineWidth',5)
xlabel('fmg[kg/s]');
ylabel('Twew[°C]');
legend('Twew1','Twew2','TwewN');

subplot(4, 1, 4)
grid on;
hold on;
plot(fmg_wek1, Tkp_wek,'LineWidth',2)
plot(fmg_wek2, Tkp_wek,'LineWidth',2)
plot(fmgn1, Tkpn,'*','LineWidth',5)
xlabel('fmg[kg/s]');
ylabel('Tkp[°C]');
legend('Tkp','TkpN');

%SCHODKI
Top_wek = [3000,7000,21000];
dqk=0.1*qkn;

figure(3)
grid on;
hold on;

for i = 1:length(Top_wek)
   Top = Top_wek(i);
   sim('pomieszczenie');
   plot(ans.t, ans.Tkp,'LineWidth',2);
   hold on;
end
xlabel('t[s]');
ylabel('Tkp[°C]');
legend('Top=3000','Top=7000','Top=21000');

dqk=0;
Top=3000;

%REAKCJA NA SKOK
simTime=400000;
dqk=0.1*qkn;
sim('pomieszczenie');

figure(4)
subplot(3, 1, 1)
grid on;
hold on;

plot(ans.t, ans.Twew1,'LineWidth',2);
plot(ans.t, ans.Twew2,'LineWidth',2);
xlabel('t[s]');
ylabel('Twew[°C]');
legend('Twew1', 'Twew2');

%

dqk=0;
dTzew=1;
sim('pomieszczenie');

subplot(3, 1, 2)
grid on;
hold on;

plot(ans.t, ans.Twew1,'LineWidth',2);
plot(ans.t, ans.Twew2,'LineWidth',2);
xlabel('t[s]');
ylabel('Twew[°C]');
legend('Twew1', 'Twew2');

simTime=5000;
sim('pomieszczenie');

subplot(3, 1, 3)
grid on;
hold on;

plot(ans.t, ans.Twew1,'LineWidth',2);
plot(ans.t, ans.Twew2,'LineWidth',2);
xlabel('t[s]');
ylabel('Twew[°C]');
legend('Twew1', 'Twew2');