% cobweb.m
%
% Grafico Cobweb para um modelo discreto de uma populaçao 
%
% O utilizador deve introduzir a equaçao que define o modelo.
% Depois de clicar na populaçao inicial, o  "cobweb" das futuras
% populaçoes e construido.
%
%
% 8/2/03

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p=0;                                              % initialize population
%                             
disp ('Introduza a fórmula que define o modelo,usando a notação "p"')
disp ('para o tamanho da espécie: (Exemplo: next_p = p+.2*p*(1-p/10))')
next_p=input ('next_p = ','s');
if isempty(next_p) next_p='p+1.8*p*(1-p/10)'; end;
eval( [next_p ';']);                              % test formula
%  
disp (' ');
disp ('Introduza os limites para "p" que serão usados na construção do gráfico:')
disp ('na construção do gráfico:')
limits=input('(por defeito é considerado  [pmin pmax]=[0 20])      [pmin pmax]= ');
if isempty(limits) limits=[0 20]; end;
pinc=(limits(2)-limits(1))/50;
x=limits(1):pinc:limits(2);
%
p=limits(1); y=eval (next_p);
for i=x(2):pinc:limits(2);                 % loop to create population vector
   p=i;                                        
   p=eval (next_p); 
   y=[y p]; 
end;
%
%
figure                                     % set up graphics
plot(x,y,x,x)
axis([limits(1),limits(2),limits(1),limits(2)]);    
xlabel ('population at time "t" ');
ylabel ('population at time "t+1"');
title (['next\_p=',next_p]);
%
continueb=1;                               % set boolean continue value to true
while continueb                            % loop until nonleft button is hit
   [p,q,button]=ginput(1);                 % get initial population
   if button==1
      %
      plot (x,y,x,x);                      %sets up graphics  
      axis([limits(1),limits(2),limits(1),limits(2)]);      
      hold on
      xlabel ('populaçao no instante "t" ');
      ylabel ('populaçao no instante "t+1"');
      title (['next_p=',next_p]);
      %
      for i=1:50;                             % loop to build cobweb sections
         w=p;                                
         p=eval (next_p);
         plot([w,w],[w,p],'k','EraseMode','none');   % plots vertical cobweb step
         pause(.1);
         if p<0; break; end;                 % eliminates negative population
         plot([w,p],[p,p],'k','EraseMode','none');   % plots horizontal cobweb step
         pause(.1);
      end;
      hold off;
   else continueb=0;                          % stop looping
   end
end
