% 10.34 - Fall 2006 HW#1
% Problem #2
% Equation Solvers
% Rob Ashcraft - Sept. 7th, 2006

% accompanying files (included in this function): fmin_eqn.m and fzero_eqn.m

function Problem2
% This function serves to set up the problem and constants used in solving
% the RK equation of state; generates a plot of the final solution vs. initial guess

% clear the workspace and command window
clear; clc;

% define certain global variables that are available to other functions
% without explicitly passing them (a global statement must also occur in
% the function where these are needed)
global psi omega sigma eps Tc Pc P R w;

% define the variable EOS parameters and physical constants
psi = 0.42748;
omega = 0.08664;
sigma = 1;
eps = 0;
Tc = 647.1;     % K
Pc = 217.7;     % atm
P = 1;          % atm
R = 0.08206;    % atm-L/mole-K
w = 0.152;      % acentric factor
T = 300;        % temp in K

b = omega*R*Tc/Pc;    % b parameter in the EOS (used in the plot)

V_0_vec = linspace(1,50,50);  %vector of initial guesses in L/mole

% solver options.  These can be very important with complex problems.  This
% particular option turns off the output display, so it does not have to
% print it to the command window (particularly annoying with iterating over fsolve)
options = optimset('Display','off'); 

% solve the equation using the three solvers specified
for i=1:length(V_0_vec)
    V_fzero_300K_gas(i) = fzero(@fzero_eqn,V_0_vec(i),options,T);

    V_fsolve_300K_gas(i) = fsolve(@fzero_eqn,V_0_vec(i),options,T);

    V_fminsearch_300K_gas(i) = fminsearch(@fmin_eqn,V_0_vec(i),options,T);
end

disp(['Vapor Volume Root (fzero, V_guess=50) (L/mole):  ',num2str(V_fzero_300K_gas(length(V_0_vec)))]);
disp(['Vapor Volume Root (fsolve, V_guess=50) (L/mole):  ',num2str(V_fsolve_300K_gas(length(V_0_vec)))]);
disp(['Vapor Volume Root (fminsearch, V_guess=50) (L/mole):  ',num2str(V_fminsearch_300K_gas(length(V_0_vec)))]);

% plot the solver solution as a function of initial guess
figure;
plot(V_0_vec,V_fzero_300K_gas,'.',V_0_vec,V_fsolve_300K_gas,'x',V_0_vec,V_fminsearch_300K_gas,'o',[0,50],[b,b],'--')
xlabel('Volume Initial Guess [L/mole]'); ylabel('Resulting Volume [L/mole]'); 
legend('fzero Result','fsolve Result','fminsearch Result','EOS b-value (lower physical limit)');



%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Include other functions in this file.
% This function defines the equation that should be equal to zero as a
% function of Volume and Temperature;  this is called when using the
% "fzero" or "fsolve" solver function

function resid = fzero_eqn(V,T)
%inputs:
%   V = specific volume in L/mole
%   T = temperature of the system in Kelvin
%outputs:
%   resid = residual of the EOS equation, which should -> 0 when solved
%           correctly

% this defines the global variables available to this function
global psi omega sigma eps Tc Pc P R w;

Tr = T/Tc;   %reduced temperature

% Define the EOS parameters
alpha = Tr^(-1/2);

a = psi*alpha*(R*Tc)^2/Pc;

b = omega*R*Tc/Pc;

% Equation to solve (set to zero).  Notice the P had to be moved to the RHS
% so that the equation would equal zero.  This is a recurring theme with
% many solvers
resid = -P + R*T/(V-b) - a/((V+eps*b)*(V+sigma*b));

return;



%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
% This function is used by fminsearch to find the optimum value of the
% variable that minimizes an objective function, Esqr in this case.  

function Esqr = fmin_eqn(V,T)
%inputs:
%   V = specific volume in L/mole
%   T = temperature of the system in Kelvin
%outputs:
%   Esqr = the square of the residual of the EOS equation, which should be 
%          minimized to 0 when solved correctly

global psi omega sigma eps Tc Pc P R w;

Tr = T/Tc;

alpha = Tr^(-1/2);

a = psi*alpha*(R*Tc)^2/Pc;

b = omega*R*Tc/Pc;

resid = -P + R*T/(V-b) - a/((V+eps*b)*(V+sigma*b));

% in this case, the plain residual cannot be used since a minimization is
% being performed, and not trying to set something to zero.  The easy way
% to solve this is to simply square the residual to make it always
% positive.  If you had multiple equations, you could use the sum of
% squares of errors as the residual in a similar manner.

Esqr = resid^2;

return;