%PROJEKT 4.1
%Wykres Molliera-Ramzina
%dla uk³adu powietrze-para wodna

function Mollier
clear all; close all; clc;
global P0 fi Y 
% warning off all
sY=1;
si=1;
hv0=hvp(0);
%-------------------------------------------------
P0=101.325;     % ciœnienie otoczenia [kPa]
fip=0.1;        % pocz¹tkowa wilgotnoœæ wzglêdna
dfi=0.1;        % krok wilgotnoœci wzglêdnej
Yp=1e-9;        % pocz. wilg. bezwzglêdna [kg/kg]
nY=1000;        % liczba kroków wilg. bezwzglêdnej [kg/kg]
Yk=.1;          % koñc. wilg. bezwzglêdna [kg/kg]
Tp=10;          % pocz¹tkowa izoterma [C]
dT=10;          % krok izotermy [C]
Tk=200;         % koñcowa izoterma[C]
ip=40;          % pocz¹tkowa entlapia [kJ/kg]
di=20;          % krok entalpii [kJ/kg]
ik=500;         % koñcowa entalpia [kJ/kg]
opis=1;         % w³¹cz opis =1 wy³¹cz opis = 0
fs=8;           % rozmiar czcionki opisu wykresu
%-------------------------------------------------
% wektor wilgotnoœci wzglêdnej
fz=(fip:dfi:1)';
% wektor wilgotnoœci bezwzglêdnej
Yz=linspace(Yp,Yk,nY);
% obliczam temperaturê wrzenia cieczy pod zadanym ciœnieniem
Tb=fzero(@psTb,100);
g1=(cp(Tb).*Tb*(sY./si).*(Yz./(si./cg(Tb)))+cg(Tb).*Tb./si)';
% obliczam linie wilgotnoœci wzglêdnej
for i=1:length(fz)
    fi=fz(i);
    t0=-80;
    for j=1:length(Yz)
        Y=Yz(j);
        yt(j,i)=fzero(@cis,t0);
        t0=floor(yt(j,i));
        yfi(j,i)=cp(yt(j,i)).*yt(j,i).*(sY./si).* (Y./sY)+yt(j,i)./(si./cg(yt(j,i)));
    end
end
%----------Budowa wykresu---------------------------------
figure(1)
% procedura rysuj¹ca linie wilgotnoœci wzglêdnej
[~,kg]=size(yfi);
for i=1:kg
    ypg=yfi(:,i);
    nrg=find(ypg <= g1);
    plot(Yz(nrg),ypg(nrg),'g')
    hold on
    [w1,~]=size(nrg);
    str = num2str(fz(i));
    if opis == 1 
text(Yz(nrg(w1)),yfi(nrg(w1),i),str, 'FontSize',fs, ...
        'HorizontalAlignment','center', 'EdgeColor','g',  'BackgroundColor',[1 1 1]);
    end
end
[w1,k1]=size(yfi);
% obliczam izotermy dla powietrza wilgotnego
Tz=(Tp:dT:Tk)';
for i=1:length(Tz)
    T=Tz(i);
    for j=1:length(Yz)
        Y=Yz(j);
        yt(j,i)=cp(T).*T*(sY./si).*(Y./(si./cg(T)))+ cg(T).*T./si;
    end
end
% procedura rysuj¹ca i opisuj¹ca izotermy
[~,k2]=size(yt);
for i=1:k2
    ypz=yt(:,i);
    nr=find(ypz >= yfi(:,k1));
    [wn,~]=size(nr);
    plot(Yz(nr),ypz(nr),'r','LineWidth',1)
    str = num2str(Tz(i));    
    if opis == 1
        yy=ypz(nr(round(0.5*wn))); %jeœli pionowa wspó³rzêdna tagu
        if yy<yfi(w1,1)            %nie wychodzi za ramkê
        text(Yz(1,nr(round(0.5*wn))),yy,str, 'FontSize',fs,...
      'HorizontalAlignment','center', 'EdgeColor','r','BackgroundColor',[1 1 1]);
        end
    end
end
plot(Yz,cp(Tb).*Tb*(sY./si).*(Yz./(si./cg(Tb)))+cg(Tb).*Tb./si,'r-.')
hold on
% obliczam entalpie dla powietrza wilgotnego
iz=(ip:di:ik)';
for k=1:length(iz)
    i=iz(k);
    for j=1:length(Yz)
        Y=Yz(j);
        yi(j,i)=-hv0.*(sY./si).*Y+i./si;
    end
end
%procedura rysuj¹ca izoentalpy
[~,k3]=size(yi);
for i=1:k3
    yiz=yi(:,i);
    nr=find(yiz >= yfi(:,k1));
    plot(Yz(nr),yiz(nr),'b')
end
axis([0 Yk 0 yfi(w1,1)])
pres=[num2str(P0) ' kPa'];
title(['Wykres Molliera. Uk³. woda-powietrze. P=' pres])
xlabel('Y [kg/kg]'); ylabel('i [kJ/kg]')
text(0.85*Yk,20,'t=const','EdgeColor','r')
text(0.85*Yk,10,'\phi=const','EdgeColor','g');
set(gca,'XGrid','on','fontsize',6)

%Funkcje opisuj¹ce w³aœciwoœci
%Aby wykonaæ wykres dla innego uk³adu gaz-ciecz nale¿y
%jedynie wymieniæ te funkcje.
%-------------------------------------------------
function y=cis(T)
% przeliczenie wilgotnoœci wzglêdnej na bezwzglêdn¹
global P0 fi Y
M=18.02; Mg=28.9645;
y=((M./Mg).*((fi*ps(T))./(P0-fi.*ps(T))))-Y;
%-------------------------------------------------
function y=psTb(T)
% funkcja obliczaj¹ca temperaturê wrzenia w [C]
global P0
y=ps(T)-P0;
%-------------------------------------------------
function y=ps(T)
% prê¿noœæ pary [kPa] od temperatury [C]
a1=16.376953; a2=3878.8223; a3=229.862;  
y=exp(a1-a2/(a3+(T)));
%-------------------------------------------------
function y=hvp(T)
% ciep³o parowania [kJ/kg] od temperatury [C]
a1=52.053; a2=647.13; a3=0.321;  
y=(a1.*(1-(T+273.15)./a2).^a3).*(1000./18.015);
%-------------------------------------------------
function y=cp(T)
% ciep³o w³aœciwe pary [kJ/(kg*K)] od temperatury [C]
a1=33.933; a2=-8.4186e-3; a3=2.9906e-5; a4=-1.7825e-8; a5=3.6934e-12;	
y=(a1+a2.*(T+273.15)+a3.*(T+273.15).^2+a4.* ...
  (T+273.15).^3+a5.*(T+273.15).^4).*(1./18.015);
%-------------------------------------------------
function y=cg(T)
% ciep³o w³aœciwe powietrza [kJ/(kg*K)] od temperatury [C]
a1=1.02287;	a2=-0.0512e-3; a3=0.181871e-6; a4=-0.05122e-9;
y=(a1+a2.*(T+273.15)+a3.*(T+273.15).^2+a4.*(T+273.15).^3);
%(C)-2014 Pakowski Z. Adamski R.
%"Podstawy MATLABa w in¿ynierii procesowej",Wyd.P£ £ódŸ 
%ISBN 978-83-7283-630-4
