clear,clf
T1d    = linspace(300,1000,1050) + 273.15;% K
P1d    = linspace(10,55,1060);             % kbar
Cname  = {'Si','O'};                    % Chemical elements in the system
phases = {'q','coe'};                   % phases to use in calculations
[names,nphs,molm,td,data_name] = thermo_data(Cname,phases); % get thermodynamic data
% calculate gibbs energy of the pure phases
g          = Gref(T1d,P1d,td) - SdT_integral(T1d,P1d,td) + VdP_integral(T1d,P1d,td) ...
    +                              dG_landau(T1d,P1d,td) +     dG_order(T1d,P1d,td);
% Gibbs minimization
[gmin,ind] = min(g,[],3);
% Plotting
pcolor(T1d,P1d,ind),colorbar,shading flat,xlabel('T(K)'),ylabel('P(kb)'),title('my first phase diagram') % contour delta g
for i = 1:length(phases),gtext(phases{i},'Color','k','BackgroundColor','w'),end % this allows you to put labels manually
colorbar off,print -dpng -r150 Exercise2_1 % this prints a png file of the figure
%% Functions
function [phases,nphs,molm,td,data_name,Cname] = thermo_data(Cname,phases,aq_yes)
ndat = 25;nc = 19;np = 289;
fid = fopen('tc-ds633.bin','r');td = fread(fid,'double');fclose(fid);td = reshape(td,np,ndat);
fid = fopen('nphs-ds633.bin','r');nphs = fread(fid,'double');fclose(fid);nphs = reshape(nphs,np,nc);
fid = fopen('names-ds633.txt','r'); for ip = 1:np,names{ip,1} = fscanf(fid,'%s\n',1);end;fclose(fid);
fid = fopen('elements-ds633.txt','r'); for ip = 1:nc,elements{ip,1} = fscanf(fid,'%s\n',1);end;fclose(fid);
fid = fopen('td-colnames-ds633.txt','r'); for ip = 1:ndat,data_name{ip,1} = fscanf(fid,'%s\n',1);end;fclose(fid);
fid = fopen('elemass-ds633.bin','r');molm = fread(fid,'double');fclose(fid);
if ~isempty(Cname)
    for i_c = 1:length(elements)
        c_curr = find(strcmp(Cname,elements(i_c)));
        if ~isempty(c_curr)
            c_ind(i_c) = c_curr;
        else
            c_ind(i_c) = -1;
        end
    end
    ele_sys = find(c_ind>0);
    if sum(c_ind==-1)>0
        inc_phs_id  = find(sum(abs(nphs(:,c_ind==-1)),2)==0);
    end
    if exist('phases','var') && ~isempty(phases)
        for ip = 1:length(phases)
            if sum(strcmp(names(inc_phs_id),phases(ip)))==1
                p_id(ip) = inc_phs_id(strcmp(names(inc_phs_id),phases(ip)));
            else
                ele_phs = find(nphs(strcmp(names,phases(ip)),:)>0);
                ele_not_present = '';
                for i_ele = 1:length(ele_phs)
                    if sum(ele_sys==ele_phs(i_ele))==0
                        ele_not_present = [ele_not_present ', ' elements{ele_phs(i_ele)}];
                    end
                end
                error(['missing elements: ' ele_not_present])
            end
        end
    else
        p_id   = inc_phs_id;
        phases = names(inc_phs_id);
    end
    nphs                   = nphs(p_id,c_ind>0);
    nphs(:,c_ind(c_ind>0)) = nphs;
    molm                   = molm(c_ind>0);
    molm(c_ind(c_ind>0))   = molm;
    td                     = td(p_id,:);
    aq_id = ~isnan(td(:,21));
    if ~exist('aq_yes','var')
        aq_yes = 0;
    end
    if aq_yes == 0
        nphs(aq_id,:) = [];
        td(aq_id,:)   = [];
        phases(aq_id) = [];
    end
else
    for ip = 1:length(phases)
        if sum(strcmp(names,phases(ip)))==1
            p_id(ip) = find(strcmp(names,phases(ip)));
        else
            error('unknown phase')
        end
    end
    td    = td(p_id,:);
    nphs  = nphs(p_id,:);
    Cname = elements(sum(nphs,1)>0);
    molm  = molm(sum(nphs,1)>0);
    nphs  = nphs(:,sum(nphs,1)>0);
end
end
%% Gref
function Gr       = Gref(T1d,P1d,td)
Gr = zeros(numel(P1d),numel(T1d),size(td,1));
for ip = 1:size(td,1)
    Gr(:,:,ip) = ones(numel(P1d),numel(T1d))*(td(ip,1) - 298.15*td(ip,2))*1e3;
end
end
%% SdT integral
function SdT_int  = SdT_integral(T1d,P1d,td)
[P2d,T2d]     = ndgrid(P1d,T1d);T = T2d(:);P = P2d(:);
SdT_int = zeros(numel(T),size(td,1));
for ip = 1:size(td,1)
    Tr = 25 + 273.15;Sr = td(ip,2);a = td(ip,4);b = td(ip,5);c = td(ip,6);d = td(ip,7);    
    SdT_int(:,ip) =                      (T.*log(T/Tr) + Tr - T)*a ...
        +                              1/2*((T-Tr).^2)*b ...
        +                    ((T - Tr).^2./(2*T*Tr^2))*c ...
        +    (-4*sqrt(T) + 2*sqrt(Tr) + 2*T./sqrt(Tr))*d ...
        +                                       (T-Tr)*Sr;
end
SdT_int = reshape(SdT_int,length(P1d),length(T1d),size(td,1))*1e3;% convert to Joules
end
%% VdP integral
function VdP_int  = VdP_integral(T1d,P1d,td)
Pr = 1e-3;
Tr = 25 + 273.15;
[P2d,T2d] = ndgrid(P1d,T1d);T = T2d(:);P = P2d(:);
VdP_int = zeros(numel(P),size(td,1));
for ip = 1:size(td,1)
Hr = td(ip,1);Sr = td(ip,2);V1_298 = td(ip,3); a0   = td(ip,8);
k0     = td(ip,9);k0p    = td(ip,10);k0pp = td(ip,11);
dk0dt  = td(ip,22); sum_nphs = td(ip,23); isL = td(ip,24);fl_id = td(ip,25);
Gr     = (Hr-Tr*Sr)*1e3;
% VdP
theta =  10636/((Sr*1e3)/sum_nphs + 6.44);
u     = theta./T;
u0    = theta./Tr;
ksi0  = (u0.^2.*exp(u0))./((exp(u0)-1).^2);
if isL == 1
    pth   = 0;
    kT    = k0 + dk0dt.*(T-Tr);
    V1_T  = V1_298.*exp(a0.*(T-Tr));
else
    pth   = a0.*k0.*(theta./ksi0).*(1./(exp(u)-1)-1./(exp(u0)-1));
    kT    = k0;
    V1_T  = V1_298;
end
ta    = (1+k0p)./(1+k0p+kT*k0pp);
tb    = k0p./kT - k0pp/(1 + k0p);
tc    = (1 + k0p + kT .* k0pp)./(k0p + k0p.^2 - kT .* k0pp);
VdP = V1_T .* ((P - Pr).*(1 - ta) + ta.*(-(1 + tb .* (P - pth)).^(1 - tc) ...
    +           (1 + tb .* (Pr - pth)).^(1 - tc))./(tb.*(tc - 1)))...
    ./         ((1 - ta) + ta.*(1 + tb * Pr).^(-tc));
VdP(isnan(VdP)) = 0;
VdP = VdP*1e3;
% Fluids
if fl_id==1 || fl_id==2 || ~isnan(td(ip,21))
    SdT_int    = SdT_integral(T1d,P1d,td(ip,:));
    Pbar       = P*1e3; % convert to pressure in bar
    if ~isnan(td(21)),fl_id = 1;end
    if fl_id == 1
        molarmass = 18.015e3;
    elseif fl_id ==2
        molarmass = 44.0095e-3;
    end
    options = optimset('Display','off');
    % Find minimum G and corresponding rho
    rho_guess = [1.2 0.08 1.2/1000];
    gH2O      = zeros(size(T,1),length(rho_guess));
    VdPall    = zeros(size(T,1),length(rho_guess));
    VdP       = zeros(size(T));
    rho       = zeros(size(T));
    rho_kgm3  = zeros(size(T));
    for k = 1:length(rho_guess)
        rho_w = zeros(size(T));
        rho1  = rho_guess(k);
        for iT = 1:length(T)
            rho_w(iT) = fzero(@(rho) PS94(rho,T(iT),fl_id)-Pbar(iT)./83.14462./T(iT),rho1,options);
        end
        rho_kgm3(:,k)        = molarmass*rho_w;
        [~,~,~,lnf]          = PS94(rho_w,T,fl_id);
        VdPall(:,k)          = 8.314462*lnf(:).*T(:);
        gH2O(:,k)            = Gr - reshape(SdT_int,numel(T),1) + VdPall(:,k);
        gH2O(rho_w<0,k)      = 1000;
        gH2O(isnan(rho_w),k) = 1000;
    end
    % Find minimum
    [~,id] = min(gH2O,[],2);
    for k = 1:length(rho_guess)
        rho(id==k) = rho_kgm3(id==k,k);
        VdP(id==k) = VdPall(id==k,k);
    end
    V = 1./(rho/molarmass)/10; % J/bar/mol
elseif fl_id == 5 % O2
    VdP = 0;
elseif fl_id == 3 || fl_id == 4 || fl_id == 6 || fl_id == 7 || fl_id == 8    
    R_kJ  = 8.314/1e3;
    a0 = 5.45963e-5;
    a1 = -8.63920e-6;
    b0 = 9.18301e-4;
    c0 = -3.30558e-5;
    c1 = 2.30524e-6;
    d0 = 6.93054e-7;
    d1 = -8.38293e-8;
    if fl_id == 3 % CO (not implemented in THERMOCALC or Perple_X
        Tc = 132.9;
        Pc = 0.035;
    end
    if fl_id == 4 % CH4
        Tc = 190.6;
        Pc = 0.0460;
    end
    if fl_id == 6 % H2
        Tc = 41.2;
        Pc = 0.0211;
    end
    if fl_id == 7 % S2
        Tc = 1313.01;%??for S from the properties of gases and liquids Poling et al
        Pc = 182e-3;%??for S from the properties of gases and liquids Poling et al
    end
    if fl_id == 8 % H2S
        Tc = 373.40;%From the properties of gases and liquids Poling et al
        Pc = 89.63e-3;%From the properties of gases and liquids Poling et al
    end    
    a = a0*Tc^(5/2)/Pc + a1*Tc^(3/2)/Pc*T;
    b = b0*Tc/Pc;
    c = c0*Tc/(Pc^(3/2)) + c1/Pc^(3/2)*T;
    d = d0*Tc/Pc^2 + d1/Pc^2*T;
    VdP = (R_kJ*T.*log(1000*P) + b.*P + a./(b.*sqrt(T)).*( log(R_kJ*T + b.*P) - log(R_kJ*T + 2*b.*P) ) + 2/3*c.*P.*sqrt(P) + d/2.*P.^2)*1e3;
end
VdP_int(:,ip) = VdP;
end
VdP_int = reshape(VdP_int,numel(P1d),numel(T1d),size(td,1));
end
%% Landau
function g_landau = dG_landau(T1d,P1d,td)
[P2d,T2d] = ndgrid(P1d,T1d);T = T2d(:);P = P2d(:);
g_landau = zeros(numel(P),size(td,1));
% Landau corrected (see technical note)
Tr = 25 + 273.15;od = 0;
for ip = 1:size(td,1)
    Tc0 = td(ip,12);Smax = td(ip,13);Vmax   = td(ip,14);
    q2  = zeros(size(T,1),size(P,2));
    if ~isnan(Smax)
        q20       = sqrt(1-Tr./Tc0);
        tcp       = Tc0 + P.*(Vmax./Smax);
        temp_q2   = sqrt((tcp-T)./Tc0);
        q2(T<tcp) = temp_q2(T<tcp);
        g_exc     = Smax.*(Tc0.*(q20.*(1 - 1/3.*q20.^2) + 1/3*q2.^3) - q2.*tcp) - T.*Smax.*(q20 - q2) + P.*Vmax.*q20;
        if od == 1 % in case of ordered endmember ?
            g_exc     = Tc0.*Smax.*(-2/3 + q20.*(1 - q20.^2/3)) - Tc0.*Smax.*(q20 - 1) + P.*Vmax.*(q20 - 1);
        elseif od == 2 % in case of disordered endmember ?
            g_exc     = Tc0.*Smax.*q20.*(1 - q20.^2/3) - T.*Smax.*q20 + P.*Vmax.*q20;
        end
    else
        g_exc = 0;
    end
    g_landau(:,ip) = g_exc*1e3;
end
g_landau = reshape(g_landau,numel(P1d),numel(T1d),size(td,1));
end
%% Order
function g_order  = dG_order(T1d,P1d,td)
[P2d,T2d] = ndgrid(P1d,T1d);T = T2d(:);P = P2d(:);
g_order = zeros(numel(P),size(td,1));
for ip = 1:size(td,1)
    sfdh   = td(ip,15);sfdv = td(ip,16);sfw    = td(ip,17);sfwv = td(ip,18);sfn = td(ip,19);sffac = td(ip,20);od = 0;
    %% Symmetric formalism (see Holland and Powell 1996)
    if ~isnan(sfn)
        R = 8.3144/1e3;
        if od == 0 % equilibrium
            for iT = 1:length(T(:))
                q(iT,1) = symm_form(T(iT),sfdh+P(iT)*sfdv,sfw + P(iT)*sfwv,sfn,sffac);
            end
            h_sf = sfdh + P.*sfdv + q.*(sfw-sfdh+P*(sfwv-sfdv))-q.^2.*(sfw+P*sfwv); % (1-Q)*sfdh + Q*(1-Q)*sfw
            s_sf = -sffac*R*((1+sfn*q).*log((1+sfn*q)/(sfn+1))+sfn*(1-q).*log((1-q)/(sfn+1)) ...
                +  sfn*(1-q).*log(sfn*(1-q)/(sfn+1))+sfn*(sfn+q).*log((sfn+q)/(sfn+1)) )/(sfn+1);
            g_sf = h_sf - T.*s_sf;
        elseif od == 2 % disordered
            if sffac<0
                s_sf = sffac * R * (log(1/(sfn+1)) + sfn*log(sfn/(sfn+1)))*(1/sffac-sfn)/(sfn+1);
            else
                s_sf = sffac * R * (log(1/(sfn+1)) + sfn*log(sfn/(sfn+1)));
            end
            g_sf = sfdh + T.*s_sf + P.*sfdv;
        elseif od == 1 % ordered
            g_sf = 0;
        end
    else
        g_sf = 0;
    end
    g_order(:,ip) = g_sf*1e3;
end
g_order = reshape(g_order,numel(P1d),numel(T1d),size(td,1));
end
%% Bragg Williams
function f = bragg_will(dH,W,Q,fac,n,T)
R = 8.3144/1e3;
% From paper:
K = ((n-n*Q).*(1-Q))./((1+n*Q).*(n+Q));
f = dH + W.*(2.*Q-1) + fac.*n./(n+1).*R.*T.*log(K);
end
%% Symmetric formalism
function ihq = symm_form(T,dH,W,n,fac)

small = eps;
itlim = 50;

q1    = small;
q2    = 1-small;

f1 = bragg_will(dH,W,q1,fac,n,T);
f2 = bragg_will(dH,W,q2,fac,n,T);

if f2>0
    ihq = 1-small;
elseif f1<0
    ihq = small;
else
    qk = 0.5;
    fk = bragg_will(dH,W,qk,fac,n,T);
    for it = 1:itlim
        if fk>0
            q1 = qk;
        else
            q2 = qk;
        end
        qk = (q1+q2)/2;
        fk = bragg_will(dH,W,qk,fac,n,T);
        ihq = qk;
    end
end
end
%% Pitzer and Sterner EOS
function [P_RT,Pbar,A_nRT,lnf,dP_RTdr] = PS94(r,T,phs_id)
if phs_id ==1
    C = [
        0               0              0.24657688e+6  0.51359951e+2  0              0;             %1
        0               0              0.58638965e+0 -0.28646939e-2  0.31375577e-4  0;             %2
        0               0             -0.62783840e+1  0.14791599e-1  0.35779579e-3  0.15432925e-7; %3
        0               0              0             -0.42719878e+0 -0.16325155e-4  0;             %4
        0               0              0.56654978e+4 -0.16580167e+2  0.76560762e-1  0;             %5
        0               0              0              0.10917883e+0  0              0;             %6
        0.38878656e+13 -0.13494878e+9  0.30916564e+6  0.75591105e+1  0              0;             %7
        0               0             -0.65537898e+5  0.18810675e+3  0              0;             %8
        -0.14182435e+14  0.18165390e+9 -0.19769068e+6 -0.23530318e+2  0              0;             %9
        0               0              0.92093375e+5  0.12246777e+3  0              0];            %10
elseif phs_id == 2
    C = [
        0               0              0.18261340e+7   0.79224365e+2  0              0            ; %1
        0               0              0               0.66560660e-4  0.57152798e-5  0.30222363e-9; %2
        0               0              0               0.59957845e-2  0.71669631e-4 0.62416103e-8; %3
        0               0              -0.13270279e+1 -0.15210731e+0  0.53654244e-3  -0.71115142e-7; %4
        0               0              0.12456776e+0   0.49045367e+1  0.98220560e-2  0.55962121e-5;  %5
        0               0              0               0.75522299e+0  0              0;             %6
        -0.39344644e+12 0.90918237e+8  0.42776716e+6  -0.22347856e+2  0              0;             %7
        0               0              0.40282608e+3   0.11971627e+3  0              0;             %8
        0               0.22995650e+8 -0.78971817e+5  -0.63376456e+2  0              0;             %9
        0               0              0.95029765e+5   0.18038071e+2  0              0];            %10
end
% calculation of coefficients Pitzer & Sterner (1994)
c1  = C(1,1)*T.^(-4)  + C(1,2)*T.^(-2)  + C(1,3)*T.^(-1)  + C(1,4)  + C(1,5)*T  + C(1,6)*T.^2;
c2  = C(2,1)*T.^(-4)  + C(2,2)*T.^(-2)  + C(2,3)*T.^(-1)  + C(2,4)  + C(2,5)*T  + C(2,6)*T.^2;
c3  = C(3,1)*T.^(-4)  + C(3,2)*T.^(-2)  + C(3,3)*T.^(-1)  + C(3,4)  + C(3,5)*T  + C(3,6)*T.^2;
c4  = C(4,1)*T.^(-4)  + C(4,2)*T.^(-2)  + C(4,3)*T.^(-1)  + C(4,4)  + C(4,5)*T  + C(4,6)*T.^2;
c5  = C(5,1)*T.^(-4)  + C(5,2)*T.^(-2)  + C(5,3)*T.^(-1)  + C(5,4)  + C(5,5)*T  + C(5,6)*T.^2;
c6  = C(6,1)*T.^(-4)  + C(6,2)*T.^(-2)  + C(6,3)*T.^(-1)  + C(6,4)  + C(6,5)*T  + C(6,6)*T.^2;
c7  = C(7,1)*T.^(-4)  + C(7,2)*T.^(-2)  + C(7,3)*T.^(-1)  + C(7,4)  + C(7,5)*T  + C(7,6)*T.^2;
c8  = C(8,1)*T.^(-4)  + C(8,2)*T.^(-2)  + C(8,3)*T.^(-1)  + C(8,4)  + C(8,5)*T  + C(8,6)*T.^2;
c9  = C(9,1)*T.^(-4)  + C(9,2)*T.^(-2)  + C(9,3)*T.^(-1)  + C(9,4)  + C(9,5)*T  + C(9,6)*T.^2;
c10 = C(10,1)*T.^(-4) + C(10,2)*T.^(-2) + C(10,3)*T.^(-1) + C(10,4) + C(10,5)*T + C(10,6)*T.^2;
% ------------------------------------------------------------------------
P_RT =  r + c1.*r.^2 - r.^2.*((c3 + 2*c4.*r + 3*c5.*r.^2 + 4*c6.*r.^3)./...
    (c2 + c3.*r + c4.*r.^2 + c5.*r.^3 + c6.*r.^4).^2) ...
    + c7.*r.^2.*exp(-c8.*r) + c9.*r.^2.*exp(-c10.*r);
Pbar = P_RT.*8.314462.*T*1e6/1e5;

dP_RTdr = 1 + 2.*c1.*r - 2.*r.*(4.*c6.*r.^3 + 3.*c5.*r.^2 + 2.*c4.*r + c3)./(c6.*r.^4 + c5.*r.^3 + c4.*r.^2 + c3.*r + c2).^2 ...
    - r.^2.*(12.*c6.*r.^2 + 6.*c5.*r + 2.*c4)./(c6.*r.^4 + c5.*r.^3 + c4.*r.^2 + c3.*r + c2).^2 ...
    + 2.*r.^2.*(4.*c6.*r.^3 + 3.*c5.*r.^2 + 2.*c4.*r + c3).^2./(c6.*r.^4 + c5.*r.^3 + c4.*r.^2 + c3.*r + c2).^3 ...
    + 2.*c7.*r.*exp(-c8.*r) - c7.*r.^2.*c8.*exp(-c8.*r) + 2.*c9.*r.*exp(-c10.*r) - c9.*r.^2.*c10.*exp(-c10.*r);
A_nRT = 0*log(r) + c1.*r + ...
    (1./(c2 + c3.*r +c4.*r.^2 + c5.*r.^3 + c6.*r.^4) ...
    - 1./c2 ) - (c7./c8).*(exp(-c8.*r)-1) - (c9./c10) ...
    .* (exp(-c10.*r)-1);

lnf = log(r) + A_nRT + Pbar./(r.*83.14462.*T) + log(83.14462*T) -1;
end