%Parametry ************************************************************
Qgn = 1000;
Twewn = 20;
Tzewn = -20;
Tpn = 10;
a = 0.25; %Kcpw = a * Kcw

Kcw = Qgn / ((Twewn - Tzewn) + a * (Twewn  - Tpn));
Kcpw = (Qgn - Kcw*(Twewn - Tzewn))/(Twewn - Tpn);
Kcp = (Kcpw*(Twewn - Tpn))/(Tpn - Tzewn);

%Zakres_zmiennych na charakterystykach *******************************
q_vector = linspace(0,2000,1000);
Tzew_vector = linspace(-40,40,1000);
Tzew_chosen_points = [-20, -10, 0,10,20];
q_chosen_points = [0,100,500,1000,1500,2000];

%Twew(Qgn) ************************************************************
xlabel('Moc grzejnika (W)');
ylabel('Temperatura wewnętrzna(Celc)');
title('Twew(Qg) przy zadanym Tzew');
grid;
hold on;
for i = 1:length(Tzew_chosen_points)
    Twew = Twew_values(q_vector,Tzew_chosen_points(i),Kcp,Kcpw,Kcw);
    leg = sprintf('Tzew = %d',Tzew_chosen_points(i));
    plot(q_vector,Twew,'DisplayName',leg);
end
legend;
Twew_n = Twew_values(Qgn,Tzewn,Kcp,Kcpw,Kcw);
scatter(Qgn,Twew_n,'DisplayName','Punkt Nominalny');

%Twew(Tzew) ************************************************************
figure();
xlabel('Temperatura zewnętrzna(Celc)');
ylabel('Temperatura wewnętrzna(Celc)');
title('Twew(Tzew) przy zadanym Qg');
grid;
hold on;
b = Kcp + Kcpw;
for i = 1:length(q_chosen_points)
    Twew = Twew_values(q_chosen_points(i),Tzew_vector,Kcp,Kcpw,Kcw);
    leg = sprintf('Qg = %d',q_chosen_points(i));
    plot(Tzew_vector,Twew,'DisplayName',leg);

end
legend;
Twew_n = Twew_values(Qgn,Tzewn,Kcp,Kcpw,Kcw);
scatter(Tzewn,Twew_n,'DisplayName','Punkt Nominalny');

%Tp(Qg) *************************************************************
figure();
xlabel('Moc grzejnika (W)');
ylabel('Temperatura poddasza (Celc)');
title('Tp(Qg) przy zadanym Tzew');
grid;
hold on;
for i = 1:length(Tzew_chosen_points)
    Tp = Tp_values(q_vector,Tzew_chosen_points(i),Kcp,Kcpw,Kcw);
    leg = sprintf('Tzew = %d',Tzew_chosen_points(i));
    plot(q_vector,Tp,'DisplayName',leg);
end
legend;
Tpn = Tp_values(Qgn,Tzewn,Kcp,Kcpw,Kcw);
scatter(Qgn,Tpn,'DisplayName','Punkt Nominalny');

%Tp(Qg) *************************************************************
figure();
xlabel('Temperatura zewnętrzna (Celc)');
ylabel('Temperatura poddasza (Celc)');
title('Tp(Tzew) przy zadanym Qg');
grid;
hold on;
for i = 1:length(q_chosen_points)
    Tp = Tp_values(q_chosen_points(i),Tzew_vector,Kcp,Kcpw,Kcw);
    leg = sprintf('Qg = %d',q_chosen_points(i));
    plot(Tzew_vector,Tp,'DisplayName',leg);
end
legend;
Tpn = Tp_values(Qgn,Tzewn,Kcp,Kcpw,Kcw);
scatter(Tzewn,Tpn,'DisplayName','Punkt Nominalny');

function Twew = Twew_values(q,Tzew,Kcp,Kcpw,Kcw)
b = Kcp + Kcpw;
    Twew = (q.*b + b.*Kcw.*Tzew + Kcp*Kcpw.*Tzew)/((Kcw+Kcpw).*b -Kcpw.^2);
end

function Tp = Tp_values(q,Tzew,Kcp,Kcwp,Kcw)
    Tp = (q - Tzew.*((Kcw/Kcwp) - Kcw + 1)) / (Kcw + Kcw*Kcp/Kcwp + Kcwp);
end