% FordamperOver.m
% Finner entalpidifferansen over fordamper
% Utgangspunkt: 3P CO2 SC IC
% h = T*s + my*N

close all

p1 = 48.30e5;           % [Pa] Gitt (før ventil)
temp = 273.15-30       % [K]   "
overheting = 51-33,19;
V1 = 4e-5;              % [m3] Startverdi, væskevolum
N = [1];                % En komponent (ren co2)
diffPR = [];
diffSW = [];
fordampertemp =[];

for T1 = temp
    T2 = 273.15-51;             % [K]  Konst T ved fordampning
    fordampertemp(end+1) = T2;
    % PENG-ROBINSON
    % Før fordamper:
    v1  = V1;
    dv1 = v1;
    % Newtoniterasjoner til V ved ønsket trykk finnes:
    while abs(dv1)>1e-6*v1
        a1 = co2_pr([T1,v1,N]');
        dv1 = (p1 + a1.g(2))/(-a1.H(2,2));
        v1 = v1 + dv1;
    end
    %disp('PR TILFELLE B: Overhetet damp')
    h1 =-a1.g(1)*T1 + a1.g(3)*N;   % Entalpi før fordamper

    % Etter fordamper, tørrmettet:
    pH_co2_pr;
    v2  = V;
    a2  = co2_pr([T2,v2,N]');
    p2 = -a2.g(2);
    % Overheting, konst p2:
    t2 = T2 + overheting;      % [K] Overhetet damp (Tevap=ca -10C ved 25bar)
    dv2= v2;
    while abs(dv2)>1e-6*v2
        a2 = co2_pr([t2,v2,N]');
        dv2 = (p2 + a2.g(2))/(-a2.H(2,2));
        v2 = v2 + dv2;
    end
    h2 = -a2.g(1)*t2 + a2.g(3)*N;
    % Differanse:
    Dh = h2 - h1
    diffPR(end+1) = Dh;

    % SPAN-WAGNER
    % Før fordamper:
    dV1 = V1;
    while abs(dV1)>1e-6*V1
        A1 = co2_sw([T1,V1,N]');
        dV1 = (p1 + A1.g(2))/(-A1.H(2,2));
        V1 = V1 + dV1;
    end
    H1 =-A1.g(1)*T1 + A1.g(3)*N;   % Entalpi før fordamper

    % Etter fordamper, tørrmettet:
    pH_co2_sw;
    V2 = V;
    A2  = co2_sw([T2,V2,N]');
    P2 = -A2.g(2);
    % Overheting:
    T2 = T2 + overheting;
    dV2= V2;
    while abs(dV2)>1e-6*V2
        A2 = co2_sw([T2,V2,N]');
        dV2 = (P2 + A2.g(2))/(-A2.H(2,2));
        V2 = V2 + dV2; 
    end
    H2 = -A2.g(1)*T2 + A2.g(3)*N;
    DH = H2 - H1
    diffSW(end+1) = DH;
end

avvik = (diffSW-diffPR)./diffSW*100

% subplot(2,1,1)
% diff = [diffPR; diffSW];
% plot(fordampertemp,diff)
% title('Entalpidifferanse over fordamper, overheting: 20 K')
% xlabel('Temperatur i fordamper [K]')
% ylabel('Entalpidifferanse [J/mol]')
% legend('Peng Robinson', 'Span Wagner')
% grid
% 
% subplot(2,1,2)
% plot(fordampertemp,avvik,'red')
% legend('Avvik for Peng Robinson relativt til Span Wagner')
% xlabel('Temperatur i fordamper [K]')
% ylabel('Prosentvis avvik [%]')
% grid