% problem 4: a through e
clc, clear

% data
hi = [6,5,4,3,2,1,0];
Ai = [1.17,0.97,0.67,0.45,0.32,0.18,0].*10.^4;

% 4a: 5th degree least squares polynomial regression
regression_coef = polyfit(hi,Ai,5);
p5_regression = @(h) polyval(regression_coef, h);

% 4b: plot least squares regression & data points
f1 = figure();
h_vals = linspace(hi(1), hi(end), 100);
scatter(hi, Ai, "o", "LineWidth", 3, "MarkerEdgeColor", "r", "MarkerFaceColor", "r");
hold on;
plot(h_vals, p5_regression(h_vals), 'LineWidth', 3, 'Color', 'b');
set(gca, "FontSize", 12);
grid;
title("\textbf{5th Degree Polynomial Regression}",...
    "FontSize", 16, "Interpreter", "latex");
xlabel("$h$-values [m]", "FontSize", 16, "Interpreter", "latex");
ylabel("$A$-values [m\textsuperscript{2}]", "FontSize", 16, "Interpreter", "latex");
hold off;

% save figure
%saveas(f1, "4b_5th_degree_poly_regression","fig");
%saveas(f1, "4b_5th_degree_poly_regression","epsc");

% 4c: plug into ODE and solve via ode45
% define parameters & inital condition
h0 = 6;
d = 0.25;
e = 1;
g = 9.81;
time_span = [0,19].*3600;

h_dot = @(t,h) -pi.*d.^2/(4.*p5_regression(h)).*sqrt(2.*g.*(h+e));

% solve using ode45
[t_sol, h_sol] = ode45(h_dot, time_span, h0);

% plot
f2 = figure();
plot(t_sol, h_sol, 'Linewidth', 3, 'Color', 'r');
set(gca, "FontSize", 12);
grid;
title("\textbf{ODE45 Solution for Pond Drainage}",...
    "FontSize", 16, "Interpreter", "latex");
xlabel("time $t$ [s]", "FontSize", 16, "Interpreter", "latex");
ylabel("Depth $h(t)$ [m]", "FontSize", 16, "Interpreter", "latex");

% save figure
%saveas(f2, "4c_ode45_solution","fig");
%saveas(f2, "4c_ode45_solution","epsc");


% 4d: find t* s.t. h(t*) is about 0.01
% define threshold
epsilon = 0.01;

% find index where h_sol first hits epislon
index = find(h_sol <= epsilon, 1, 'first');

% corresponding t*
t_star = t_sol(index);

% display result
fprintf('The pond is almost empty at t* = %.2f seconds (%.2f hours).\n', t_star, t_star/3600);


% 4e: analytically solve the nonlinear ODE
% construct A(h) as symbolic expression
syms h(t) 
A_expr = 0;
degree = length(regression_coef) - 1;
for i = 1:length(regression_coef)
    A_expr = A_expr + regression_coef(i)*h(t)^degree;
    degree = degree - 1;
end

% define the ODE & initial condition
ode = diff(h,t) == -pi*d^2/(4*A_expr)*sqrt(2*g*(h+e));
cond = h(0) == 6;

% try to solve
hSol(t) = dsolve(ode, cond);

% display solution
disp(hSol)