## Copyright (C) 2017 Marina Campolo
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn {Function File} {@var{retval} =} SP (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

## Author: Marina Campolo <marina@boopsie>
## Created: 2017-04-26

function dydt=WQequations (t,y)
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
dydt=zeros(9,1);
N_tot=y(5)+y(7);
mu=mu_max*N_tot/(N_tot+KN)*y(9)/(y(9)+KP);
# y(1)= DO
dydt(1)=k2*(Ostar-y(1))-k1*y(2)+(alfa3*mu-alfa4*rho)*y(3)-alfa5*beta1*y(5)...
   -alfa6*beta2*y(6);
# y(2)=COD
dydt(2)=-k3*y(2)-k1*y(2);
# y(3)=alghe
dydt(3)=mu*y(3)-rho*y(3)-sigma1/depth*y(3);
# y(4)=ORG-N
dydt(4)=alfa1*rho*y(3)-beta3*y(4)-sigma4*y(4);
# y(5)=NH3
dydt(5)=beta3*y(4)-beta1*y(5)-alfa1*F*mu*y(3);
# y(6)=NO2
dydt(6)=beta1*y(5)-beta2*y(6);
# y(7)=NO3
dydt(7)=beta2*y(6)-alfa1*mu*(1.-F)*y(3);
# y(8)=Org-P
dydt(8)=-beta4*y(8)-sigma5*y(8)+alfa2*rho*y(3);
# y(9)=Diss-P
dydt(9)=beta4*y(8)-alfa2*mu*y(3);
endfunction
