Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
With MATLAB, I calculated the concentrations of Algae, Zooplankton, Oxygen, & Carbon in a river as they interact and are consumed.
%% Matthew Widjaja
% Environmental Modeling -- Program 1
% Parking Lot & Lake
clear t c
format short
global volumeV areaA baseFlowQ sediFlowM sediVelocityVs
%% Obtain Data
% Lets the user select from Assignment 1 Data or Alternate Data.
stepSizeN = input('Declare a Stepsize (ex. 0.1) = ');
sediDiameterD = input('Declare Diameter of Sediment (m) = ');
sediVelocityVs = (9.81 * ((2650/1000)-1) * sediDiameterD^2) ...
/ (18 * .000001103);
dataUse = input('Should Assignment 1 Data be used? (y/n) = ','s');
if dataUse == 'y'
volumeV = 200;
areaA = 20;
baseFlowQ = .2;
stormTime = 3600;
finalTime = 14400;
sediFlowM = 50;
else
volumeV = input('Declare Volume (m3) = ');
areaA = input('Declare Plan Area (m2) = ');
baseFlowQ = input('Declare Flow through Lake (m3/s) = ');
stormTime = input('Declare Duration of Storm (secs) = ');
finalTime = stormTime * 4;
sediFlowM = input('Delcare Sediment Flow during Storm (mg/s) = ');
end
%% Euler Variables
% Automates Variables for Eulers & Starting Parameters
initialTime = 1;
finalTime = finalTime + 1;
stepCountH = (finalTime-initialTime)/stepSizeN;
t(1) = 0;
c(1) = 0;
%% Special Conditions
% If sediDiameterD = 0, we set sediFlowM = 0.
if sediDiameterD <= 0;
sediFlowM = 0;
end
%% Test Steady State
% This uses a Formula to determine Theoratical Steady State
steadyStateC = (sediFlowM/volumeV) /...
((baseFlowQ/volumeV) + (sediVelocityVs/(volumeV/areaA)));
fprintf('Theoratical Steady State = %g\n',steadyStateC);
%% Euler Method
% This performs Eulers Method using the function in f.m
for i=initialTime:stepCountH;
if i <= stormTime/stepSizeN
c(i+1) = c(i) + (stepSizeN * f(c(i)));
t(i+1) = t(i) + stepSizeN;
else
sediFlowM = 0;
c(i+1) = c(i) + (stepSizeN * f(c(i)));
t(i+1) = t(i) + stepSizeN;
end
end
%% Data Presentation
% This produces a graph & confirms that the script is complete.
plot(t,c);
title('Sediment Concentration in a Lake');
xlabel('Time (Sec)');
ylabel('Concentration (mg/m3)');
fprintf('Process Completed\n\n');
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment