clear all

% Plot constant pressure heat capacity vs temperature for chlorine gas as well as absolute entropy.
% 
% Executes properly on MATLAB version 9.4.0.813654 (R2018a).
% Executes mostly properly on OCTAVE version 4.2.1.
% Joshua C. Vaughan 1/3/2021

mp = 172.12; % melting point
bp = 239.05; % boling point

% NIST data from Chase 1998 (Shomate equation)
% https://webbook.nist.gov/cgi/cbook.cgi?ID=C7782505&Mask=7
A = 33.05;
B = 12.23;
C = -12.07;
D = 4.39;
E = -0.16;

T_gas  = bp:1:500;

C_gas = A + B*(T_gas/1000) + C*(T_gas/1000).^2 + D*(T_gas/1000).^3 + E./((T_gas/1000).^2);

% Solid and liquid data from Giauque and Powell, JACS 1939. giauque_chlorine_1939
% Note, JV estimated heat capacity: lowest three values for solid, and
% added one data point for each of solid/liquid at the melting point and for
% liquid/gas at the boiling point. 
T_sol = [0 3 8 14.05 15.21 17.40 19.81 22.93 26.37 29.88 33.94 37.68 42.37 ...
    47.14 58.59 65.42 70.50 75.25 79.71 84.06 88.66 93.47 98.06 ...
    102.92 108.06 112.99 118.16 123.53 128.67 134.06 139.67 145.18 ... 
    150.39 155.45 160.31 164.99 mp]; % in K

T_liq = [mp 178.87 183.75 188.40 193.21 194.08 201.50 208.63 215.90 222.82 229.73 236.77 bp]; % in K

C_sol = 4.18*[0 0.06 0.2 0.810 0.916 1.331 1.842 2.467 3.192 3.946 4.804 5.421 6.018 6.675 ...
    7.879 8.382 8.720 8.973 9.201 9.427 9.626 9.827 10.03 10.22 10.38 ...
    10.57 10.81 11.00 11.24 11.47 11.72 11.92 12.24 12.41 12.69 12.93 13.2]; % in Joule/(K mol)

C_liq = 4.18*[16.0 16.01 16.02 15.99 15.99 15.95 15.94 15.92 15.85 15.81 15.79 15.73 15.7]; % in Joule/(K mol)
    
% create plot of heat capacity vs temperature.
ax_lw = 2;
pl_lw = 3;
ax_fs = 30;
la_fs = 30;
ti_fs = 30;
fontname = 'Arial';

figure(1); clf

plot(T_sol,C_sol,'r','linewidth',pl_lw)
hold on
plot(T_liq,C_liq,'g','linewidth',pl_lw)
plot(T_gas,C_gas,'b','linewidth',pl_lw)
hold off

xlabel('Temperature [K]','fontsize',la_fs)
%ylabel('{\itC_P} [J/(K mol)]','fontsize',la_fs)
ylabel('$\overline{C}_P$ \textrm{[J/(K mol)]}','interpreter','latex','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 400])
ylim([0 70])
set(gca,'LooseInset',get(gca,'TightInset'));

hold on
plot([mp mp],[1 80],'k--','linewidth',pl_lw)
plot([bp bp],[1 80],'k--','linewidth',pl_lw)
hold off

ypos = 15;
an_fs = 30;
text(90,ypos,'solid','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','r')
text(206,ypos,'liq.','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','g')
text(320,ypos,'gas','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','b')

% now do calculation for absolute entropy in five parts. values, below from Wikipedia.

% Interpolate heat capacity at smaller increments for better numerical integration.
num_pts = 200;
T_sol2 = linspace(T_sol(1),T_sol(end),num_pts);
C_sol2 = interp1(T_sol,C_sol,T_sol2);
T_liq2 = linspace(T_liq(1),T_liq(end),num_pts);
C_liq2 = interp1(T_liq,C_liq,T_liq2);
T_gas2 = linspace(T_gas(1),T_gas(end),num_pts);
C_gas2 = interp1(T_gas,C_gas,T_gas2);


deltaH_fus = 6406;      % J/mol
deltaH_vap = 20410;     % J/mol

S1 = [0 cumsum(C_sol(2:end).*diff(T_sol)./T_sol(2:end))];
S2 = S1(end) + deltaH_fus/mp + [0 cumsum(C_liq(2:end).*diff(T_liq)./T_liq(2:end))];
S3 = S2(end) + deltaH_vap/bp + [0 cumsum(C_gas(2:end).*diff(T_gas)./T_gas(2:end))];

S = [S1 S2 S3];
T = [T_sol T_liq T_gas];

figure(2); clf
hold on
plot([mp mp],[1 300],'k--','linewidth',pl_lw)
plot([bp bp],[1 300],'k--','linewidth',pl_lw)
plot(T,S,'m','linewidth',pl_lw)
hold off
xlabel('Temperature [K]','fontsize',la_fs)
ylabel('$S$ \textrm{[J/(K mol)]}','interpreter','latex','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 400])
ylim([0 250])
set(gca,'LooseInset',get(gca,'TightInset'));
set(gca,'box','on')

ypos = 175;
an_fs = 30;
text(90,ypos,'solid','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','r')
text(206,ypos,'liq.','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','g')
text(320,ypos,'gas','fontsize',an_fs,'horizontalalignment','center','fontweight','bold','color','b')
