function TD_Phase_6(action)
%LORENZ Plot the orbit around the Lorenz chaotic attractor.
%   This demo animates the integration of  the    
%   three coupled nonlinear differential equations
%   that define the "Lorenz Attractor", a chaotic 
%   system first described by Edward Lorenz of    
%   the Massachusetts Institute of Technology.    
%
%   As the integration proceeds you will see a    
%   point moving around in a curious orbit in     
%   3-D space known as a strange attractor. The   
%   orbit is bounded, but not periodic and not    
%   convergent (hence the word "strange").        
%
%   Use the "Start" and "Stop" buttons to control 
%   the animation.                                

%   Adapted for Demo by Ned Gulley, 6-21-93; jae Roh, 10-15-96
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.9 $  $Date: 1997/11/21 23:26:25 $

% Possible actions:
% initialize
% close
global RHO

% Information regarding the play status will be held in
% the axis user data according to the following table:
play= 1;
stop=-1;

global RHO

if nargin<1,
    action='initialize';
end;

if strcmp(action,'initialize'),
    
    RHO = 24;
    oldFigNumber=watchon;

    figNumber=figure( ...
        'Name','Lorenz Attractor', ...
        'NumberTitle','off', ...
        'Visible','off');
    colordef(figNumber,'black')
    axes( ...
        'Units','normalized', ...
        'Position',[0.05 0.10 0.75 0.95], ...
        'Visible','off');

    text(0,0,'Press the "Start" button to see the Lorenz demo', ...
        'HorizontalAlignment','center');
    axis([-1 1 -1 1]);

    %===================================
    % Information for all buttons
    labelColor=[0.8 0.8 0.8];
    yInitPos=0.90;
    xPos=0.75;
    btnLen=0.20;
    btnWid=0.10;
    % Spacing between the button and the next command's label
    spacing=0.05;

    %====================================
    % The CONSOLE frame
    frmBorder=0.02;
    yPos=0.05-frmBorder;
    frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];
    h=uicontrol( ...
        'Style','frame', ...
        'Units','normalized', ...
        'Position',frmPos, ...
        'BackgroundColor',[0.50 0.50 0.50]);
    
    %====================================
    % The START button
    btnNumber=1;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Start';
    cmdStr='start';
    callbackStr='TD_Phase_6(''start'');';

    % Generic button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    startHndl=uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Interruptible','on', ...
        'Callback',callbackStr);
    %====================================
    % The STOP button
    btnNumber=2;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Stop';
    % Setting userdata to -1 (=stop) will stop the demo.
    callbackStr='set(gca,''Userdata'',-1)';
    
    % Generic  button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    stopHndl=uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'Enable','off', ...
        'String',labelStr, ...
        'Callback',callbackStr);
    

    
%====================================
    % RHo slider
    btnNumber=3;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing)
    labelStr='RHO';
    cmdStr='Rho';
    callbackStr='TD_Phase_6(''Rho'');';
    % Generic  button information
    btnPos=[xPos yPos-spacing btnLen btnWid/2];
    RhoHndl =uicontrol('Style','slider','Min',20,'Max',28,'SliderStep',[.005 .05],...
         'Units','normalized', ...
        'Position',btnPos, ...
        'Enable','on', ...
        'String',labelStr, ...
        'Value',24, ...
        'Callback',callbackStr);
     
    % RHo affiche
    btnNumber=4;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing)
    
    % Generic  button information
    btnPos=[xPos yPos-spacing btnLen btnWid/2];
    AffHndl =uicontrol('Style','text',...
         'Units','normalized', ...
        'Position',btnPos, ...
        'Enable','on', ...
        'String',RHO);
    



%====================================
    % The INFO button
    labelStr='Info';
    callbackStr='TD_Phase_6(''info'')';
    infoHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.20 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);

    %====================================
    % The CLOSE button
    labelStr='Close';
    callbackStr='close(gcf)';
    closeHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.05 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);
    

     % Uncover the figure
    hndlList=[startHndl stopHndl infoHndl closeHndl RhoHndl AffHndl];
    set(figNumber,'Visible','on', ...
        'UserData',hndlList);
     
    watchoff(oldFigNumber);
    figure(figNumber);
    
elseif strcmp(action,'Rho'),
    axHndl=gca;
    figNumber=gcf;
    hndlList=get(figNumber,'UserData');
    startHndl=hndlList(1);
    stopHndl=hndlList(2);
    infoHndl=hndlList(3);
    closeHndl=hndlList(4);
    RhoHndl=hndlList(5);
    AffHndl=hndlList(6);
    
    RHO = get(RhoHndl,'Value')
    set(AffHndl,'String',RHO)
    
    
elseif strcmp(action,'start'),
    axHndl=gca;
    figNumber=gcf;
    hndlList=get(figNumber,'UserData');
    startHndl=hndlList(1);
    stopHndl=hndlList(2);
    infoHndl=hndlList(3);
    closeHndl=hndlList(4);
    RhoHndl=hndlList(5);
    
    set([startHndl closeHndl infoHndl],'Enable','off');
    set([stopHndl RhoHndl],'Enable','on');

    % ====== Start of Demo
    set(figNumber,'Backingstore','off');
    % The graphics axis limits are set to values known 
    % to contain the solution.
    set(axHndl, ...
        'XLim',[-10 20],'YLim',[-50 20],'ZLim',[-10 50], ... 
        'Userdata',play, ...
        'XTick',[],'YTick',[],'ZTick',[], ...
        'Drawmode','fast', ...
        'Visible','on', ...
        'NextPlot','add', ...
        'Userdata',play, ...
        'View',[-25,10]);  % -37.5 30
    xlabel('X');
    ylabel('Y');
    zlabel('Z');

    % The values of the global parameters are
    global SIGMA RHO BETA
    SIGMA = 10.;
    %RHO = input('Rho = ');
    BETA = 8./3.;
 
    % The orbit ranges chaotically back and forth around two different points,
    % or attractors.  It is bounded, but not periodic and not convergent.
    % The numerical integration, and the display of the evolving solution,
    % are handled by the function ODE23P.
 
    FunFcn='lorenzeq';
    Temps = [];
    Coord=[];
    % The initial conditions below will produce good results
    %y0 = [20 5 -5];
    % Random initial conditions
    y0(1)=rand*30+5;
    y0(2)=rand*35-30;
    y0(3)=rand*40-5;
    t0=0;
    tfinal=150;
    pow = 1/3;
    tol = 0.0001;
 
    t = t0;
    hmax = (tfinal - t)/100;
    hmin = (tfinal - t)/200000;
    h = (tfinal - t)/2500;
    y = y0(:);
    tau = tol * max(norm(y,'inf'),1);
 
    % Save L steps and plot like a comet tail.
    L = 200;
    Y = y*ones(1,L);
 
    cla;
    head = line( ...
        'color','r', ...
        'Marker','.', ...
        'markersize',25, ...
        'erase','xor', ...
        'xdata',y(1),'ydata',y(2),'zdata',y(3));
    body = line( ...
        'color','y', ...
        'LineStyle','-', ...
        'erase','none', ...
        'xdata',[],'ydata',[],'zdata',[]);
    tail=line( ...
        'color','g', ...
        'LineStyle','-', ...
        'erase','none', ...
        'xdata',[],'ydata',[],'zdata',[]);
    afac=line( ...
        'color','k', ...
        'LineStyle','-', ...
        'erase','none', ...
        'visible','on',...
        'xdata',[],'ydata',[],'zdata',[]);
 
    % The main loop
    while (get(axHndl,'Userdata')==play) & (h >= hmin)
        if t + h > tfinal, h = tfinal - t; end
        % Compute the slopes
        s1 = feval(FunFcn, t, y);
        s2 = feval(FunFcn, t+h, y+h*s1);
        s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4);
 
        % Estimate the error and the acceptable error
        delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
        tau = tol*max(norm(y,'inf'),1.0);
 
        % Update the solution only if the error is acceptable
        ts = t;
        ys = y;
        if delta <= tau
            t = t + h;
            y = y + h*(s1 + 4*s3 + s2)/6;
            %Temps = [Temps, t];
            %Coord = [Coord, y];
            %save Loren Temps Coord

            % Update the plot
            Y = [y Y(:,1:L-1)];
            set(head,'xdata',Y(1,1),'ydata',Y(2,1),'zdata',Y(3,1))
            set(body,'xdata',Y(1,1:12),'ydata',Y(2,1:12),'zdata',Y(3,1:12))
            set(tail,'xdata',Y(1,13:L-2),'ydata',Y(2,13:L-2),'zdata',Y(3,13:L-2))
            set(afac,'xdata',Y(1,L-2:L),'ydata',Y(2,L-2:L),'zdata',Y(3,L-2:L))
            
            drawnow;
        end
 
        % Update the step size
        if delta ~= 0.0
            h = min(hmax, 0.9*h*(tau/delta)^pow);
        end

    end;    % Main loop ...
    % ====== End of Demo
    set([startHndl closeHndl infoHndl RhoHndl],'Enable','on');
    set(stopHndl,'Enable','off');

elseif strcmp(action,'info');
    helpwin(mfilename);

end;    % if strcmp(action, ...

function ydot = lorenzeq(t,y)
%LORENZEQ Equation of the Lorenz chaotic attractor.
%   ydot = lorenzeq(t,y).
%   The differential equation is written in almost linear form.


global SIGMA RHO BETA

%A = [ -BETA    0     y(2)  
%         0  -SIGMA   SIGMA 
%      -y(2)   RHO    -1  ];
%RHO
A = [ -SIGMA    SIGMA     0  
         RHO   -1  -y(1) 
      y(2)   0    -8/3  ];

ydot = A*y;