% Find minimum DO level, time to minimum DO, and time to recovery.

% DO calculated based on the Streeter-Phelps Equation.



% Problem parameters

L = 30;     % pollutant load

c0 = 7;     % initial DO, mg/L 

cs = 8;     % saturation DO, mg/L

KD = .35;   % decay rate, /day

KR = .35;   % removal rate, /day

KA = .8;    % reaeration rate, /day



% Use brute-force algorithm:  set bounds for time (independent

% variable) and use uniform step size.



% Modeler-chosen parameters

deltaT= 0.5;        % time step, day

lBound= 0;          % lower bound for time (independent variable)

uBound= 12;         % upper bound for time (independent variable)

tolerance= 0.1;     % acceptable difference from target concentration



minSoFar= realmax;  % min DO conc. found so far (bogus initialization)

tmin= -1;           % time to min DO seen so far (bogus init.)



tt= 0:deltaT:uBound;

for t= 0:deltaT:uBound

    

    % Calculate c(t), print intermediate conc. for even t

    c = cs - L*KD/(KA-KR)*(exp(-KR*t)-exp(-KA*t)) - ...

        (cs-c0)*exp(-KA*t);

    if (mod(t,2)==0)

        fprintf('c(%.0f) = %.4f mg/L\n', t, c)

    end

    

    % Find min DO and time to min DO

    if (c<minSoFar)

        minSoFar= c;

        tmin= t;

    end

    

    % Find recovery time

    if (t>tmin && abs(c-c0)<tolerance)

        % condition t>tmin assures that c is climbing towards c0

        % and not dropping from c0.

        recoveryDay= t;

    end

end



fprintf('Lowest oxygen level, %0.2f mg/L, occurs after %.2f days.\n',...

        minSoFar, tmin)

fprintf('Oxygen level recovers to initial level after %.2f days.\n',...

        recoveryDay)