%This is a function for resource dynamics of nospec species and one resource
function res1_plot_para(specpara,supply,loss,paratochange,inival,endval,stepval,animation)
% species parameters:
% 1 growth: growth per unit biomass (u)
% 2 uptake: uptake per unit growth (w)
% 3 mortality (d)
% 4 initial size (n0)

%resource parameters:
% 5          supply: resource supply (S)
% 6          loss: rate at which resource is removed (a)

%further parameters:
% paratochange: number of the para to be changed (see above)
% inival: initial parameter value
% endval: final parameter value
% stepval: parameter change step
% animation - whether to animate

paratochange = fix(paratochange);
if (paratochange < 1 | paratochange > 6 | paratochange == 4)
   'Bad parameter code'
   return;
end

paraname = ['growth   ';'uptake   ';'mortality';'         ';'supply   ';'loss     ']';

para(1) = specpara(1,1); 
para(2) = specpara(1,2); 
para(3) = specpara(1,3); 
para(5) = supply;
para(6) = loss;

para(paratochange) = inival;
Nstarini = (para(5)- para(6)*para(3)/para(1)) / (para(2) * para(3));
Rstarini = para(3)/para(1);
Rfreeini = para(5)/para(6);
para(paratochange) = endval;
Nstarend = (para(5)- para(6)*para(3)/para(1)) / (para(2) * para(3));
Rstarend = para(3)/para(1);
Rfreeend = para(5)/para(6);

maxx = 1.2*max([Rstarini Rstarend Rfreeini Rfreeend]);
maxy = 1.2*max([50 Nstarini Nstarend]);
nobin = 100;
   
for value = inival:stepval:endval
   para(paratochange) = value;

   hold on;
   Nstar = (para(5)- para(6)*para(3)/para(1)) / (para(2) * para(3));
   Rstar = para(3)/para(1);
   Rfree = para(5)/para(6);
   text = sprintf('%s %7.2g  Rfree %7.2f  Rstar %7.2f   Nstar  %7.2f', paraname(:,paratochange), value, Rfree,Rstar,Nstar);
   plot([Rstar, Rstar],[0 maxy]);
   axis([0, maxx, 0, maxy])
   xlabel('Resource');
   ylabel('Species');
   title(text);
   
   xhyp=[];
   yhyp=[];
   for t=para(5)/(para(6)*nobin):para(5)/(para(6)*nobin):para(5)/para(6);
      xhyp=[xhyp t];
      yhyp=[yhyp (para(5)- para(6)*t)/(para(2)*para(1)*t)];
   end
   plot(xhyp,yhyp,'g');
   if Rstar>0
      plot(Rstar,Nstar,'r+');
   end
   
   if (animation)
      pause(1);
   end
   
end