% Examination of the behavior of a second-order linear constant
% coefficient equation, animating a figure showing a mass-spring
% system and the corresponding solution to the differential equation.
%
% parts by Justin Dunmyre
% remainder and rewriting by Gavin LaRose, 9/27/2012

% (c)2012 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 Spring_Animation()

%-------------------------------------------------------------------------------
% configuration variables

% our differential equation is m x'' + c x' + k x = 0, with initial
%    conditions x(0) = x0 and x'(0) = xp0
    m = 1;
    c = 1;
    k = 8; 

% we give a list of initial conditions; only one of x0 or xp0 should
%    be a list of more than one item (if both are lists, only x0 will
%    be varied).  we pause between the graph with each of the specified
%    values of x0 or xp0
    x0  = [ 1 2 3 4 5 ];
    xp0 = [ 0 ];
    
% range of t values to solve for
    t_range = [ 0 10 ];
% number of points at which to generate animation frames; increasing
%    this smooths out the animation but makes it slower
    num_points = 150;
% number of frames to skip by when updating display; this controls the
%    speed of the animation.
    frame_skip = 1;
    
% the color of the background of the display
    bggray = [ 0.6 0.6 0.6 ];
    
% the color of the solution curve
    sol_color  = 'r';
% and the color used for the traces of the previous solutions
    fade_color = [ 0.2 0.2 0.2 ];

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

% first we make sure that the lists of x0 and xp0 values are what we expect
    if length( x0 ) > 1
        xp0 = [ xp0(1) ];
    end

% initialization of spring diagram data
    springGraphInitialized = 0;
    springCoilh = 0;
    springDoth  = 0;
    dashpoth    = 0;
    d1=0;d2=0;d3=0;d4=0;d5=0;
% and solution data
    solTail = 0;

% the maximum displacement is determined by x0 and xp0 and the roots
%    of the characteristic equation.  we assume that real(roots) < 0
%    (that is, we have sensible damping).  this is a bit messy.
    char_roots = roots( [ m c k ] );
    
    if char_roots(1) == char_roots(2)
        % then we have repeated roots, x = a e^(rt) + b t e^(rt),
        %    and x has an extremum at
        r = char_roots(1);
        if length(xp0) == 1
            [t_check, ind] = xp0(1) / (r*(r*x0 - xp0(1)));
            x0_check  = x0(ind);
            xp0_check = xp0(1);
        else
            [t_check, ind] = xp0 ./ (r*(r*x0(1) - xp0));
            x0_check  = x0(1);
            xp0_check = xp0(ind);
        end
        if max(t_check) < 0
            spring_max   = round(max(abs(x0))) + 1;
        else
            amp_check = abs(x0_check*exp(r*t_check) + ...
                            t_check*exp(r*t_check)*(xp0_check - r*x0_check));
            spring_max = round(max(horzcat( abs(x0), amp_check ))) + 1;
        end
        
    elseif imag(char_roots(1)) ~= 0
        % imaginary roots, x = a e^(rt) cos(w t) + b e^(rt) sin(w t),
        %    and the max amplitude is a^2 + b^2
        r = real(char_roots(1));
        w = abs(imag(char_roots(1)));
        if length(xp0) == 1
            amp_vals = sqrt( x0.^2 + (-r*x0 + xp0(1)).^2/w^2 );
        else
            amp_vals = sqrt( x0(1)^2 + (-r*x0(1) + xp0).^2/w^2 );
        end
        spring_max = round(max( amp_vals )) + 1;
        
    else
        % distinct real roots, x = a e^(r1 t) + b e^(r2 t); the max
        %    is either at t=0 or at an extremum
        r1 = min(char_roots);
        r2 = max(char_roots);
        if length(xp0) == 1
            if max(r1*x0 - xp0(1)) < 0
                spring_max = round(max(abs( (xp0(1)-r1*x0) + ...
                                            (r2*x0-xp0(1)) )/(r2-r1) ) ) + 1;
            else
                arg = r1*(r2*x0 - xp0(1))/(r2*(r1*x0 - xp0(1)));
                t_check = (arg > 0) .* log(arg)./(r2 - r1);
                amp_vals = abs( (exp(r2*t_check).*(-r1*x0 + xp0(1)) + ...
                                 exp(r1*t_check).*(r2*x0 - xp0(1))) ) ./ ...
                    (r2 - r1);
                zero_vals = abs( (-r1*x0 + xp0(1)) + (r2*x0 - xp0(1)) ) ./ ...
                    (r2 - r1);
                spring_max = round(max(horzcat( amp_vals, zero_vals ))) + 1;
            end
        else
            if max(r1*x0(1) - xp0) < 0
                spring_max = round(max(abs( (xp0-r1*x0(1)) + ...
                                            (r2*x0(1)-xp0) )/(r2-r1) ) ) + 1;
            else
                arg = r1*(r2*x0(1) - xp0)/(r2*(r1*x0(1) - xp0));
                t_check = (arg > 0) .* log(arg)./(r2 - r1);
                amp_vals = abs( (exp(r2*t_check).*(-r1*x0(1) + xp0) + ...
                                 exp(r1*t_check).*(r2*x0(1) - xp0)) ) ./ ...
                    (r2 - r1);
                zero_vals = abs( (-r1*x0(1) + xp0) + (r2*x0(1) - xp0) ) ./ ...
                    (r2 - r1);
                spring_max = round(max(horzcat( amp_vals, zero_vals ))) + 1;
            end
        end
    end
    spring_range = [ -1*spring_max spring_max ];
    
% set up the figure
    scrsz = get( 0,'ScreenSize' );
    figh  = figure( 'OuterPosition', [1 scrsz(4)/2 scrsz(3)/3 scrsz(4)/3], ...
                    'Name', 'Spring Simulation' );
    axes( 'Parent', figh, 'Color', bggray);
    
    spring_sub_plot = subplot(2,1,1,'parent',figh);
    soln_sub_plot   = subplot(2,1,2,'parent',figh); 

    hold(spring_sub_plot);
    hold(soln_sub_plot);
    
    axis( spring_sub_plot, [ spring_range(1) spring_range(2) -1 1 ] );
    axis( soln_sub_plot, [ t_range(1) t_range(2) spring_range(1) ...
                           spring_range(2) ] );
    
    plot( zeros(1,100), linspace(-1,1,100), 'k:', 'Parent', spring_sub_plot, ...
          'LineWidth', 2 );
    plot( linspace( t_range(1), t_range(2), 100), zeros(1,100), 'k:' );
    
    xlabel( 't', 'FontSize', 14, 'Parent', soln_sub_plot );
    ylabel( 'x', 'FontSize', 14, 'Parent', soln_sub_plot );
    title( 'Spring Animation', 'FontSize', 16, 'Parent', spring_sub_plot );
    title( 'Solution Graph', 'FontSize', 16, 'Parent', soln_sub_plot );

% do plot for each set of initial conditions
    for i = 1:length(x0)
        for j = 1:length(xp0)
            
            if(length(solTail)>1)
                set( solTail(end), 'Color', fade_color, ...
                                  'LineStyle','-', 'LineWidth', 1.5 );
                set( solDot(end), 'Visible', 'off' );
            end

            x_0  = x0(i);
            xp_0 = xp0(j);
            
            plotSystem( x_0, m, c, k );

            [t, springSol] = ode45( @(t,x) [x(2); (-k*x(1)-c*x(2))/m], ...
                linspace( t_range(1), t_range(2), num_points ), ...
                [ x_0; xp_0 ] );
        
            springX = springSol(:,1);

            solDot(1) = plot( t(1), springX(1), [sol_color '.'], ...
                              'MarkerSize', 15, 'Parent', soln_sub_plot );
            solTail(1) = plot( t(1:2), springX(1:2), [sol_color '-'], ...
                               'LineWidth', 1, 'Parent', soln_sub_plot );
            getframe();   

            for jj=2:frame_skip:length(t)
                set( solTail(jj-1), 'Visible', 'off' );
                set( solDot(jj-1), 'Visible', 'off' );

                solTail(jj) = plot( t(1:jj), springX(1:jj), ...
                                    [sol_color '-'], 'LineWidth', 3, ...
                                    'Parent', soln_sub_plot );
                solDot(jj) = plot( t(jj),springX(jj),[sol_color '.'], ...
                                   'MarkerSize', 40, 'Parent', soln_sub_plot );

                plotSystem( springX(jj), m, c, k );
                getframe();
            end
            
            if length(xp0) > 1 & j < length(xp0)
                pause
            end
        end
        if length(x0) > 1 & i < length(x0)
            pause
        end
    end
    
    function plotSystem( x, m, c, k )
        
        if ( springGraphInitialized > 0 )
            delete(springCoilh);
            delete(springDoth);
            % for some reason we can't delete these
            set( d1, 'Visible', 'off' );
            set( d2, 'Visible', 'off' );
            set( d3, 'Visible', 'off' );
            set( d4, 'Visible', 'off' );
            set( d5, 'Visible', 'off' );
            % delete( d1,d2,d3,d4,d5 );
        else
            springGraphInitialized = 1;
        end
        
        orig    = spring_range(1);
        numCoil = 10;
        amp     = .4;
        coilx   = linspace(orig,x,10000);
        tempx   = (coilx-orig)*(pi-orig)/(x-orig) + orig;
        coily   = amp*sin(numCoil*tempx);

        springCoilh = plot( coilx, coily, 'b-', 'Parent', spring_sub_plot );
        
        if ( c>0 )
            % if( springGraphInitialized>0 )
            %     delete(d1,d2,d3,d4,d5); 
            %     % d1=0;d2=0;d3=0;d4=0;d5=0;
            % end
            origin     = x;
            stemLength = 1;
            potAmp     = .4;
            totalx     = spring_max;
            realOrigin = -2;
            endstem    = origin+stemLength;

            d1 = plot( [origin, endstem], [0 0], 'k-', 'LineWidth', 2, ...
                       'Parent', spring_sub_plot );
            d2 = plot( [endstem endstem], [-potAmp potAmp], 'k:', ...
                       'LineWidth', 2, 'Parent', spring_sub_plot );
            d3 = plot( [realOrigin+stemLength/2, totalx], [potAmp potAmp], ...
                       'k-', 'Parent', spring_sub_plot );
            d4 = plot( [realOrigin+stemLength/2, totalx], ...
                       [-potAmp -potAmp], 'k-', 'Parent', spring_sub_plot );
            d5 = plot( [totalx totalx], [-potAmp, potAmp], 'k-', ...
                       'Parent', spring_sub_plot );        
        % else
            % if( springGraphInitialized>0 )
            %     delete( d1,d2,d3,d4,d5 );
            %     % d1=0;d2=0;d3=0;d4=0;d5=0;
            % end
        end
        
        springDoth = plot( x, 0, 'r.', 'MarkerSize', 80, ...
                           'Parent', spring_sub_plot );
    end
end

% all done
%-------------------------------------------------------------------------------
