# modello di variazione della qualità dell'acqua sistema 0D
clear all
close all
pkg load odepkg
# costituenti di qualità dell'acqua simulati
# y(1)= DO
# y(2)=COD
# y(3)=alghe
# y(4)=ORG-N
# y(5)=NH3
# y(6)=NO2
# y(7)=NO3
# y(8)=Org-P
# y(9)=Diss-P
global k1 k2 k3 Ostar
global alfa1 alfa2 alfa3 alfa4 alfa5 alfa6
global mu_max rho KN KP
global sigma1 sigma4 sigma5 depth
global beta1 beta2 beta3 beta4
global F
# decay/settling rates
k1=1.;  % [0.02:3.5] 1/day COD degradation
k2=5.; % [0:100] 1/day oxygen reareation
k3=0.1; % [-3.6:3.6] 1/day COD settling
# stochiometric coefficients
alfa1=0.07 ; % [0.07:0.09] mg N/mg A
alfa2=0.01 ; % [0.01:0.02] mg P/mg A
alfa3=1.4; % [1.4:1.8] mg O/mg A oxygen production by algae
alfa4=1.6; % [1.6:2.3] mg O/mg A oxygen uptake by algae
alfa5=3.; % [3:4] mg O/mg/NH3
alfa6=1. ; % [1:1.14] mg O/mg NO2
# algal growth
mu_max=1.; % [1:3] 1/day max algal growth rate
rho=0.05; % [0.05:0.5] 1/day respiration rate
# saturation half constant
KN=0.01; % [0.01:0.3] mg N/L half saturation NH3+NO3
KP=0.001; % [0.001:0.05] mg P/L half saturation diss P
# settling rates
sigma1=0.174 ; % [0.174:1.83] m/day algae settling
sigma4=0.001 ; % [0.001:0.1] 1/day ORG-N settling
sigma5=0.001; % [0.001:0.1] 1/day ORG-P settling
# oxidation rates
beta1=0.10; % [0.1:1] 1/day NH3-> NO2
beta2=0.2; % [0.2:2] 1/day NO2->NO3
beta3=0.02; %[0.02:0.4] 1/day org_N -> NH3
beta4=0.01; % [0.01:0.7] 1/day org-P -> dissolved P
# ossigeno a saturazion
Ostar=8.; % mg/L @25°C 1atm
# altezza della corrente
depth=1.; % m
# fattore di preferenza delle alghe per NH3 o NO3
F=0.5; % [0:1]
# condizioni iniziali
y0(1)= 5.;
y0(2)=50.;
y0(3)=0.1;
y0(4)=0.001;
y0(5)=0.05;
y0(6)=0.01;
y0(7)=0.02;
y0(8)=0.05;
y0(9)=0.01;
# tempo di simulazione [day]
time=20.;
t=linspace(0,time,100);
options=odeset('AbsTol',1e-20,'RelTol',1e-5,'NonNegative',1:9,...
'InitialStep',0.05,'MaxStep',0.2);
[t,y]=ode45(@WQequations,[0 time],[y0],options);
figure(1)
plot(t,y(:,1),'r-',t,y(:,2),'b')
legend('DO','COD')
xlabel('time')
ylabel('BOD, DO')
figure(2)
plot(t,y(:,4),'r-',t,y(:,5),'b',t,y(:,6),'g',t,y(:,7),'p')
legend('ORG-N','NH3','NO2','NO3')
xlabel('time')
ylabel('concentrazioni')
figure(3)
plot(t,y(:,3),'r-')
legend('algae')
xlabel('time')
ylabel('concentrazioni')
figure(4)
plot(t,y(:,8),'b',t,y(:,9),'g')
legend('org-P','diss-P')
xlabel('time')
ylabel('concentrazioni')



