%% ME-446 Homework 6 - Problem 3
clear all
close all
clc

% Define a shorter alias for the CoolProp PropsSI Python function
% See http://www.coolprop.org/coolprop/HighLevelAPI.html#propssi-function
PropsSI = @py.CoolProp.CoolProp.PropsSI;

% Define figure properties
opts.Colors     = get(groot,'defaultAxesColorOrder');
opts.width      = 8;
opts.height     = 6;
opts.fontType   = 'Arial';
opts.fontSize   = 8;

% Parameters
T_l = 300;                                  % [K] Liquid temperature
R = 8.3145;                                 % [J/(mol.K)] Universal gas constant
M = 18.015e-3;                              % [kg/mol] Water molar mass 
R_s = R/M;                                  % [J/(kg.K)] Specific gaz constant
rho_l = PropsSI('D','T',300,'Q',0,'Water'); % [kg/m³] Water liquid density
P_sat = PropsSI('P','T',300,'Q',0,'Water'); % [Pa] Water saturation pressure at 300K
sigma = PropsSI('I','T',300,'Q',0,'Water'); % [N/m] Surface tension

% Range of pore size from 20 nm to 0.2 mm
r = 2*logspace(-8, -4, 100)';

% Equilibrium for vapor pressure in the pore (see Part A)
equilibrium = @(P_v, r) log(P_v/P_sat) - 1/(rho_l*R_s*T_l)*(P_v-P_sat-2*sigma/r);

% Initial guess for the pressure
P_v_0 = 3e-3;

% Solving
P_v = arrayfun(@(r) fsolve(@(P_v) equilibrium(P_v, r), P_v_0), r);

% Liquid pressure from Y-L
P_l = P_v - (2*sigma)./r;

% Plot
fig1 = figure();

% Plot Pv/Psat as a function of r
subplot(2,1,1);
semilogx(r, P_v/P_sat)
xlim([20e-9, 0.2e-3])
xlabel('Pore radius of curvature r [m]')
ylabel('P_{v}/P_{sat} [-]') 

% Plot Pl as a function of r
subplot(2,1,2);
semilogx(r, P_l*1e-5)  % Convert to bar
xlim([20e-9, 0.2e-3])
xlabel('Pore radius of curvature r [m]')
ylabel('P_{l} [bar]') 

% Figure properties
fig1.Units               = 'centimeters';
fig1.Position(3)         = opts.width;
fig1.Position(4)         = opts.height;
set(fig1.Children, ...
    'FontName',     opts.fontType , ...
    'FontSize',     opts.fontSize);

% Save the plot
exportgraphics(gcf, 'HW6_P3_Plot.pdf', 'ContentType', 'vector');