% Hp_co2.m
% Kalkulerer p som funksjon av H ved ulike T

format short e;
warning('off');                 % Hvis warning: vil ikke synes

i  = 2:3;                       % Volum og moltall (plassering i matr)
n  = 3;                         % Moltall
t  = 220:20:440;                        % Temperatur
v = logspace(log10(3.9e-5),log10(1e-2),100);     % Volumer
x0 = [NaN NaN 1]';            % Initiell tilstandsvektor
e1 = [1,zeros(size(n))];        % Beskrankningsmatrise (stabilitetssøk)
enthalpy  = [];                 % Beregnede entalpier
pressure  = [];                 % Beregnede trykk
maxit     = 20;                 % Max antall iterasjoner(stabilitetssøk)

for T = t
    phs = 1;                      % Antar 1 fase
    h   = [];                     % Beregnede entalpier
    p   = [];                     % Beregnede trykk
    for V = v
        if phs==1
            x        = x0;                        % Default tilstandsvektor
            x(1:2)   = [T;V];                     % Legger inn T og V
            A        = co2_pr(x);                   % Beregner gradient og Hessian
            h(end+1) =-A.g(1)*T + A.g(n)'*x(n);   % Entalpi
            p(end+1) =-A.g(2);                    % Trykk
            H0   = A.H;                           % Trenger denne Hessiske hvis
                                                  % en mer stabil fase finnes
            g0   = A.g;                           % Tangentplan
            y    = [x(1:2);x(n)/10000];           % Crude estimate for vap phase. 
                                                  % Fortynner massen slik av vi helt 
                                                  % sikkert ikke har væskefase lenger
            dy   = y;                             % Initialize update vector
            iter = 0;                             % Restrict number of iterations 
                                                  % (stability calculation may  
                                                  % oscillate for unstable phases)
            while iter<maxit & norm(dy)>1e-6      % Do a stability check: Sjekker at
                                                  % vi ikke har væske
                iter = iter + 1;
                A    = co2_pr(y);
                tmp  = inv([e1',A.H(i,i);0,e1])*[g0(i)-A.g(i);0]; % Løsningsvektor
                dy   = [0;tmp(i)];
                s    = min(1,-0.8/min([dy./y;-0.3*dy./y]));
                y    = y + s*dy;
            end
            if iter==maxit
                disp(['Stability test did not converge (T=',num2str(T),',V=',num2str(V),')'])
            end
            if tmp(1)>1e-6*p(end)                  % Assume positive pressure 
                                                   % i.e. p(end)>0 in this test
                alpha  = min(0.1,tmp(1)*V/(y(i)'*H0(i,i)*y(i))); % Assume alpha>0 
                                                   % (maybe<0 for unstable phases)
                y(i)   = y(i)*alpha;               % Scale V and n in vapor phase
                x(i)   = x(i) - y(i);              % Apply V and n balances for 
                                                   % liquid phase
                dy     = y;                        % Initialize update vector
                while norm(dy)>1e-6                % Entrer tofase-flash
                    Vap = co2_pr(y);               % Forventer å finne to faser
                    Liq = co2_pr(x);
                    dy  = [0;inv(Liq.H(i,i)+Vap.H(i,i))*(Liq.g(i)-Vap.g(i))];
                    s   = min(1,-0.8/min([dy./y;-dy./x]));
                    y   = y + s*dy;
                    x   = x - s*dy;
                end
                if any(isnan([x,y]))               % Hvis vi har NaN --> noe er galt
                    h(end) = NaN;
                    p(end) = NaN;
                    disp(['Phase is metastable but flash did not converge (T=',num2str(T),',V=',num2str(V),')'])
                else
                    phs    = 2;                    % Assume two phases from now on
                    h(end) = -(Liq.g(1)+Vap.g(1))*T + Vap.g(n)'*(x(n)+y(n));  
                                                   % Same chem.pot. in both phases
                    p(end) = - Vap.g(2);           % Same pressure in both phases
                    disp(['Vapor-liquid 2-phase (T=',num2str(T),',V=',num2str(V),'),V/L=',num2str(sum(y(n))/sum(x(n)))])
                end
            end
        else
            y(1:2) = [T;V-x(2)];
            x(1:2) = [T;V-y(2)];
            dy     = y;
            while norm(dy)>1e-6 & norm(x(i))>1e-6
                Vap = co2_pr(y);              
                Liq = co2_pr(x);
                dy  = [0;inv(Liq.H(i,i)+Vap.H(i,i))*(Liq.g(i)-Vap.g(i))];
                s   = min(1,-0.8/min([dy./y;-dy./x]));
                y   = y + s*dy;
                x   = x - s*dy;
            end
            h(end+1) = -(Liq.g(1)+Vap.g(1))*T + Vap.g(n)'*(x(n)+y(n));  
                                                   % Same chem.pot. in both phases
            p(end+1) = - Vap.g(2);                 % Same pressure in both phases
            if norm(x(i))<1e-6
                phs    = 1;                        % Assume one (vapor) phase
                                                   % from now on
                x(i)   = x(i) + y(i);              % Apply V and n balances to 
                                                   % vapor phase
                h(end) =-Vap.g(1)*T + Vap.g(n)'*x(n);
                p(end) =-Vap.g(2);
                disp(['Single (vapor) phase (T=',num2str(T),',V=',num2str(V),')'])
            end
        end
    end
    enthalpy(:,end+1) = h';
    pressure(:,end+1) = p';
end
semilogy(enthalpy,pressure)
%axis([min(enthalpy), max(enthalpy)])
%title('Hp-diagram for CO2')
xlabel('Entalpi [J/mol]')
ylabel('Trykk [Pa]')
legend('220K', '240K', '260K', '280K', '300K', '320K', '340K', '360K', '380K', '400K', '420K', '440K')
return
