clc, clear all

% Vstupní hodnoty

p=[0,0,0,0]; 
s=[0,0,0,0];
T=[0,0,0,0]; 
h=[0,0,0,0]; 
v=[0,0,0,0];
p(2)=input('Tlak v parogenerátoru (Pressure in a steam generator) [MPa]:');
T(3)=input('Teplota přehřáté páry na vstupu do turbíny (Temperature at enter in the turbine) [°C]:');  
p(4)=input('Tlak v kondenzátoru (Pressure in a condenser) [kPa]:');

%
% Výpočet stavových veličin

p(4)=p(4)/1000; % převod do stejných jednotek (kPa -> MPa) 
p=10*p; % korekce jednotek do XSteam (MPa -> bar) 

% využití XSteam pro určení potřebných energetických veličin 
p(1)=p(4); % izobarická kondenzace
s(1)=XSteam('sL_p',p(1));
T(1)=XSteam('Tsat_p',p(1)); 
h(1)=XSteam('hL_p',p(1)); 

s(2)=s(1); % isoentropická komprese
T(2)=XSteam('T_ps',p(2),s(2)); 
h(2)=XSteam('h_ps',p(2),s(2)); 

p(3)=p(2); % izobarický ohřev 
s(3)=XSteam('s_pT',p(3),T(3));
h(3)=XSteam('h_pT',p(3),T(3));

s(4)=s(3); % isotermická expanze
T(4)=XSteam('T_ps',p(4),s(4));
h(4)=XSteam('h_ps',p(4),s(4));

%
% Termodynamická analýza

% měrná práce čerpadla 
aN=(h(2)-h(1)); 
% dodané teplo v parogenerátoru
qH=h(3)-h(2); 
% měrná práce turbíny 
aT=h(3)-h(4); 
% odvedené teplo v kondenzátoru
qC=h(4)-h(1);

% termická účinnost cyklu 
nt=(1-((h(4)-h(1))/(h(3)-h(2))))*100;

%
% Výstupní hodnoty

fprintf('Termická účinnost cyklu: %.2f %%\n',nt);

fprintf('Měrná práce čerpadla: %.2f kJ/kg\n',aN); 
fprintf('Měrné přivedené teplo: %.2f kJ/kg\n',qH);
fprintf('Měrná práce turbíny: %.2f kJ/kg\n',aT);
fprintf('Měrné odvedené teplo: %.2f kJ/kg\n',qC); 

%
% Vykreslení cyklu v T-s diagramu

% vytvoření okna grafu, název a popisy os
figure; hold on;
if  p(2)>=220.6
    title('T-s diagram nadkritického Rankinova cyklu') 
else
    title('T-s diagram základního Rankinova cyklu')
end
xlabel('s [kJ/kg*K]') 
ylabel('T [^\circC]')

% vykreslení mezních křivek 
Tmz = linspace(0,400,400); 
shorni = arrayfun(@(x) XSteam('sV_T',x),Tmz); 
sdolni = arrayfun(@(x) XSteam('sL_T',x),Tmz); 
plot(shorni,Tmz,'k-')
plot(sdolni,Tmz,'k-') 

% vykreslení jednotlivých dějů 
% děj 1-2 izoentropický
plot([s(1) s(2)],[T(1) T(2)],'r-') 
% děj 2-3 izobarický
T23=linspace(T(2),T(3),600); 
s23 = arrayfun(@(x) XSteam('s_PT',p(2),x),T23); 
plot(s23,T23,'r-') 
% děj 3-4 izoentropický
plot([s(3) s(4)],[T(3) T(4)],'r-') 
% děj 4-1 izobarický
s41 = linspace(s(4),s(1)); 
t41 = arrayfun(@(x) XSteam('T_ps',p(1),x),s41); 
plot(s41,t41,'r-')

% popis stavů
plot(s,T,'o','MarkerSize',3); 
text(s(1)-0.1,T(1)-15,'1'); 
text(s(2)-0.1,T(2)+20,'2'); 
text(s(3)+0.1,T(3),'3'); 
text(s(4),T(4)-10,'4'); 