%% ME-446 Homework 5
clear all
close all
clc

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

%% Problem 1 - Part B
delta_p_star = linspace(0, 0.2, 100)';
sigma1 = 1;
sigma2 = 0.1;
omega = 32*pi/(32+9*pi);

% Defining j* as a function of delta_p_star for both cases
j_star_schrage = @(s, p) 2*s/(2-s) * p;
j_star_ytrehus = @(s, p) (omega*s)/(s+(1-s)*omega) * p;

% Plot the functions
fig1 = figure();

% Schrage (Mass flux)
plot(delta_p_star, j_star_schrage(sigma1, delta_p_star),...
    "Color", "k", "LineStyle", "-", "DisplayName", "Schrage, \sigma="+sprintf("%.1f", sigma1))
hold on
plot(delta_p_star, j_star_schrage(sigma2, delta_p_star),...
    "Color", "k", "LineStyle", "--", "DisplayName", "Schrage, \sigma="+sprintf("%.1f", sigma2))


% Ytrehus (Momentum and energy)
plot(delta_p_star, j_star_ytrehus(sigma1, delta_p_star), ...
    "Color", "r", "LineStyle", "-", "DisplayName", "Ytrehus, \sigma="+sprintf("%.1f", sigma1))
plot(delta_p_star, j_star_ytrehus(sigma2, delta_p_star),...
    "Color", "r", "LineStyle", "--", "DisplayName", "Ytrehus, \sigma="+sprintf("%.1f", sigma2))

l = legend("Location", "northwest");
xlabel('\Delta P^*')
ylabel('j^*') 

% Legend properties
l.FontSize = 7;
l.Box = 'off';
l.ItemTokenSize = [15,9];
l.NumColumns = 1;

% 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, 'HW5_P1_Plot.pdf', 'ContentType', 'vector');

%% Problem 2 - Part B
alpha = linspace(0.6, 0.9, 100)';   % Array of absorptivity [-]
q_sun = 1000;                       % Solar heat flux [W/m²]
h_loss = 15;                        % Heat loss heat transfer coefficient [W/(m².K)]
t_inf = 25;                         % Ambient temperature [°C]
hfg = 2200e3;                       % Enthalpy of vaporization of water [J/kg]
l = 10e-2;                          % Absorber characteristic dimension [m]
sh = 5;                             % Mass transfer Sherwood number [-]
p = 1;                              % Ambient pressure [atm]
phi = 0.3;                          % Ambient relative humidity [-]

% Emprical equations for ρ_sat and D_va (both function of temperature)
rho_sat = @(t) (6.335 + 0.6718*t - 2.0887*10^(-2)*t.^2 + 7.3095*10^(-4)*t.^3) / 1e3;
d_va = @(t) 1.87*10^(-10)*((t+273.15).^2.072/p);

% Energy balance as a function of α and T_surf
q_solar = @(a) q_sun * a;
q_loss = @(t) h_loss * (t - t_inf);
j_vd = @(t) (sh*d_va((t_inf+t)/2))/l .* (rho_sat(t)-phi*rho_sat(t_inf));
q_evap = @(t) j_vd(t) * hfg;
energy_balance = @(t, a) q_solar(a) - q_evap(t) - q_loss(t);

% Initial guess for the temperature
t_surf_0 = 30;
t_surf = arrayfun(@(a) fsolve(@(t) energy_balance(t, a), t_surf_0), alpha);

% Plot mass flux
fig2 = figure();
yyaxis left
plot(alpha, j_vd(t_surf)*1e3)
hold on

xlim([min(alpha), max(alpha)])
xlabel('\alpha [-]')
ylabel('Evaporation mass flux [g/m^2.s]') 

yyaxis right
plot(alpha, t_surf)

ylabel('Absorber temperature [°C]') 

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

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