//TP04.sci

function y= fauxLogistique(t,x)
    y = x.*(1 + x) - 4.*x./(1+x);
endfunction

function y = fourche(t,x)
    // N = nombre de données initales x;
    // L = nombre de paramètres r
    // x,r = matrices de dimensions N*L ;
    y = r.*x.*(1-x.^2/2)-5*(x-x.^3/6);
endfunction

function y = pendule(t,x)
    // x = matrice de taille N*2 ; N nombre de conditions initiales
    y = [x(:,2),-sin(x(:,1) - gam*x(:,2))];
endfunction

function y = modelePeche(t,x)
    y = r.*x.*(1-x./q) - (1-exp(-x.^4));
endfunction

function y = modelePecheBis(t,x)
    r = r0*cos(2*%pi*t/T) + r1;
    y = r.*x.*(1-x./q) - (1-exp(-x.^4));
endfunction

function test01()
    clf()
    x = 0:0.5:2;
    t = 0:0.5:5;
    t_prime = ones(1,length(t));
    x_prime = fauxLogistique(t,x);
    [T_PRIME,X_PRIME] = meshgrid(t_prime,x_prime);
    champ(t,x,T_PRIME',X_PRIME',arfact=1.7)
    
    axes = gca();
    axes.thickness = 3;
    axes.font_size = 5.5;
    xstring(4.5,-0.4,"$t$");
    xstring(0.1,2.5,"$x(t)$");
    axes.data_bounds = [0,-0.5;5,3];
    
    t0 = 0;
    tt = 0:0.01:5;
    x0 = [0;1];
    xt = ode("adams",x0,t0,tt,fauxLogistique);
    plot(tt',xt','-b')
    x0 = [0.9;0.995;0.9996];
    xt = ode("adams",x0,t0,tt,fauxLogistique);
    plot(tt',xt','-r')
    
    t0 = 0;
    x0 = 1.1;
    tt = 0:0.01:1.04;
    xt = ode("adams",x0,t0,tt,fauxLogistique);
    plot(tt',xt','-r')
    x0 = 1.001;
    tt = 0:0.01:3.3;
    xt = ode("adams",x0,t0,tt,fauxLogistique);
    plot(tt',xt','-r')
    x0 = 1.0001;
    tt = 0:0.01:4.45;
    xt = ode("adams",x0,t0,tt,fauxLogistique);
    plot(tt',xt','-r')
  
    xs2pdf(gcf(),'../Images/imgTP04_01.pdf')
    xs2gif(gcf(),'../Images/imgTP04_01.gif')
endfunction

function test02()
    xmin = -3; xmax = 3;
    ymin = -12; ymax = 12;
    x = xmin:0.01:xmax;
    clf()
    for r = [1,5,15]
        y =  fourche(0,x);  
        plot(x,y,'-b','LineWidth',3)
    end
    plot([xmin,xmax],[0,0],'-k')
    plot(0,0,'.r',"MarkSize",10)
    axes = gca();
    axes.thickness = 3;
    axes.font_size = 5;
    xstring(xmax-0.5,ymin+1,"$\text{x}$");
    xstring(xmin,ymax-3,"$x''$")
    xstring(2.1,-4,"$r=1$")
    xstring(-2.5,7,"$r=5$")
    xstring(-2,-5,"$r=15$")
    xstring(0,-3,"$x_*$")
    legend("$x'' = f_r(x)$");
    axes.data_bounds = [xmin,ymin;xmax,ymax];
    
    xs2pdf(gcf(),'../Images/imgTP04_02.pdf')
    xs2gif(gcf(),'../Images/imgTP04_02.gif')
endfunction

function test03()
    clf();
    x0 = [-0.1,0.1]';
    r0 = 4:0.001:6;
    N = length(x0);
    L = length(r0);
    XX0 = repmat(x0,1,L);
    r = repmat(r0,N,1);
    t0 = 0;
    TT = 500;
    stacksize('max');
    XX = ode(XX0,t0,TT,fourche);
    plot(r0',XX',"Linestyle","none", "MarkStyle", ".", "MarkSize", 1)
    axes = gca();
    axes.thickness = 3;
    axes.font_size = 5;
    xstring(5.8,-0.75,"$r$")
    xstring(4.1,0.6,"$x_*(r)$")
    
    xs2pdf(gcf(),'../Images/imgTP04_03.pdf')
    xs2gif(gcf(),'../Images/imgTP04_03.gif')
    xs2png(gcf(),'../Images/imgTP04_03.png')
endfunction

function test04()
    gam = 0;
    x0 = [ 1,1; ...
           0.5,0.5; ...
           0.25,0.25; ...
           1+2*%pi,1; ...
           0.5+2*%pi,0.5; ...
           0.25+2*%pi,0.25; ...
           1-2*%pi,1; ...
           0.5-2*%pi,0.5; ...
           0.25-2*%pi,0.25; ...
           -3*%pi,0.01; ...
           3*%pi,-0.01; ...
           -10,2; ...
           10,-2; ...
           -10,3.5; ...
           10,-3.5 ];
    t0 = 0;
    TT = 0:0.1:40;
    L = length(TT);
    xx = ode(x0,t0,TT,pendule);
    clf()
    plot(xx(:,1:2:(2*L-1))',xx(:,2:2:(2*L))','thickness',3)
    xlabel("x")
    ylabel("y")
    axes = gca();
    axes.isoview = "on";
    axes.data_bounds = [-10,-5;10,5];
    axes.margins = [0.08, 0.08, 0.01, 0.01];
    
    xs2pdf(gcf(),'../Images/imgTP04_04.pdf')
    xs2gif(gcf(),'../Images/imgTP04_04.gif')
    xs2png(gcf(),'../Images/imgTP04_04.png')
    xs2eps(gcf(),'../Images/imgTP04_04.eps')
endfunction

function test05()
    gam = 0.5;
    x0 = [ 1,1; ...
           0.5,0.5; ...
           0.25,0.25; ...
           1+2*%pi,1; ...
           0.5+2*%pi,0.5; ...
           0.25+2*%pi,0.25; ...
           1-2*%pi,1; ...
           0.5-2*%pi,0.5; ...
           0.25-2*%pi,0.25; ...
           -3*%pi,0.01; ...
           3*%pi,-0.01; ...
           -10,2; ...
           10,-2; ...
           -10,3.5; ...
           10,-3.5 ];
    t0 = 0;
    TT = 0:0.1:100;
    L = length(TT);
    xx = ode(x0,t0,TT,pendule);
    clf()
    plot(xx(:,1:2:(2*L-1))',xx(:,2:2:(2*L))','thickness',3)
    axes = gca();
    axes.isoview = "on";
    axes.thickness = 2;
    axes.font_size = 4;
    axes.x_label.text = "x";
    axes.x_label.font_size = 5;
    axes.y_label.text = "y";
    axes.y_label.font_size = 5;
    axes.data_bounds = [-10,-5;10,5];
    axes.margins = [0.08, 0.08, 0.01, 0.01];
    
    xs2pdf(gcf(),'../Images/imgTP04_05.pdf')
    xs2gif(gcf(),'../Images/imgTP04_05.gif')
    xs2png(gcf(),'../Images/imgTP04_05.png')
endfunction

//function test06()
//-- Vous avez s\^urement constat\'e qu'on ne peut pas obtenir de comportements p\'eriodiques avec des \'equations diff\'erentielles du premier ordre. On consid\`ere ici l'\'equation du pendule amorti
//\[
//\theta'' + \gamma \theta' + k\sin\theta = 0, \quad\textrm{(on prendra $k=1$.)}.
//\]
//\noindent C'est une \'equation du second ordre que l'on r\'ecrit sous la fome d'un syst\`eme
//\[
//\begin{cases}
//x' = y, \\
//y' = -\sin x - \gamma y.
//\end{cases}
//\]
//\noindent On demande de tracer quelques trajectoires dans l'espace des phases $(x,y)$,  lorsque le frottement $\gamma$ est nul et lorsqu'il est positif. On montrera que $x=n\pi$ et $y=0$ sont les points stationnaires. On constatera alors l'existence de solutions p\'eriodiques proches des points stables $x=2n\pi$, $y=0$ dans le cas non amorti. On devra trouver la figure 
//\begin{figure}[h]
//\centering
//\vspace{-0.5cm}
//\includegraphics[width=\textwidth]{../Images/imgTP04_04} 
//\vspace{-1cm}
//\caption{Trajectoire du pendule non amorti}
//\label{imgTP04_04}
//\end{figure}
//\begin{figure}[h]
//\centering
//\includegraphics[width=\textwidth]{../Images/imgTP04_05} 
//\vspace{-1cm}
//\caption{Trajectoire du pendule amorti}
//\label{imgTP04_05}
//\end{figure}
//endfunction

function test07()
    clf();
    xmin = 0; xmax = 8;
    ymin = -1; ymax = 1;
    xx = xmin:0.01:xmax;
    tt = 0:0.01:10;
    q = 8; r = 0.3; plot(xx,modelePeche(tt,xx),'-b','LineWidth',3)
    q = 8; r = 0.5; plot(xx,modelePeche(tt,xx),'-b','LineWidth',3)
    q = 8; r = 0.87; plot(xx,modelePeche(tt,xx),'-b','LineWidth',3)
    q = 8; r = 1; plot(xx,modelePeche(tt,xx),'-b','LineWidth',3)
    
    plot(6.85,0,'.r','MarkSize',10)
    plot(0.65,0,'.r','MarkSize',10)
    
    plot([xmin,xmax],[0,0],'-k','LineWidth',3)
    
    a = gca();
    a.thickness = 3;
    a.font_size = 5.5;
    a.data_bounds = [xmin,ymin;xmax,ymax]
    xstring(7,-0.95,"$x$");
    xstring(0.1,0.7,"$x''$");
    xstring(3.2,-0.6,"$r=0.3$");
    xstring(3.2,-0.2,"$r=0.5$");
    xstring(3.2,0.45,"$r=0.87$");
    xstring(3.2,0.8,"$r=1$");
    xstring(0.1,-0.25,"$x_*^-$");
    xstring(7,0,"$x_*^+$");
    
    xs2pdf(gcf(),'../Images/imgTP04_06.pdf')
endfunction

function test08()
    clf();
    qq = 0:0.1:8;
    tt = 0:0.01:10;
    xmin = 0; xmax = 8;
    rmin = 0; rmax = 1;
    zmin = -1; zmax = 1;
    xx = xmin:0.1:xmax;
    rr = rmin:0.1:rmax;
    [r,x] = meshgrid(rr,xx);
    q = 8;
    x_prime = modelePeche(tt,x);
    f=gcf(); f.color_map = graycolormap(128);
    //f=gcf(); f.color_map = bonecolormap(128);
    surf(x,r,x_prime,'FaceColor','interp');
    e = gce();
    e.color_mode = 1;
    a = gca();
    a.font_size = 4;
    a.data_bounds = [xmin,rmin,zmin;xmax,rmax,zmax];
    xlabel("$x$","fontsize",6);
    ylabel("$r$","fontsize",6);
    zlabel("$x''$","fontsize",6);
    
    xs2pdf(gcf(),'../Images/imgTP04_07.pdf')
    xs2png(gcf(),'../Images/imgTP04_07.png')
endfunction

function test09()
    clf()
    q = 8;
    r0 = 0.3; r1 = 0.7;
    x0 = 4; t0 = 0;
    tt = 0:0.01:100;
    
    
    subplot(2,3,1)
    T = 3;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,4,"$T=3$")
    subplot(2,3,2)
    T =5;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,4,"$T=5$")
    subplot(2,3,3)
    T = 10;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,4,"$T=10$")
    subplot(2,3,4)
    T = 20;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,4,"$T=20$")
    subplot(2,3,5)
    T = 30;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,4,"$T=30$")
    subplot(2,3,6)
    T = 50;
    xx = ode("adams",x0,t0,tt,modelePecheBis);
    plot(tt,xx,'LineWidth',3)
    a = gca();
    a.font_size = 4;
    a.thickness = 3;
    xstring(30,3,"$T=50$")
    
    xs2pdf(gcf(),'../Images/imgTP04_08.pdf')
endfunction

//test01()
//test02()
test03()
//test04()
//test05()
//test07()
//test08()
//test09()