clear all

% Plot pressure of CS2-dimethoxymethane binary mixture to illustrate 
% Henry's law. 
%
% Original data from Zawidzki Zeitschrift fur Physicalische Chemie 1900.
% https://doi.org/10.1515/zpch-1900-3514
% 
% Executes properly on MATLAB version 9.4.0.813654 (R2018a).
% Executes mostly properly on OCTAVE version 4.2.1.
% Joshua C. Vaughan 1/31/2021

Xcs2 = [0 4.89 5.04 10.30 10.59 16.40 16.63 27.10 27.28 34.70 34.90...
    39.01 39.07 45.36 45.42 45.47 49.38 49.46 53.62 53.93 60.48 60.71...
    67.80 68.27 73.30 73.77 79.04 79.50 83.98 84.45 85.73 91.08 ...
    91.53 95.54 95.99 100]/100;

Xdmm = 1-Xcs2;

Pdmm = [0 54.5 55.6 109.3 112.9 159.5 163.2 234.8 237.0 277.6 279.4 ...
    297.8 299.5 324.8 325.6 327.3 340.4 340.2 357.2 357.2 380.0 381.9 ...
    404.8 407.0 423.4 424.3 440.8 442.6 456.4 458.1 462.2 481.8 483.4 ...
    501.0 500.0 514.5]*133.322; % Units of Pa

Pcs2 = [587.7 558.3 558.3 529.1 526.5 500.4 497.1 451.2 449.1 419.0 ...
    417.4 402.7 400.8 378.0 377.9 376.0 361.4 360.8 344.1 342.2 315.9 ...
    313.3 281.2 277.8 253.0 250.1 222.5 217.4 189.8 184.9 172.5 124.2 ...
    117.1 65.1 62.3 0]*133.322; % Units of Pa

Ptot = Pcs2 + Pdmm;

Xac = 1-Xcs2;

% Raoult's law prediction
n = 2^8;
Xa = linspace(0,1,n);
Pa_pure = Pdmm(end);
Pb_pure = Pcs2(1);

Pa = Xa*Pa_pure;
Pb = (1-Xa)*Pb_pure;
Ptot_raoult = Pa+Pb;

% create plot

ax_lw = 2;
pl_lw = 2;
ax_fs = 30;
la_fs = 28;
ti_fs = 30;
lg_fs = 20;
fontname = 'Arial';

figure(1); clf
plot(Xcs2,Pcs2/1e4,'bo-',Xcs2,Pdmm/1e4,'ro-','linewidth',pl_lw,'markersize',4)
hold on
plot([0 0.25],[0 4],'k--','linewidth',1.5*pl_lw) % Put this line in by eye.
plot([0.82 1],[4 0],'k--','linewidth',1.5*pl_lw) % Put this line in by eye.
plot(Xa,Pa/1e4,'r-',Xa,Pb/1e4,'b-','linewidth',1)
hold off


xlabel('{\itX}_{CS_2}','fontsize',la_fs)
ylabel('Pressure [\times10^4 Pa]','fontsize',la_fs)
set(gca,'fontsize',la_fs)
set(gca,'linewidth',ax_lw)
set(gca,'tickdir','out')
set(gcf,'color',[1 1 1])
set(gca,'fontname',fontname)
xlim([0 1])
ylim([0  1.1*max([Pa_pure Pb_pure])/1e4])
set(gca,'xtick',[0, 1])
%set(gca,'ytick',[0, Pb_pure, 1])
%set(gca,'yticklabel',{'0','{\itP}_B^*','{\itP}_A^*'})
set(gca,'LooseInset',get(gca,'TightInset'));

% add legend
h = legend('{\itP}_{DMM}','{\itP}_{CS2}','Henrys Law');
set(h,'fontsize',lg_fs);
set(h,'linewidth',1.5);
set(h,'position',[0.4972 0.2901 0.3200 0.2680])
set(h,'box','off')