%% Clean up
clc
clear all
close all

%% Constants
CelsiusToKelvin = 273.15;

%% Given (in bar and C)
p1 = 100;
p2 = 7;
p4 = 0.06;
T1 = 480;
T3 = 480;

p3 = p2;
p5 = p4;
p6 = p1;

eta12 = 0.8;
eta34 = 0.8;
eta56 = 0.8;

%% Point 1
% This is already done to give you an idea of how to proceed.
h1 = XSteam('h_pt',p1,T1);
s1 = XSteam('s_pt',p1,T1);

%% Process 1-2
h2is = XSteam('h_ps',##,##);
h2 = h1 + (h2is-h1)*eta12; 
s2 = XSteam('s_ph',p2,h2);
T2 = XSteam(##,##,##);

%% Point 3
h3 = XSteam(##,##,##);
s3 = XSteam(##,##,##);

%% Process 3-4
h4is = XSteam(##,##,##);
h4 = ##
s4 = XSteam(##,##,##);
T4 = XSteam(##,##,##);

%% Point 5
h5 = XSteam('hL_p',p5);
s5 = XSteam(##,##);
T5 = XSteam(##,##);
v5 = XSteam(##,##);

%% Process 5-6
h6is = h5 + v5*(p6-p5)*10^5/1000;
h6 = ##;

%% Efficiency
eta_th1 = ((h1-h2)+(h3-h4))/((h1-h6)+(h3-h2)-(h5-h6))
eta_th2 = ((h1-h2)+(h3-h4)+(h5-h6))/((h1-h6)+(h3-h2))

%% YOU DO NOT NEED TO CHANGE ANYTHING IN THIS SECTION
% This simply draws the Mollier Diagram based on the computed solution.
T6 = XSteam('T_ph',p6,h6);
s6 = s5 + 4.187*log((T6+CelsiusToKelvin)/(T5+CelsiusToKelvin));

TMin = 0;
TMax = 2000;
nIt = 4000;
for it = 1:nIt
    T = TMin+it/nIt*(TMax-TMin);
    if ( isnan(XSteam('hL_T',T)) == 0 )
      hL(it) = XSteam('hL_T',T);
      sL(it) = XSteam('sL_T',T);
      hV(it) = XSteam('hV_T',T);
      sV(it) = XSteam('sV_T',T);
    end
end

plot([sL flip(sV)],[hL flip(hV)], 'k')
hold on

sMin = 0;
sMax = 10;
nIt = 4000;
ds = (sMax-sMin)/nIt;
for it = 1:nIt
    s(it) = sMin+it/nIt*(sMax-sMin);
    hiso1(it) = XSteam('h_ps',p2,s(it));
end

plot(s,hiso1,'-.b')
hold on

sMin = 0;
sMax = 10;
nIt = 4000;
ds = (sMax-sMin)/nIt;
for it = 1:nIt
    s(it) = sMin+it/nIt*(sMax-sMin);
    hiso2(it) = XSteam('h_ps',p1,s(it));
end

plot(s,hiso2,'-.b')
hold on

sMin = 0;
sMax = 10;
nIt = 4000;
ds = (sMax-sMin)/nIt;
for it = 1:nIt
    s(it) = sMin+it/nIt*(sMax-sMin);
    hiso3(it) = XSteam('h_ps',p5,s(it));
end

plot(s,hiso3,'-.b')
hold on

plot(s1,h1,'or')
hold on
plot(s2,h2,'or')
hold on
plot(s3,h3,'or')
hold on
plot(s4,h4,'or')
hold on
plot(s5,h5,'or')
hold on
plot(s6,h6,'or')

axis([0 9 0 3700])
xlabel('entropy [kJ/kg/K]')
ylabel('enthalpy [kJ/kg]')
legend('x=0/1','isobars','Location','NorthWest')