% co2_pr.m

function [S] = co2_pr(x)

[S] = Helmholtz_anonymous_23829684(x);


function [S] = Helmholtz_anonymous_23829684(x)

[S] = StandardState_anonymous_23817264(x);


function [S] = StandardState_anonymous_23817264(x)

[S]      = EquationOfState_anonymous_25942204(x);
i        = [3];
[S_1]    = MuT_cp_dippr_23808960(x(1));
S.g(1)   = S.g(1) + S_1.dmudT'*x(i);
S.g(i)   = S.g(i) + S_1.mu;
S.H(1,1) = S.H(1,1) + S_1.d2mudTdT'*x(i);
S.H(i,1) = S.H(i,1) + S_1.dmudT;
S.H(1,i) = S.H(i,1)';


function [S] = EquationOfState_anonymous_25942204(x)

[S]   = EquationOfState_anonymous_25857832(x);
[S_1] = ModTVN_ideal_idealgas_25941088(x);
S.g   = S.g + S_1.g;
S.H   = S.H + S_1.H;


function [S] = EquationOfState_anonymous_25857832(x)

[S] = ModTVN_pr_25856752(x);


function [S] = ModTVN_pr_25856752(x)

R                          = 8.314511984;       % [J/Kmol] Gasskonstant
d_1                        = 2.414213562;
d_2                        = 0.4142135624;
t_c                        = [304.21];          % [K] Kritisk punkt
p_c                        = [7383000.0];       % [Pa]
Omega_b                    = 0.0777960739;      % ok
T                          = x(1);
V                          = x(2);
i                          = [3];
n                          = x(i);
e                          = [1];               % Enhetsvektor (hjelpevektor)
c                          = 1/(d_1 + d_2);
NR                         = sum(n)*R;
RT                         = R*T;
NRT                        = NR*T;
dBdn                       = Omega_b*R*(t_c./p_c);
B                          = n'*dBdn;
C_1circ                    = (1/(V + d_1*B) - 1/(V - d_2*B))/B;
C_2circ                    = (-d_1/(V + d_1*B)^2 - d_2/(V - d_2*B)^2 - C_1circ)/B;
C_lnd                      = log((V + d_1*B)/(V - d_2*B));
C_lnv                      = log(V/(V - B));
C_1                        = (d_1/(V + d_1*B) + d_2/(V - d_2*B) - C_lnd/B)/B;
C_2                        = (-(d_1/(V + d_1*B))^2 + (d_2/(V - d_2*B))^2 - 2*C_1)/B;
[d2AdTdT,d2AdndT,d2AdndnT] = ModTVN_pr_a_soave_25855840(T,n,R,t_c,p_c);
dAdT                       = d2AdndT'*n/2;
dAdn                       = d2AdndnT*n;
A                          = dAdn'*n/2;
g_1                        = NR*C_lnv - c*dAdT*C_lnd/B;
g_2                        = -NRT*B/(V*(V - B)) + A/((V + d_1*B)*(V - d_2*B));
g_i                        = RT*C_lnv*e + (NRT/(V - B) - c*A*C_1)*dBdn - c*C_lnd/B*dAdn;
H_11                       = -c*d2AdTdT*C_lnd/B;
H_21                       = -B*NR/(V*(V - B)) + dAdT/((V + d_1*B)*(V - d_2*B));
H_i1                       = R*C_lnv*e + (NR/(V - B) - c*dAdT*C_1)*dBdn - c*C_lnd/B*d2AdndT;
H_22                       = -NRT/V^2 + NRT/(V - B)^2 + c*A*(1/(V + d_1*B)^2 - 1/(V - d_2*B)^2)/B;
H_i2                       = -RT*B/(V*(V - B))*e - (NRT/(V - B)^2 + c*A*C_2circ)*dBdn - c*C_1circ*dAdn;
H_ii                       = RT/(V - B)*(e*dBdn' + dBdn*e') + (NRT/(V - B)^2 - c*A*C_2)*dBdn*dBdn' - c*C_1*(dAdn*dBdn' + dBdn*dAdn') - c*C_lnd/B*d2AdndnT;
S.g                        = [g_1;g_2;g_i];
S.H                        = [H_11,H_21',H_i1';H_21,H_22,H_i2';H_i1,H_i2,H_ii];


function [d2AdTdT,d2AdndT,d2AdndnT] = ModTVN_pr_a_soave_25855840(T,n,R,t_c,p_c)

e        = [1];
Omega_a  = 0.4572355289;    % ok
K        = [0];             % Interaksjonsparameter. Null ved ren komp
r_c      = t_c./sqrt(p_c);
F        = 2*Omega_a*R^2*(r_c*r_c'.*(e*e' - K));
[m]      = ModTVN_pr_a_soave_m_pr_25854832;
k        = e + m.*(e - sqrt(T./t_c));
k_t      = -1/(2*T)*(m - k + e);
k_tt     = -1/(2*T)*k_t;
d2AdTdT  = (n.*k_tt)'*F*(k.*n) + (n.*k_t)'*F*(k_t.*n);
d2AdndT  = k_t.*(F*(k.*n)) + k.*(F*(k_t.*n));
d2AdndnT = F.*(k*k');


function [m] = ModTVN_pr_a_soave_m_pr_25854832

omega_az = [0.223621];      % Asentrisk faktor
%omega_az = [0.239];
%omega_az = [0.024];
e        = [1];
m        = 0.37464*e + 1.54226*omega_az + -0.2699*omega_az.^2; % ok


function [S] = ModTVN_ideal_idealgas_25941088(x)

R     = 8.314511984;
pcirc = [101325.0];
xcirc = [1.0];
T     = x(1);
V     = x(2);
i     = [3];
n     = x(i);
NR    = sum(n)*R;
NRT   = NR*T;
e     = [1];
p     = R*log(R*T*(n./pcirc)/V);
g_1   = n'*p;
g_2   = -NRT/V;
g_i   = T*p;
H_11  = NR/T;
H_21  = -NR/V;
H_i1  = p + R*e;
H_22  = NRT/V^2;
H_i2  = -R*(T/V)*e;
H_ii  = R*T*diag(e./n);
S.g   = [g_1;g_2;g_i];
S.H   = [H_11,H_21',H_i1';H_21,H_22,H_i2';H_i1,H_i2,H_ii];


function [S] = MuT_cp_dippr_23808960(T)

t_0        = [298.15];
t_min      = [50.0];
t_max      = [5000.0];
C          = [29.37,34.54,1428.0,26.4,588.0];
c_1        = C(:,1);
c_2        = C(:,2);
c_3        = C(:,3);
c_4        = C(:,4);
c_5        = C(:,5);
c_p        = c_1 + c_2.*c_3.^2./(T^2*sinh(c_3/T).^2) + c_4.*c_5.^2./(T^2*cosh(c_5/T).^2);
h          = c_1*T + c_2.*c_3.*cosh(c_3/T)./sinh(c_3/T) - c_4.*c_5.*sinh(c_5/T)./cosh(c_5/T) - c_1.*t_0 - c_2.*c_3.*cosh(c_3./t_0)./sinh(c_3./t_0) + c_4.*c_5.*sinh(c_5./t_0)./cosh(c_5./t_0);
s          = c_1*log(T) + c_4.*(log(cosh(c_5/T)) - c_5.*sinh(c_5/T)./cosh(c_5/T)/T) + c_2.*(c_3.*cosh(c_3/T)./sinh(c_3/T)/T - log(sinh(c_3/T))) - c_1.*log(t_0) - c_4.*(log(cosh(c_5./t_0)) - c_5.*sinh(c_5./t_0)./cosh(c_5./t_0)./t_0) - c_2.*(c_3.*cosh(c_3./t_0)./sinh(c_3./t_0)./t_0 - log(sinh(c_3./t_0)));
S.mu       = h - T*s;
S.dmudT    = -s;
S.d2mudTdT = -(c_p/T);
i          = [1];
[S_1]      = MuT_hs_h0_23796204(T);
S.mu(i)    = S.mu(i) + S_1.mu;
S.dmudT(i) = S.dmudT(i) + S_1.dmudT;


function [S] = MuT_hs_h0_23796204(T)

hcirc   = [-393510.0];              % Dannelsesentalpi
[S]     = MuT_hs_s0_23790012(T);
S.mu    = S.mu + hcirc;             
S.dmudT = S.dmudT + [0];


function [S] = MuT_hs_s0_23790012(T)

scirc   = [213.677];
S.mu    = -(T*scirc);
S.dmudT = -scirc;
