% Green-Ampt formula
% exercise 7.4.3 (Mays book)

% Given initial conditions:
se          = 0.2;              % Initial effective saturation [-]
% For a silty clay soil (from table):
K           = 0.05;             % Hydraulic conductivity [cm/h]
theta_e     = 0.423;            % Effective saturation [-]
psi         = 29.22;            % Wetting front soil suction head [cm]
n           = 0.479;            % Porosity [-]

Delta_theta = (1-se)*theta_e;   % Difference in moisture across the wetting front

% Time
t=linspace(0,6,61);             % Time [h]

% Numerical iterations
tol0   = 10;                    % Initial value for the difference between subsequent approximations
thresh = 10^-5;                 % Threshold on the tolerance (note how results change by modifying this)
F      = zeros(size(t));        % Preallocation (Cumulative infilitration F)
f      = zeros(size(t));            % preallocation (Infiltration rate f)

for i=1:length(t)
    tol  = tol0; 
    Fold = K*t(i);                % A good trial value is a good trial value is F=Kt
    while abs(tol)>thresh
        Fnew = K*t(i)+psi*Delta_theta*log(1+Fold/(psi*Delta_theta));
        tol  = Fnew-Fold;
        Fold = Fnew;
    end
    F(i) = Fnew;
    f(i) = K*((psi*Delta_theta)/F(i)+1);
end

% Plot results
figure(1)
plot(t(2:end),F(2:end),'-k','LineWidth',1.5)
hold on
plot(t(2:end),f(2:end),'--r','LineWidth',1.5)
xlabel('Time [hr]','Interpreter','latex')
ylabel('$F$ [cm], $f$ [cm/hr]','Interpreter','latex')
legend('$F$','$f$','Interpreter','latex')