clear,clf
T1d    = linspace(300,1000,20) + 273.15;% K
P1d    = linspace(1,100,21);             % 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);
% Reactions
v  = null(nphs','r'); % find coefficients of reaction
dg = zeros(size(g,1),size(g,2));
for iphs = 1:length(phases)
    dg = dg + v(iphs)*g(:,:,iphs); % do products - reactants
end
% Plotting
contour(T1d,P1d,dg,[0,0],'k'),xlabel('T(K)'),ylabel('P(kb)'),title('my first reaction') 
%% 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