LorentzAttractor

classdef LorentzAttractor < handle
    %LORENTZATTRACTOR : simulation of the Lorentz Attractor chaotic system
    %   The Lorentz attractor describes a model of the atmosphere that can
    %   become chaotic in certain parameter regimes.  This class
    %   encapsulates the simulation parameters and has methods to run the
    %   simulation and plot the results.
    
    %   Copyright 2008-2009 The MathWorks, Inc.
    %   $Revision: 35 $  $Date: 2009-05-29 15:27:34 +0100 (Fri, 29 May 2009) $
    
    properties
        Sigma
        Beta
        Rho
        MaxTime
        SampleCount
        StateValues
        TimeValues
        X0
    end
    
    methods
        function obj = LorentzAttractor(varargin)
            % Parse the input arguments using an input parser object
            p = inputParser;
            % sigma, beta and rho are the model parameters
            p.addParamValue('sigma',10);
            p.addParamValue('beta',8/3);
            p.addParamValue('rho',28);
            % Maximum simulation time
            p.addParamValue('tmax',25);
            % Number of points to return from the simulation result
            p.addParamValue('npts',1e4);
            % Starting state
            p.addParamValue('x0',[0 1 1.05]');
            % Parse the input arguments
            p.parse(varargin{:});
            % Store the parameters as properties of the object
            obj.Sigma = p.Results.sigma;
            obj.Beta = p.Results.beta;
            obj.Rho = p.Results.rho;
            obj.MaxTime = p.Results.tmax;
            obj.SampleCount = p.Results.npts;
            obj.X0 = p.Results.x0;
        end
        
        function simulate(obj)
            % Perform the simulation
            [obj.TimeValues,obj.StateValues] = ode45(...
                @obj.propagate,[0 obj.MaxTime],obj.X0);
        end
        
        function y = propagate(obj,t,x) %#ok<INUSL>
            % Equations of motion.  Time argument is unused for this system
            % as the equations of motion are time independent, however the
            % ODE solver requires equations of motion of the form 
            % f = f(t,x)
            y = zeros(size(x));
            y(1) = obj.Sigma.*(x(2) - x(1));
            y(2) = x(1).*(obj.Rho - x(3)) - x(2);
            y(3) = x(1).*x(2) - obj.Beta .* x(3);
        end
        
        function plot3(obj,varargin)
            % Plot the system evolution in three dimensions.  If the
            % simulation has not been performed yet then perform it here.
            if isempty(obj.StateValues)
                obj.simulate();
            end
            % Plot the evolution path
            plot3(obj.StateValues(:,1),...
                obj.StateValues(:,2),...
                obj.StateValues(:,3),varargin{:});
            % Show the starting point
            startingPoint = line(obj.X0(1),...
                obj.X0(2),...
                obj.X0(3));
            set(startingPoint,'LineStyle','none','Marker','o',...
                'MarkerSize',8,'MarkerFaceColor','k');
            % Set a nice view angle    
            view(45, 45);
        end
        
        function plot2D(obj,columnIndexes,varargin)
            % Plot a 2D slice of the evolution, simulating the motion if
            % the simulation has not been run already.
            % columnIndexes is a 2 element vector of the dimensions to plot
            % ie x=1, y=2, z=3
            if isempty(obj.StateValues)
                obj.simulate();
            end
            
            % Default to plotting y vs x if columnIndexes is unspecified
            if nargin<2
                columnIndexes = [1 2];
            end
            
            % Do the plot
            plot(obj.X0(columnIndexes(1)),...
                obj.X0(columnIndexes(2)),...
                'k.',...
                obj.StateValues(:,columnIndexes(1)),...
                obj.StateValues(:,columnIndexes(2)),...
                varargin{:});            
        end
        
        function plotXY(obj,varargin)
            % Plot the evolution in the XY plane
            obj.plot2D([1,2],varargin{:}); 
        end
        
        function plotYZ(obj,varargin)
            % Plot the evolution in the YZ plane
            obj.plot2D([2,3],varargin{:}); 
        end
        
        function plotZX(obj,varargin)
            % Plot the evolution in the ZX plane
           obj.plot2D([3,1],varargin{:});           
        end
    end
    
end