clear all
close all

mum=[.6 .9]; %1/d
kr=[0.005 0.01]; %mmolN/m3
mp=[.1 .1]; %1/d
grm=[5 5 ]; %m3/(mmolN d)
ass=[.6 .6]; % unitless
mz=[.1 .1];  % 1/d
%%%%%%%
mucrit=(mum(2)*kr(1))/kr(2) + (mp(1)*(kr(2)-kr(1)))/(kr(2));
%%%%%%
Sr=[1e-4:1e-4:1.5e-2];  %mmolN/(m3 d)
%%%%%
AA=(mp(1)/grm(1)) * (mum(1)*kr(2) - mum(2)*kr(1) - mp(1)*(kr(2)-kr(1)))/ ...
              (kr(1)*mum(2)+mp(1)*(kr(2)-kr(1)));
CC=AA; %???
srcrit(1)=(mp(1)*mz(1))/(ass(1)*grm(1));
srcrit(2)=srcrit(1)+(mz(1)*AA)/ass(1);
srcrit(3)=srcrit(1)+(mp(2)*mz(2))/(ass(2)*grm(2))+(mz(1)*CC)/ass(1);
srmu1crit(1)=100;
srmu1crit(2)=(mum(1)*mz(1))/(ass(1)*grm(1));
srmu1crit(3)=100;
srmu1crit(4)=100;
z1crit(1)=0;
z1crit(2)=(mum(1)+mp(1))/grm(1);
%%%%%
% REGION 1
clear fi, fi=find(Sr<srcrit(1));
Rs(fi)=(kr(1)*mp(1))./(mum(1)-mp(1));
B1(fi)=Sr(fi)/mp(1);
Z1(fi)=0;
B2(fi)=0;
Z2(fi)=0;
Rscheck(fi)=Rs(fi);

% REGION 2
clear fi, fi=find(Sr>srcrit(1) & Sr<=srcrit(2) & Sr<srmu1crit(2));
Rs(fi)=(kr(1)*ass(1)*grm(1)*Sr(fi))./(mum(1)*mz(1)-ass(1)*grm(1)*Sr(fi));
B1(fi)=mz(1)/(ass(1)*grm(1));
Z1(fi)=(Sr(fi)*ass(1))/mz(1) - mp(1)/grm(1);
B2(fi)=0;
Z2(fi)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Z1(fi)+mp(1)))./(mum(1)-(grm(1)*Z1(fi)+mp(1)));

clear f1i, f1i=find(Z1>z1crit(2));

% REGION 3
clear fi, fi=find(Sr>srcrit(2) & Sr<=srcrit(3));
Rs(fi)=(kr(2)*mp(2))/(mum(2)-mp(2));
B1(fi)=mz(1)/(ass(1)*grm(1));
Z1(fi)=AA;
B2(fi)=Sr(fi)/mp(2) - (mz(2)*Z1(fi))/(ass(2)*mp(2)) - mz(2)/(ass(1)*grm(1));
Z2(fi)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Z1(fi)+mp(1)))./(mum(1)-(grm(1)*Z1(fi)+mp(1)));

% REGION 4
clear fi, fi=find(Sr>srcrit(3));
B1(fi)=mz(1)/(ass(1)*grm(1));
B2(fi)=mz(2)/(ass(2)*grm(2));
Z1(fi)=CC;
Z2(fi)=(ass(2)*Sr(fi))/mz(2)-mp(1)/grm(1)-mp(2)/grm(2) - CC;
Rs(fi)=(kr(2)*(grm(2)*Z2(fi)+mp(2)))./(mum(2)-(grm(2)*Z2(fi)+mp(2)));
Rscheck(fi)=(kr(1)*(grm(1)*Z1(fi)+mp(1)))./(mum(1)-(grm(1)*Z1(fi)+mp(1)));

%%%%%%%%%%%%%%%%%%
figure
subplot(3,1,1)
plot(Sr,Rs,'b-','LineWidth',2);
hold on
plot(Sr,Rscheck,'r--','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end % for

subplot(3,1,2)
plot(Sr,B1,'k-','LineWidth',2);
hold on
plot(Sr,B2,'r-','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end 

subplot(3,1,3)
plot(Sr,Z1,'k-','LineWidth',2);
hold on
plot(Sr,Z2,'r-','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end 

