% nonlinear pendulum demo
% this shows the phase portraits for the nonlinear pendulum,
%    theta'' + c theta' + k sin(theta) = 0
% around the two critical points (0,0) and (pi, 0), and then combines
% these in the phase plane.
%
% by Gavin LaRose, 3/29/13
%
% (c)2013 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.
%
function Nonlinear_Pendulum()
%--------------------------------------------------------------------------
% configuration variables

    % values of the parameters c and k in the equation
    c = 2;  % values 0 and 2 are nice
    k = 3;  % 3 works well
    
    ltraj_col = 'b';
    ntraj_col = 'm';
    
    % set this to 0 to not have pauses between showing the different
    %    phase portraits
    include_pauses = 1;
    
    bg_grey = [ .6 .6 .6 ];

%--------------------------------------------------------------------------
% initialization

    fwin   = figure( 'Name', 'An Elegant Portrait' );
    
    lpport1 = subplot( 2,2,1, 'Parent', fwin, 'Color', bg_grey );
    lpport2 = subplot( 2,2,2, 'Parent', fwin, 'Color', bg_grey );
    nlpport = subplot( 2,2,[3,4], 'Parent', fwin, 'Color', bg_grey );
    hold(lpport1);
    hold(lpport2);
    hold(nlpport);
    
    % put in some axes
    plot( [-1 1], [0 0], ':w', 'LineWidth', 2, 'Parent', lpport1 );
    plot( [0 0], [-1 1], ':w', 'LineWidth', 2, 'Parent', lpport1 );
    plot( [-1 1], [0 0], ':w', 'LineWidth', 2, 'Parent', lpport2 );
    plot( [0 0], [-1 1], ':w', 'LineWidth', 2, 'Parent', lpport2 );
    plot( [-4 7], [0 0], ':w', 'LineWidth', 2, 'Parent', nlpport );
    plot( [0 0], [-2 2], ':w', 'LineWidth', 2, 'Parent', nlpport );
    axis( lpport1, [ -1 1 -1 1 ] );
    axis( lpport2, [ -1 1 -1 1 ] );
    axis( nlpport, [ -4 7 -2 2 ] );
    
    % add the critical points
    plot( [-pi,0,pi,2*pi],[0,0,0,0], '.k', 'MarkerSize', 25, ...
          'Parent', nlpport );
    
    title( lpport1, 'Linear portrait near (2n\pi,0)', ...
           'FontSize', 14 );
    xlabel( lpport1, 'x = \theta', 'FontSize', 14 );
    ylabel( lpport1, 'y = \theta''', 'FontSize', 14 );
    title( lpport2, 'Linear portrait near ((2n+1)\pi,0)', ...
           'FontSize', 14 );
    xlabel( lpport2, 'x = \theta', 'FontSize', 14 );
    ylabel( lpport2, 'y = \theta''', 'FontSize', 14 );
    title( nlpport, 'Nonlinear phase portrait', 'FontSize', 14 );
    xlabel( nlpport, 'x = \theta', 'FontSize', 14 );
    ylabel( nlpport, 'y = \theta''', 'FontSize', 14 );
    
    pause on;
    disp( 'Press any key to plot first linear phase portrait' );
    pause;
    
%--------------------------------------------------------------------------
% the linearized matrices

    mat1 = [ 0 1; -k -c ];
    mat2 = [ 0 1; k -c ];

%--------------------------------------------------------------------------
% graph the trajectories

% first linearized system
    plot_eigenvectors( lpport1, mat1 );
    plot_solutiontraj( lpport1, mat1, [ '-' ltraj_col ] );
    
    disp( 'Press any key to add to nl system' );
    pause;
    plot_localtraj( nlpport, [0 0], mat1, .5, ['-' ltraj_col] );
    plot_localtraj( nlpport, [2*pi 0], mat1, .5, ['-' ltraj_col] );

    disp( 'Press any key to plot second linear phase portrait' );
    pause
    plot_eigenvectors( lpport2, mat2 );
    plot_solutiontraj( lpport2, mat2, [ '-' ltraj_col ] );

    disp( 'Press any key to add to nl system' );
    pause;
    plot_localtraj( nlpport, [-pi,0], [0 1; k -c], .5, ['-' ltraj_col] );
    plot_localtraj( nlpport, [pi,0], [0 1; k -c], .5, ['-' ltraj_col] );
    
    disp( 'Press any key to add nl solutions' );
    pause;
    
    ic_list = [ ];
    if c == 0
        x_list = [ -3.5 .25 1 2 2.75 3.5 4.5 5.5 6 ];
        for j1 = 1:length(x_list)
            ic_list = vertcat( ic_list, [ x_list(j1), 0 ] );
        end
        y_list = [ 1.5 2 ];
        for j2 = 1:length(y_list)
            ic_list = vertcat( ic_list, [ -4, y_list(j2) ] );
        end
        ic_list = vertcat( ic_list, [ 4.25, -2 ] );
        ic_list = vertcat( ic_list, [ 4, -2 ] );
    else
        x_list  = -3.5:1:6.5;
        y_list  = -2:1:2;
        for j1 = 1:length(x_list)
            ic_list = vertcat( ic_list, [ x_list(j1), 2 ] );
            ic_list = vertcat( ic_list, [ x_list(j1), -2 ] );
        end
        for j2 = 2:length(y_list) - 1
            ic_list = vertcat( ic_list, [ -4, y_list(j2) ] );
            ic_list = vertcat( ic_list, [ 7, y_list(j2) ] );
        end
    end
    
    for j1 = 1:length(ic_list)
        [t_sol, x_sol] = ode45( @(t,x) [ x(2); ...
                                         -k*sin(x(1)) - c*x(2) ], ...
                                [ 0 20 ], ic_list(j1,:) );
        plot( nlpport, x_sol(:,1), x_sol(:,2), ['-' ntraj_col], ...
              'LineWidth', 2 );
    end    
    
% all done
%--------------------------------------------------------------------------
