% Program to evaluate dissolved oxygen from Streeter-Phelps equation
clear all
close all

% Assign values for coefficients
Kd=0.8 % organic matter decay rate, day^(-1)
Ka=1.0 % aeration rate, day^(-1)
cs=8   % saturation concentration dissolved oxygen, mg/L
c0=6   % upstream concentration dissolved oxygen, mg/L
L0=10  % concentration of oxygen-consuming organic matter injetced at x=0, mg/L
U=4    % stream velocity km/day
% define vector of stream locations
x=[0:10] % km 
% evaluate K1, t, and c
K1=Kd/(Ka-Kd)
% t is a vector of travel times (in days)
% corresponding to locations specified in x
t=x/U
% c is a vector of concentrations corresponding
% to locations specified in x 
c=cs-K1*L0*(exp(-Kd*t)-exp(-Ka*t))-(cs-c0)*exp(-Ka*t);
% plot c vs. x
figure
plot(x,c)
title('Dissolved Oxygen in a Stream')
xlabel('Distance, km')
ylabel('Dissolved oxygen, mg/L')
axis([0 L0 0 cs])

% Using a for loop:
for x2=1:11
    t2=(x2-1)/U;
    c2(x2)=cs-K1*L0*(exp(-Kd*t2)-exp(-Ka*t2))-(cs-c0)*exp(-Ka*t2);
end

figure
plot(x,c2)
title('Dissolved Oxygen in a Stream')
xlabel('Distance, km')
ylabel('Dissolved oxygen, mg/L')
axis([0 L0 0 cs])