% Kompressor.m
% Entalpidifferanse i kompressor, ved gitt innløpstilstand (T, p) og gitt
% utløpstrykk, isentropisk kompresjon

clear all
p01 = 6.322e5;          % [Pa] Gitt
T01 = 273.15-33.19;      % [K] Gitt
V1 = 7.7e-4;            % [m3] Startverdi, gassvolum
p02 = 13.42e5;          % [Pa] P gitt
n = [1];

% PENG-ROBINSON
% FØR KOMPRESSOR:
v1 = V1;
dv1 = v1;
% Newtoniterasjoner til V ved ønsket trykk finnes:
while abs(dv1)>1e-6*v1
    a1 = co2_pr([T01,v1,n]');
    dv1 = (p01 + a1.g(2))/(-a1.H(2,2));
    v1 = v1 + dv1;
end
h1 = -a1.g(1)*T01 + a1.g(3)*n;   % Entalpi før kompressor
disp(['PENG ROBINSON'])
disp(['Entalpi før kompressor:   H1 = ',num2str(h1),' J/mol'])

% ETTER KOMPRESSOR: 
s02 = a1.g(1);          % Isentropisk kompresjon
g02 = [-s02 p02]';

t2 = T01;               % [K] Startverdi
v2 = v1;                % [m3]    "
x  = [t2 v2];
dx = x;
while abs(dx./x)>1e-6
    a2 = co2_pr([t2,v2,n]');
    dx = ((g02 + a2.g(1:2))'/(-a2.H(1:2,1:2)));
    x  = x + dx;
    t2 = x(1);
    v2 = x(2);
end
h2 = -a2.g(1)*t2 + a2.g(3)*n;   % Entalpi etter kompressor
disp(['Entalpi etter kompressor: H2 = ',num2str(h2),' J/mol'])

% Differanse:
Dh = h2 - h1;
disp(['Differanse:            H2-H1 =    ',num2str(Dh),' J/mol'])
t2
disp(' ')

% SPAN-WAGNER
% FØR KOMPRESSOR:
dV1 = V1;
while abs(dV1)>1e-6*V1
    A1 = co2_sw([T01,V1,n]');
    dV1 = (p01 + A1.g(2))/(-A1.H(2,2));
    V1 = V1 + dV1;
end
H1 = -A1.g(1)*T01 + A1.g(3)*n;   % Entalpi før kompressor
disp(['SPAN WAGNER'])
disp(['Entalpi før kompressor:   H1 = ',num2str(H1),' J/mol'])

% ETTER KOMPRESSOR: 
S02 = A1.g(1);          % Isentropisk kompresjon
G02 = [-S02 p02]';

T2 = T01;               % [K] Startverdi
V2 = V1;                % [m3]    "
X  = [T2 V2];
dX = X;
while abs(dX./X)>1e-6
    A2 = co2_sw([T2,V2,n]');
    dX = ((G02 + A2.g(1:2))'/(-A2.H(1:2,1:2)));
    X  = X + dX;
    T2 = X(1);
    V2 = X(2);
end
H2 = -A2.g(1)*T2 + A2.g(3)*n;   % Entalpi etter kompressor
disp(['Entalpi etter kompressor: H2 = ',num2str(H2),' J/mol'])

% Differanse:
DH = H2 - H1;
disp(['Differanse:            H2-H1 =    ',num2str(DH),' J/mol'])
T2
disp(' ')

