function Lorenz9_analysis(begin,vid,rmin,rmax,step)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% This tool analyzes the 9D Lorenz system in the range of            %
% r_min < r < r_max                                                  %
% Input variables:                                                   %
% begin     -    # of step to begin with                             %
% vid       -    animate flow and temperature field? (0 - No,1 - Yes)%
% rmin      -    minimum r                                           %
% rmax      -    maximum r                                           %
% step      -    stepsize between r                                  %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% clear head                                                         %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
clc
close all

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% global declaration                                                 %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

dim = 9;

if begin == 1
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% The first step simulates the System                                %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% declaration                                                        %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% initial condition
    x0 = [1 1 1 1 1 1 1 1 1];

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% simulation call (initial_condition, r, t_end)                      %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
    for r = rmin:5:rmax
        [tv, vv] = Lorenz9_sim(x0,r,100);
    end

 end

if begin == 2
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% The second step calculates the lyapunov exponents                  %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% declaration                                                        %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

    i = 0;
    nrs = (rmax-rmin)/step+1;
    lyap = zeros(dim,nrs);
    
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% lyapunov call (Dimension, # Iterations, r, Trans. Phase, t_step)   %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
  disp('Running Lyapunov calculations')  
    for r=rmin:step:rmax
        i = i+1;
        perc = i/nrs*100;
        fprintf('Currently at %3.1f %% of calculation \n',perc)
        tic
        lyap(:,i) = Lorenz9_lyap(dim,1000,r,50,1);
        toc  
    end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% visualization                                                      %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

    r=rmin:step:rmax;
    
    figure;
    plot(r,lyap);
    title('all 9 Lyapunov Exponents');
    xlabel('r');
    ylabel('\lambda');
    legend('\lambda _1','\lambda _2','\lambda _3','\lambda _4', '\lambda _5', '\lambda _6', ...
        '\lambda _7', '\lambda _8', '\lambda _9');
    
    figure;
    plot(r,lyap(1:3,:));
    title('Leading 3 Lyapunov Exponents');
    xlabel('r');
    ylabel('\lambda');
    legend('\lambda _1','\lambda _2','\lambda _3');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% The third step analyzes the lyapunov exponents                     %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% analysis call (Lyapunov Matrix, Dimension)                         %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

    [diss, dlyap h_mu] = Lorenz9_lyan(lyap,dim);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% visualization                                                      %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %    
    
    figure;
    subplot(3,1,1);
    plot(r,dlyap);
    title('Lyapunov Dimension of the attractor');
    xlabel('r');
    ylabel('d_{lyap}');

    subplot(3,1,2);
    plot(r,diss);
    title('Dissipation rate');
    xlabel('r');
    ylabel('D');
    
    subplot(3,1,3);
    plot(r,h_mu);
    title('Entropy rate');
    xlabel('r');
    ylabel('h_{\mu}');
end 

if vid == 1    
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% Simulation of Flow and Temperature Field (v) for rmax              %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %    
         field_sim(tv,vv,dim)
end