// TP02.sci

function Nt1 = Hassell(t,Nt0)
    Nt1 = R*Nt0./((1+Nt0).^b);
endfunction

function test01()
    // imgTP01_01.pdf
    R = 3;
    b1 = 1.5;
    b2 = 2;
    b3 = 3.5;
    Nmax = 2;
    Nt0 = 0:0.001:Nmax;
    Nt1 = R*Nt0./(1+Nt0).^b1;
    Nt2 = R*Nt0./(1+Nt0).^b2;
    Nt3 = R*Nt0./(1+Nt0).^b3;
    
    clf()
    
    plot(Nt0,Nt0,'k-',Nt0,Nt1,'-r','thickness',3)
    plot(Nt0,Nt2,'g-','thickness',3)
    plot(Nt0,Nt3,'b-','thickness',3)
    axes = gca();
    axes.font_size = 5;
    axes.x_ticks = ...
        tlist(["ticks","locations","labels"],(0:0.5:Nmax)',(string(0:0.5:Nmax))');
    axes.y_ticks = ...
        tlist(["ticks","locations","labels"],(0:0.5:Nmax)',(string(0:0.5:Nmax))');
    axes.x_label.text = "R=3";
    axes.x_label.font_size = 5;
    legend(gca(),['b=1.5','b=2','b=3.5'],[0.2,1.8]);
    
    xs2png(gcf(),'../Images/imgTP02_01.png');
    
endfunction
    
function test02()
    // imgTP01_02.pdf
    R = 30;
    b1 = 1.5;
    b2 = 2;
    b3 = 3.5;
    Nmax = 10;
    Nt0 = 0:0.001:Nmax;
    Nt1 = R*Nt0./(1+Nt0).^b1;
    Nt2 = R*Nt0./(1+Nt0).^b2;
    Nt3 = R*Nt0./(1+Nt0).^b3;
    
    clf()
    
    plot(Nt0,Nt0,'-r', Nt0,Nt1, '-k','thickness',3)
    plot(Nt0,Nt2,'g-','thickness',3)
    plot(Nt0,Nt3,'b-','thickness',3)
    axes = gca();
    axes.font_size = 5;
    axes.x_ticks = ..
        tlist(["ticks","locations","labels"],(0:1:Nmax)',(string(0:1:Nmax))');
    axes.y_ticks = ..
        tlist(["ticks","locations","labels"],(0:2:Nmax+2)',(string(0:2:Nmax+2))');
    axes.x_label.text = "R=30";
    axes.x_label.font_size = 5;
    legend(gca(),['b=1.5','b=2','b=3.5'],[6.5,11]);
    
    xs2png(gcf(),'../Images/imgTP02_02.png');
    
endfunction

function test03()
    // imgTP01_03.pdf
    bmax = 10;
    b = 2.075:0.01:bmax;
    r = (1-2.0./b).^(-b);
    clf()
    plot(b,r,'-k','thickness',3);
    axes = gca();
    axes.font_size = 5;
    axes.x_label.text = 'b';
    axes.x_label.font_size = 5;
    axes.y_label.text = 'R';
    axes.y_label.font_size = 5;
    axes.data_bounds = [0,0;bmax,800];
    xstring(0.5,200,'stable')
    
    xs2pdf(gcf(),'../Images/imgTP01_03.pdf');
    xs2eps(gcf(),'../Images/imgTP01_03.eps');
endfunction

function test04()
    // imgTP01_04.pdf
    b = 3.5;
    R = 30;    
    Nmax = 5;
    Nt0 = 0:0.001:Nmax;
    Nt1 = R*Nt0./(1+Nt0).^b;
    Nt2 = R*Nt1./(1+Nt1).^b;
    Nstar = R^(1/b)-1;
    
    clf();
    plot(Nt0, Nt0, '-b', Nt0, Nt1, '-k', Nt0, Nt2, '-k', 'thickness', 4)
    plot(Nstar,Nstar,'.r', "MarkSize", 8)
    
    axes =gca();
    axes.font_size = 6;
    axes.x_ticks = ..
        tlist(["ticks","locations","labels"],(0:1:Nmax)',(string(0:1:Nmax))');
    axes.y_ticks = ..
        tlist(["ticks","locations","labels"],(0:1:Nmax)',(string(0:1:Nmax))');
    xstring(2,-0.1,'$\boldsymbol{N_t}$');
    xstring(0.1,4.3,"$\boldsymbol{N}$");
    xstring(2.7,0.5,'$\boldsymbol{\leftarrow N=N_{t+1}}$');
    xstring(0.8,3,"$\boldsymbol{N=N_{t+2}\rightarrow}$");
    xstring(Nstar-0.1,0.9,"$\boldsymbol{N_*}$");
    
    xs2pdf(gcf(),'../Images/imgTP02_04.pdf');
    xs2png(gcf(),'../Images/imgTP02_04.png');
endfunction

function test05()
    // imgTP01_05.pdf
    Nmax = 5;
    NNt0 = 0.1:0.01:Nmax;
    bb = (2.6:0.01:8)';
    lNN = length(NNt0);
    lbb = length(bb);
    Nt0 = repmat(NNt0,lbb,1);
    b = repmat(bb,1,lNN);
    R = 30;
    Nt1 = ode("discrete",Nt0,1,2000,Hassell);
    Nt1 = gsort(Nt1,'c','i');
    clf()
    stacksize('max');
    plot(b',Nt1',"Linestyle","none", "MarkStyle", ".", "MarkSize", 1)
    xlabel("$b$")
    ylabel("$N_*$")
    axes =gca();
    axes.font_size = 5;
    axes.x_label.font_size = 5;
    axes.y_label.font_size = 5;
    //xs2eps(gcf(),'../Images/imgTP01_05.eps');
    //xs2gif(gcf(),'../Images/imgTP01_05.gif');
    //xs2pdf(gcf(),'../Images/imgTP01_05.pdf');
    xs2png(gcf(),'../Images/imgTP02_05.png')
endfunction

//test01()
//test02()
//test03()
test04()
//test05()
