// TP05.sci

exec("dde.sci");

function xdot = cosRetardee(t,x)
    z = ddeGetBackwardPositions(t,x,t-T);
    xdot = -(%pi/2*T)*z;
endfunction

function  y = initCosRetardee(bt)
    y = cos(2*%pi*(bt-T)/(4*T)); 
endfunction

function  y = initCosRetardeeBis(bt)
    //y = cos(2.99999999999*%pi*bt/(2*T)); 
    //y = cos(3.0000000000000001*%pi*bt/(2*T));
    y = cos(2*%pi*bt/(T/5));
endfunction

function xdot = logistiqueRetardee(t,x)
    z = ddeGetBackwardPositions(t,x,t-T);
    xdot = x.*(1-z);
endfunction

function y = initLogistiqueRetardee(bt)
    y = x_0*ones(1,length(bt));
endfunction

function N_plus = NicholsonBailey(t,N)
    //N = matrice de dimension L*2
    //L nombre de conditions initiales
    //N(:,1) = composante x_t ; N(:,2) = composante y_t
    X_plus = N(:,1).*exp(R*(1-N(:,1))-N(:,2));
    Y_plus = GAMMA*N(:,1).*(1-exp(-N(:,2)));
    N_plus = [X_plus,Y_plus];
endfunction

function [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile)
    //calcule les paramètres manquant de Nicholson-Bailey
    //R, x_etoile, y_etoile, GAMMA sont des matrices de même dimension
    y_etoile = R.*(1-x_etoile);
    GAMMA = y_etoile./(x_etoile.*(1-exp(-y_etoile)));
endfunction

function test01() 
    T = 1;
    t_init = 0;
    t_final = 8*T;
    t = t_init:0.01:t_final;
    
    x = dde(T,initCosRetardee,t,cosRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,-1.2;8,1.2];
    xstring(7.5,-1.1,"$\Large t$")
    xstring(0.3,0.9,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_01.pdf')
endfunction

function test02() 
    T = 1;
    t_init = 0;
    t_final = 8*T;
    t = t_init:0.01:t_final;
    
    x = dde(T,initCosRetardeeBis,t,cosRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,-1.5;8,1.5];
    xstring(7.5,-1.4,"$\Large t$")
    xstring(0.3,1.1,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_02.pdf')
endfunction

function test03()
    T = %pi/2 - 0.01;
    t_init = 0;
    t_final = 100*T;
    t = t_init:0.01:t_final;
    x_0 = 0.1;
    
    x = dde(T,initLogistiqueRetardee,t,logistiqueRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,0;150,2];
    xstring(140,0.05,"$\Large t$")
    xstring(8,1.7,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_03.pdf')
endfunction

function test04()
    T = %pi/2 - 0.1;
    t_init = 0;
    t_final = 100*T;
    t = t_init:0.01:t_final;
    x_0 = 0.1;
    
    x = dde(T,initLogistiqueRetardee,t,logistiqueRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,0;150,2];
    xstring(140,0.05,"$\Large t$")
    xstring(8,1.7,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_04.pdf')
endfunction

function test05()
    T = %pi/2 + 0.15;
    t_init = 0;
    t_final = 100;
    t = t_init:0.01:t_final;
    x_0 = 0.1;
    
    x = dde(T,initLogistiqueRetardee,t,logistiqueRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,0;100,2.6];
    xstring(95,0.05,"$\Large t$")
    xstring(2,2.2,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_05.pdf')
endfunction

function test06()
    T = %pi/2 + 0.5;
    t_init = 0;
    t_final = 100;
    t = t_init:0.01:t_final;
    x_0 =0.1;
    
    x = dde(T,initLogistiqueRetardee,t,logistiqueRetardee);
    clf()
    plot(t,x,'-b','LineWidth',3)
    a = gca();
    a.tight_limits = "on";
    a.font_size = 5;
    a.thickness = 3;
    a.data_bounds = [0,0;100,3.4];
    xstring(95,0.05,"$\Large t$")
    xstring(2,3,"$\Large x(t)$")
    
    xs2pdf(gcf(),'../Images/imgTP05_06.pdf')
endfunction

function test08()
    clear
    clf
    
    x_etoile = 0.5; R = 1.5;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:30;
    //quadrillage
    lx0 = 10; ly0 = 50;
    x0 = linspace(0.1,1.2,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    plot(Nf(:,1:2:($-1))',Nf(:,2:2:$)','.k',"MarkerSize",1.4)
    x_min = 0; x_max = 1.2;
    y_min = 0; y_max = R;
    a = gca()
    a.data_bounds = [x_min,y_min;x_max,y_max];
    a.font_size = 5;
    xstring(1.1,0,"$\Large x_t$");
    xstring(0.03,0.95*R,"$\Large y_t$");
    xstring(0.5,0.95*R,["$\Large x_*=$",string(x_etoile)," ",...
                         "$\Large y_*=$",string(y_etoile)]);
    xstring(0.535,0.85*R,["$\Large r=$",string(R),"$\quad\,$",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);
    
    xs2pdf(gcf(),'../Images/imgTP05_08.pdf')
    xs2png(gcf(),'../Images/imgTP05_08.png')
endfunction

function test09()
    clear
    clf
    
    x_etoile = 0.5; R = 3;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:30;
    //quadrillage
    lx0 = 10; ly0 = 50;
    x0 = linspace(0.1,1.2,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    plot(Nf(:,1:2:($-1))',Nf(:,2:2:$)','.k',"MarkerSize",1.4)
    plot(x_etoile,y_etoile,'.r',"Markersize",10)
    x_min = 0; x_max = 1.2;
    y_min = 0; y_max = R;
    a = gca()
    a.data_bounds = [x_min,y_min;x_max,y_max];
    a.font_size = 5;
    xstring(1.1,0,"$\Large x_t$");
    xstring(0.03,0.87*R,"$\Large y_t$");
    xstring(0.5,0.87*R,["$\Large x_*=$",string(x_etoile)," ",...
                         "$\Large y_*=$",string(y_etoile)]);
    xstring(0.535,0.77*R,["$\Large r=$",string(R),"$\quad\quad\ $",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);
    
    xs2pdf(gcf(),'../Images/imgTP05_09.pdf')
    xs2png(gcf(),'../Images/imgTP05_09.png')
endfunction

function test10()
    clear
    clf
    
    x_etoile = 0.5; R = 2.2;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:500;
    //quadrillage
    lx0 = 10; ly0 = 10;
    x0 = linspace(0.1,1,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    plot(Nf(:,451:2:($-1))',Nf(:,452:2:$)','.k',"MarkerSize",1.4)
    plot(x_etoile,y_etoile,'.r',"Markersize",10)
    x_min = 0.2; x_max = 0.8;
    y_min = 0.4; y_max = 1.8;
    a = gca()
    a.data_bounds = [x_min,y_min;x_max,y_max];
    a.font_size = 5;
    xstring(0.75,0.4,"$\Large x_t$");
    xstring(0.21,1.65,"$\Large y_t$");
    xstring(0.4,1.65,["$\Large x_*=$",string(x_etoile)," ",...
                         "$\Large y_*=$",string(y_etoile)]);
    xstring(0.42,1.5,["$\Large r=$",string(R),"$\quad\,$",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);
    
    xs2pdf(gcf(),'../Images/imgTP05_10.pdf')
    xs2png(gcf(),'../Images/imgTP05_10.png')
endfunction

function test11()
    clear
    clf
    stacksize('max') 
    
    x_etoile = 0.5; R = 3.3;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:3000;
    //quadrillage
    lx0 = 20; ly0 = 20;
    x0 = linspace(0.1,1.2,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    Nf = Nf(:,5901:$);
    plot(Nf(:,1:2:($-1))',Nf(:,2:2:$)','.k',"MarkerSize",1.4)
    plot(x_etoile,y_etoile,'.r',"MarkerSize",10)
    x_min = 0.1; x_max = 1.1;
    y_min = 0.2; y_max = R;
    a = gca()
    a.tight_limits = "on";
    a.data_bounds = [x_min,y_min;x_max,y_max];
    //a.data_bounds = [0.14,1.04;0.2,1.15]
    a.font_size = 5;
    xstring(1,0.3,"$\Large x_t$");
    xstring(0.12,3,"$\Large y_t$");
    xstring(0.4,3,["$\Large x_*=$",string(x_etoile)," ",...
                         "$\Large y_*=$",string(y_etoile)]);
    xstring(0.42,2.7,["$\,\Large r=$",string(R),"$\quad\,$",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);

    
    xs2pdf(gcf(),'../Images/imgTP05_11.pdf')
    xs2png(gcf(),'../Images/imgTP05_11.png')
endfunction

function test12()
    clear
    clf
    stacksize('max') 
    
    x_etoile = 0.5; R = 2.95;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:100;
    //quadrillage
    lx0 = 20; ly0 = 20;
    x0 = linspace(0.1,1.2,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    Nf = Nf(:,51:$);
    plot(Nf(:,1:2:($-1))',Nf(:,2:2:$)','.k',"MarkerSize",1.4)
    plot(x_etoile,y_etoile,'.r',"MarkerSize",10)
    x_min = 0.1; x_max = 1.2;
    y_min = 0.2; y_max = R;
    a = gca()
    a.tight_limits = "on";
    a.data_bounds = [x_min,y_min;x_max,y_max];
    //a.data_bounds = [0.14,1.04;0.2,1.15]
    a.font_size = 5;
    xstring(1.1,0.2,"$\Large x_t$");
    xstring(0.12,2.7,"$\Large y_t$");
    xstring(0.4,2.7,["$\Large x_*=$",string(x_etoile),"$\quad$",...
                         "$\Large y_*=$",string(y_etoile)]);
    xstring(0.42,2.4,["$\ \Large r=$",string(R),"$\quad$",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);

    
    xs2pdf(gcf(),'../Images/imgTP05_12.pdf')
    xs2png(gcf(),'../Images/imgTP05_12.png')
endfunction

function test13()
    clear
    clf
    stacksize('max') 
    
    x_etoile = 0.5; R = 2.95;
    [y_etoile,GAMMA] = paramNicholsonBailey(R,x_etoile);
    t0 = 0;
    tf = 1:100;
    //quadrillage
    lx0 = 40; ly0 = 40;
    x0 = linspace(0.1,1.2,lx0);
    y0 = linspace(0.1,R,ly0);
    X0 = repmat(x0',1,ly0);
    Y0 = repmat(y0,lx0,1).*ones(lx0,ly0);
    N0 = [X0(:),Y0(:)];
    Nf = ode("discrete",N0,t0,tf,NicholsonBailey);
    Nf = Nf(:,51:$);
    plot(Nf(:,1:2:($-1))',Nf(:,2:2:$)','.k',"MarkerSize",1.4)
    //plot(x_etoile,y_etoile,'.r',"MarkerSize",10)
    x_min = 0.1; x_max = 1.2;
    y_min = 0.2; y_max = R;
    a = gca()
    a.tight_limits = "on";
    //a.data_bounds = [x_min,y_min;x_max,y_max];
    a.data_bounds = [0.16,1;0.19,1.16]
    a.font_size = 5;
    xstring(0.188,1.003,"$\Large x_t$");
    xstring(0.16,1.143,"$\Large y_t$");
    xstring(0.17,1.143,["$\ \Large r=$",string(R),"$\quad\quad\quad$",...
                         "$\Large \gamma=$",msprintf("%2.2f",GAMMA)]);

    
    xs2pdf(gcf(),'../Images/imgTP05_13.pdf')
    xs2png(gcf(),'../Images/imgTP05_13.png')
endfunction


//test01()
//test02()
//test03()
//test04()
//test05()
//test06()
//test08()
//test09()
//test10()
//test11()
//test12()
//test13()