clear all
close all

nplank=3; 
load source.dat
source=source*1e-5; %mmolN/m3/s
source=source*60*60*24;
%mum=[6 .9]; %1/d
%kr=[0.005 0.01]; %mmolN/m3
%mum=[3 2];
%kr=[0.05 0.09];
mum=[8.825e-6 9.69993e-6 1.0066e-5]; mum=mum*24*60*60;
kr=[8.2646e-3 1.406e-2 2.3928e-2];
mp=[.1 .1 .1]; %1/d
grm=[10 10 10]; %m3/(mmolN d)
%grmh=[1 1]; %1/d
grmh=[2.639e-5 2.639e-5 2.6394-5]; grmh=grmh*24*60*60;
kp=[0.0125 0.0125 0.0125]; %mmolN/m3
ass=[.6 .6 .6]; % unitless
mz=[.1 .1  .1];  % 1/d
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
mucrit=(mum(2)*kr(1))/kr(2) + (mp(1)*(kr(2)-kr(1)))/(kr(2));
%%%%%%
%Sr=[1e-4:5e-4:1e-2 1.1e-2:1e-3:2e-1]; % 1e-3+1e-3:1e-3:.5e-2];  %mmolN/(m3 d)
Sr=[1e-4:5e-4:2e-2];
Rs=zeros(size(Sr)); Bs(:,1)=Rs; Bs(:,2)=Rs; Zs(:,1)=Rs; Zs(:,2)=Rs;
Rscheck=Rs;
%Sr=source;
%%%%%
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=(mum(1)*mz(1))/(ass(1)*grm(1));
srmu2crit=(mum(2)*mz(2))/(ass(2)*grm(2));
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));
Bs(fi,1)=Sr(fi)/mp(1);
Zs(fi,1)=0;
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=Rs(fi);
loss(fi,1)=mp(1);
loss(fi,2)=0;

% REGION 2
clear fi, fi=find(Sr>srcrit(1) & Sr<=srcrit(2) & Sr<srmu1crit);
Rs(fi)=(kr(1)*ass(1)*grm(1)*Sr(fi))./(mum(1)*mz(1)-ass(1)*grm(1)*Sr(fi));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Zs(fi,1)=(Sr(fi)*ass(1))/mz(1) - mp(1)/grm(1);
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));
loss(fi,1)=grm(1)*Zs(fi,1)+mp(1);
loss(fi,2)=0;
grow(fi,1)=mum(1);
grow(fi,2)=0;
growt(fi,1)=mum(1)*Rs(fi)/(Rs(fi)+kr(1));
growt(fi,2)=0;

clear f1i, f1i=find(Zs(:,1)>z1crit(2));

% REGION 3
clear fi, fi=find(Sr>srcrit(2) & Sr<=srcrit(3) & Sr<srmu1crit & Sr<srmu2crit);
Rs(fi)=(kr(2)*mp(2))/(mum(2)-mp(2));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Zs(fi,1)=AA;
Bs(fi,2)=squeeze(Sr(fi))/mp(2);
Bs(fi,2)=Bs(fi,2)- (mz(2)*squeeze(Zs(fi,1)))/(ass(2)*mp(2)) - mz(2)/(ass(1)*grm(1));
Zs(fi,2)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));
loss(fi,1)=grm(1)*Zs(fi,1)+mp(1);
loss(fi,2)=0;
grow(fi,1)=mum(1);
grow(fi,2)=0;
growt(fi,1)=mum(1)*Rs(fi)/(Rs(fi)+kr(1));
growt(fi,2)=0;

% REGION 4
clear fi, fi=find(Sr>srcrit(3)); % & Sr<srmu1crit & Sr<srmu2crit);
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Bs(fi,2)=mz(2)/(ass(2)*grm(2));
Zs(fi,1)=(ass(2)/mz(2))*Sr(fi)-mp(1)/grm(1)+mp(2)/grm(2);
Zs(fi,2)=(ass(2)*Sr(fi))/mz(2)-mp(1)/grm(1)-mp(2)/grm(2) + CC;
Rs(fi)=(kr(2)*(grm(2)*Zs(fi,2)+mp(2)))./(mum(2)-(grm(2)*Zs(fi,2)+mp(2)));
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));
loss(fi,1)=grm(1)*Zs(fi,1)+mp(1);
loss(fi,2)=grm(2)*Zs(fi,2)+mp(2);
grow(fi,1)=mum(1);
grow(fi,2)=mum(2);
growt(fi,1)=mum(1)*Rs(fi)/(Rs(fi)+kr(1));
growt(fi,2)=mum(2)*Rs(fi)/(Rs(fi)+kr(2));

% REGION 5
clear fi, fi=find(loss(:,1)>grow(:,1) & loss(:,2)<grow(:,2));
Rs(fi)=(kr(2)*ass(2)*grm(2)*Sr(fi))./(mum(2)*mz(2)-ass(2)*grm(2)*Sr(fi));
Bs(fi,1)=0;
Bs(fi,2)=mz(2)/(ass(2)*grm(2));
Zs(fi,1)=0;
Zs(fi,2)=(Sr(fi)*ass(2))/mz(2) - mp(2)/grm(2);
Rscheck(fi)=(kr(2)*(grm(2)*Zs(fi,2)+mp(2)))./(mum(2)-(grm(2)*Zs(fi,2)+mp(2)));
loss(fi,1)=NaN;

% REGION 5
clear fi, fi=find(loss(:,2)>grow(:,2) & loss(:,1)<grow(:,1));
Rs(fi)=(kr(1)*ass(1)*grm(1)*Sr(fi))./(mum(1)*mz(1)-ass(1)*grm(1)*Sr(fi));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Zs(fi,1)=(Sr(fi)*ass(1))/mz(1) - mp(1)/grm(1);
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));
loss(fi,2)=NaN;



clear fi, fi=find(loss(:,2)>grow(:,2) & loss(:,1)>grow(:,1));
Rs(fi)=NaN;
Bs(fi,1)=0;
Zs(fi,1)=0;
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=NaN;
loss(fi,:)=NaN;




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

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

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

subplot(4,1,4)
plot(Sr,squeeze(loss(:,1)),'b-','LineWidth',2);
hold on
plot(Sr,squeeze(loss(:,2)),'r--','LineWidth',2);
plot(Sr,squeeze(growt(:,1)),'b:','LineWidth',2);
plot(Sr,squeeze(growt(:,2)),'r:','LineWidth',2);
plot([Sr(1) Sr(length(Sr))],[mum(1) mum(1)],'k:','LineWidth',2);
plot([Sr(1) Sr(length(Sr))],[mum(2) mum(2)],'r:','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% numerical soluiont %%%%%%%%%&
fisr=[1:2:40 ];
%fisr=1:length(source);
subplot(4,1,1), plot(Sr(fisr),Rs(fisr),'mo','LineWidth',2);
delt=[.5/24]; ntim=50000;
npp=0;
for is=[fisr];
  npp=npp+1;
  if (npp>2), npp=1; end
  srr=Sr(is);
  rsv(1)=1e-5; %Rs(is);
  for ip=1:nplank
   bsv(1,ip)=1e-5; %Bs(is,ip);
   zsv(1,ip)=1e-5; %Zs(is,ip);
  end % for ip
  rsvh(1)=1e-5; %Rs(is);
  for ip=1:nplank
   bsvh(1,ip)=1e-5; %Bs(is,ip);
   zsvh(1,ip)=1e-5; %Zs(is,ip);
  end % for ip
  tim(1)=0;
  for it=1:ntim
    rr=squeeze(rsv(it));
    bb=squeeze(bsv(it,:));
    zz=squeeze(zsv(it,:));
%%
    rrh=squeeze(rsvh(it));
    bbh=squeeze(bsvh(it,:));
    zzh=squeeze(zsvh(it,:));
%
    gr(it)=0; gb(it,1:nplank)=0; gz(it,1:nplank)=0;
    grh(it)=0; gbh(it,1:nplank)=0; gzh(it,1:nplank)=0;
    for ip=1:nplank,
       gr(it)=gr(it)-mum(ip)*rr*bb(ip)/(rr+kr(ip));
       gb(it,ip)=gb(it,ip)+mum(ip)*rr*bb(ip)/(rr+kr(ip)) ...
                     -grm(ip)*zz(ip)*bb(ip)  ...
                     - mp(ip)*bb(ip);
       gz(it,ip)=gz(it,ip)+ass(ip)*grm(ip)*zz(ip)*bb(ip) ...
                     - mz(ip)*zz(ip);
    end % for ip 
    gr(it)=gr(it)+srr;
%%
    for ip=1:nplank,
       grh(it)=grh(it)-mum(ip)*rrh*bbh(ip)/(rrh+kr(ip));
       gbh(it,ip)=gbh(it,ip)+mum(ip)*rrh*bbh(ip)/(rrh+kr(ip)) ...
                     -grmh(ip)*zzh(ip)*(bbh(ip)/(bbh(ip)+kp(ip)))  ...
                     - mp(ip)*bbh(ip);
       gzh(it,ip)=gzh(it,ip)+ ...
                   ass(ip)*grmh(ip)*zzh(ip)*(bbh(ip)/(bbh(ip)+kp(ip))) ...
                     - mz(ip)*zzh(ip);
    end % for ip
    grh(it)=grh(it)+srr;
%%%
    rsv(it+1)=rr+gr(it)*delt;
    bsv(it+1,:)=bb+squeeze(gb(it,:))*delt; 
    zsv(it+1,:)=zz+squeeze(gz(it,:))*delt;
%%
    rsvh(it+1)=rrh+grh(it)*delt;
    bsvh(it+1,:)=bbh+squeeze(gbh(it,:))*delt;
    zsvh(it+1,:)=zzh+squeeze(gzh(it,:))*delt;
%%
    tim(it+1)=tim(it)+delt;
  end % for it
  for ipl=1:2,
   figure(ipl)
   subplot(4,1,1), plot(Sr(is),squeeze(mean(rsv(ntim-200:ntim))),'kx','LineWidth',2);
                   hold on
   subplot(4,1,2),
    for ip=1:nplank,
      if (ip==1), cstr=['bx']; end
      if (ip==2), cstr=['rx']; end
      if (ip==3), cstr=['co']; end
      plot(Sr(is),squeeze(mean(bsv(ntim-200:1:ntim,ip))),cstr,'LineWidth',2);
                   hold on
   end % for ip
   subplot(4,1,3)
    for ip=1:nplank,
      if (ip==1), cstr=['bx']; end
      if (ip==2), cstr=['rx']; end
      if (ip==3), cstr=['co']; end
      plot(Sr(is),squeeze(mean(zsv(ntim-200:1:ntim,ip))),cstr,'LineWidth',2);
                   hold on
   end % for ip
  end % for ipl

   figure(3)
   subplot(4,1,1), plot(Sr(is),squeeze(mean(rsvh(ntim-200:ntim))),'kx','LineWidth',2);
                   hold on
   subplot(4,1,2),
    for ip=1:nplank,
      if (ip==1), cstr=['bx']; end
      if (ip==2), cstr=['rx']; end
      if (ip==3), cstr=['co']; end
      plot(Sr(is),squeeze(mean(bsvh(ntim-200:1:ntim,ip))),cstr,'LineWidth',2);
                   hold on
   end % for ip
   subplot(3,1,3)
    for ip=1:nplank,
      if (ip==1), cstr=['bx']; end
      if (ip==2), cstr=['rx']; end
      if (ip==3), cstr=['co']; end
      plot(Sr(is),squeeze(mean(zsvh(ntim-200:1:ntim,ip))),cstr,'LineWidth',2);
                   hold on
   end % for ip



%%%%%%%
  if (npp==2),
   figure
   subplot(3,1,1)
   plot(tim, rsv,'k-','LineWidth',2);
   hold on
%  plot(tim, rsvh,'k-.','LineWidth',2);
   subplot(3,1,2)
   for ip=1:nplank
    if (ip==1), cstr='b'; end
    if (ip==2), cstr='r'; end
    if (ip==3), cstr='c'; end
     plot(tim, squeeze(bsv(:,ip)),[cstr,'-'],'LineWidth',2);
     hold on
%    plot(tim, squeeze(bsvh(:,ip)),[cstr,'-.'],'LineWidth',2);
   end % for
   subplot(3,1,3)
   for ip=1:nplank
    if (ip==1), cstr='b'; end
    if (ip==2), cstr='r'; end
    if (ip==3), cstr='c'; end
     plot(tim, squeeze(zsv(:,ip)),[cstr,'-'],'LineWidth',2);
     hold on
%    plot(tim, squeeze(zsvh(:,ip)),[cstr,'-.'],'LineWidth',2);
   end % for
  end % if npp
end % for is

