% Examination of the effect of changing the input frequency, w, or
% spring force (capacitance) coefficient of the forcing term for the
% linear constant-coefficient non-homogeneous system
%    a y'' + b y' + c y = -f0 cos(w t).
% We assume that a and c are positive and b is non-negative.
%
% We graph in one subplot the solution and long-term response
% (particular solution) for the system, and in the second the maximum
% amplitude of the response. Pauses are included between different
% graphs.
%
% by Justin Dunmyre, 8/22/12
% rewritten by Gavin LaRose, 10/11/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 Response_Simple()
%-------------------------------------------------------------------------------
% configuration variables

% we are solving a y'' + b y' + c y = -f0 cos(w t)
    a = 1;
    b = 0.1;
% if we change c when graphing the response, this is updated below
    c = 1;
% initial conditions used in generating the solution: y(0) = y0,
%    y'(0) = yp0
    y0  = 1;
    yp0 = 0;

    f0 = 1;
    
% if we change w when graphing the response, this is updated below
    w  = 1;
    
% the range of times on which to plot the solutions
    t_range = [ 0 20 ];
    
% should we change w or c when exploring the response curves?
    change_var = 'w';
% the range of the variable to explore
    var_range  = [ 0 2 ];
% and the step size we use in exploring this range
    var_step   = .25;
    
% color choices
    current_full_color   = 'k';
    current_partic_color = 'g';
    current_amp_color    = 'r';
    past_soln_color      = [ 0 0 0.6 ];
    past_line_width      = 0.2;
    past_soln_style      = ':';
    background_grey      = [ 0.6 0.6 0.6 ];
    
%-------------------------------------------------------------------------------
% initialization

    scrsz = get(0,'ScreenSize');
    fig0 = figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes( 'Parent', fig0, 'Color', background_grey );
    hold on;
    
% set up subplots
    soln_sub_plot = subplot(2,1,1, 'Parent', fig0, 'Color', background_grey );
    resp_sub_plot = subplot(2,1,2, 'Parent', fig0, 'Color', background_grey );
    hold( soln_sub_plot );
    hold( resp_sub_plot );

% label these
    function c=num2charspc(v)
        if v == 1
            c = ' ';
        else
            c = [ num2str(v) ' ' ];
        end
    end
    if b == 0
        term2 = ' ';
    else
        term2 = [ '+ ' num2charspc(b) 'y\prime' ];
    end
    if c == 0
        term3 = ' ';
    else
        term3 = [ '+ ' num2charspc(c) 'y' ];
    end
    if change_var == 'c'
        term3 = [ '+ c y' ];
        fterm = [ num2str(f0) ' cos(' num2str(w) ' t)' ];
    else
        fterm = [ num2str(f0) ' cos(w t)' ];
    end

    xlabel( 't', 'FontSize', 16, 'Parent', soln_sub_plot );
    ylabel( 'y', 'FontSize', 16, 'Parent', soln_sub_plot );
    xlabel( change_var, 'FontSize', 16, 'Parent', resp_sub_plot );
    ylabel( 'long-term amp', 'FontSize', 16, 'Parent', resp_sub_plot );

    title( [ 'Solution and particular solution to ' ...
             num2charspc(a) 'y\prime\prime' term2 term3 ' = -' fterm ], ...
           'FontSize', 16, 'Parent', soln_sub_plot );
    title( [ 'Long-term response amplitude as a function of ' change_var ], ...
           'FontSize', 16, 'Parent', resp_sub_plot );
    
% the maximum amplitude is 
    var_vals = linspace( var_range(1), var_range(2), 1000 );
    if change_var == 'w'
        amplist = amp( a, b, c, var_vals );
    else
        amplist = amp( a, b, var_vals, w );
    end
    max_amp = max( amplist );
% constrain the maximum to retain some detail
    if max_amp > 10*f0
        max_amp = 10*f0;
    end

% we can then set the axes for the subplots
    axis( soln_sub_plot, [ t_range(1) t_range(2) -max_amp-1 max_amp+1 ] );
    axis( resp_sub_plot, [ 0 var_range(2) 0 max_amp+1 ] );
    
% the vector of variable values we consider, and t values on which to
%    graph the solutions
    var_vec = var_range(1):var_step:var_range(2);
    t_vals  = linspace( t_range(1), t_range(2), 1000 );

% initialize containers for plot handles
    soln_plot_handles = zeros(length(var_vec),4);
    amp_plot_handles  = zeros(length(var_vec),1);
    
% initialize containers for full and particular solutions and the
%    amplitude of the particular solutions
    partic_solutions = zeros( length(var_vec), length(t_vals) );
    full_solutions   = zeros( length(var_vec), length(t_vals) );
    amplitudes       = zeros( length(var_vec), 1 );
    
% get these solutions for each value in the range of parameter values
%    we consider
    for j=1:length(var_vec)
        evalc( [ change_var '=' num2str(var_vec(j)) ] );
        partic_solutions(j,:) = partic_soln( a, b, c, w, t_vals );
        full_solutions(j,:)   = full_soln( a, b, c, w, y0, yp0, t_vals );
        amplitudes(j)         = amp( a, b, c, w );
    end 
    
%-------------------------------------------------------------------------------
% do graphing

% we plot an initial graph and pause, then go through the first couple
%    of graphs in order, pausing after each.  after that we animate the
%    rest

% initial graph
    pause on;
    j = round(length(var_vec)/2);
    soln_plot_handles(j,:) = ...
        plot( t_vals, full_solutions(j,:), ['-' current_full_color], ...
              t_vals, partic_solutions(j,:), ['-' current_partic_color], ...
              [ 0 t_range(2) ], [ amplitudes(j) amplitudes(j) ], ...
                  [ '-' current_amp_color ], ...
              [ 0 t_range(2) ], [ -amplitudes(j) -amplitudes(j) ], ...
                  [ '-' current_amp_color ], 'MarkerSize', 15, ...
              'LineWidth', 2, 'Parent', soln_sub_plot );
    sol_legend_handle = legend( soln_plot_handles(j,1:3), ...
                                { 'full solution y', ...
                                  'particular solution y_p', ...
                                  'long-term amplitude' } );
    if amplitudes(j) ~= 0
        amp_plot_handles(j) = ...
            plot( var_vec(j), amplitudes(j), [ '.' current_amp_color ], ...
                  'MarkerSize', 30 );
        amp_legend_handle = legend( amp_plot_handles(j), ...
                                    [ change_var ' = ' num2str(var_vec(j)) ] );
    end
        
    disp( 'Press any key to continue' );
    pause;
    
    if amplitudes(j) ~= 0
        delete( amp_legend_handle );
        set( amp_plot_handles(j), 'Color', past_soln_color );
    end
    delete( sol_legend_handle );
    set( soln_plot_handles(j,:), 'LineWidth', past_line_width, ...
                      'Color', past_soln_color, 'LineStyle', past_soln_style );
    
    for j=1:min(3, length(var_vec))
        soln_plot_handles(j,:) = ...
            plot( t_vals, full_solutions(j,:), ['-' current_full_color], ...
                  t_vals, partic_solutions(j,:), ['-' current_partic_color], ...
                  [ 0 t_range(2) ], [ amplitudes(j) amplitudes(j) ], ...
                  [ '-' current_amp_color ], ...
                  [ 0 t_range(2) ], [ -amplitudes(j) -amplitudes(j) ], ...
                  [ '-' current_amp_color ], 'MarkerSize', 15, ...
                  'LineWidth', 2, 'Parent', soln_sub_plot );
        sol_legend_handle = legend( soln_plot_handles(j,1:3), ...
                                    { 'full solution y', ...
                                      'particular solution y_p', ...
                                      'long-term amplitude' } );
        if amplitudes(j) ~= 0
            amp_plot_handles(j) = ...
                plot( var_vec(j), amplitudes(j), [ '.' current_amp_color ], ...
                      'MarkerSize', 30 );
            amp_legend_handle = legend( amp_plot_handles(j), ...
                                        [ change_var ' = ' num2str(var_vec(j)) ] );
        end
        
        disp( 'Press any key to continue' );
        pause;
        
        delete( sol_legend_handle );
        set( soln_plot_handles(j,:), 'LineWidth', past_line_width, ...
             'Color', past_soln_color, 'LineStyle', past_soln_style );
        if amplitudes(j) ~= 0
            delete( amp_legend_handle );
            set( amp_plot_handles(j), 'Color', past_soln_color );
        end
    end
    
% finally, go through and plot the rest with minimal pauses
    for j=4:length(var_vec)
        soln_plot_handles(j,:) = ...
            plot( t_vals, full_solutions(j,:), ['-' current_full_color], ...
                  t_vals, partic_solutions(j,:), ['-' current_partic_color], ...
                  [ 0 t_range(2) ], [ amplitudes(j) amplitudes(j) ], ...
                  [ '-' current_amp_color ], ...
                  [ 0 t_range(2) ], [ -amplitudes(j) -amplitudes(j) ], ...
                  [ '-' current_amp_color ], 'MarkerSize', 15, ...
                  'LineWidth', 2, 'Parent', soln_sub_plot );
        sol_legend_handle = legend( soln_plot_handles(j,1:3), ...
                                    { 'full solution y', ...
                                      'particular solution y_p', ...
                                      'long-term amplitude' } );
        if amplitudes(j) ~= 0
            amp_plot_handles(j) = ...
                plot( var_vec(j), amplitudes(j), [ '.' current_amp_color ], ...
                      'MarkerSize', 30 );
            amp_legend_handle = legend( amp_plot_handles(j), ...
                                        [ change_var ' = ' num2str(var_vec(j)) ] );
        end
        pause(.2);
        if j ~= length(var_vec)
            set( soln_plot_handles(j,:), 'LineWidth', past_line_width, ...
                 'Color', past_soln_color, 'LineStyle', past_soln_style );
            delete( sol_legend_handle );
        end
        if amplitudes(j) ~= 0
            set( amp_plot_handles(j), 'Color', past_soln_color );
            delete( amp_legend_handle );
        end
    end;
    
% replace zero amplitudes, which indicate that we're at a singularity,
%    with something that gives a moderately reasonable graph
    zeroindices = find(amplist==0);
% there should be at most one of these
    if length(zeroindices) == 1
        amplist(zeroindices(1)) = 2*max_amp;
    end
    plot( var_vals, amplist, '-k', 'LineWidth', 2 );
    
%-------------------------------------------------------------------------------
% helper functions to find amplitude and solutions

% the amplitude of the particular solution (or zero, if it's unbounded)
%    this allows for c or w to be a vector
    function out = amp( a, b, c, w )
	if b > 0
            out = sqrt( f0^2./(c.^2 + b^2*w.^2 - 2*a*c.*w.^2 + a^2*w.^4) );
        else
            out = abs( f0.*(1 - ((c - a*w.^2)==0)) ./ ...
                       ((c - a*w.^2) + (c - a*w.^2 == 0)) );
 	end
    end

% the particular solution to the equation; all parameters are scalar
%    save the range of t values
    function out = partic_soln( a, b, c, w, tvals )
        if b > 0
            den = (c - a*w^2)^2 + (b*w)^2;
            a0 = (c - a*w^2)*(-1*f0) / den;
            b0 = (-1*b*w*f0) / den;

            out = a0*cos(w*tvals) + b0*sin(w*tvals);
        else
            if c - a*w^2 ~= 0
                out = (-1*f0/(c - a*w^2))*cos(w*tvals);
            
            else
                out = (-1*f0/(2*sqrt(c*a)))*tvals.*sin(w*tvals);
            end
        end
    end

% the full solution to the problem; all parameters are scalar save the
%    range of t values
    function out = full_soln( a, b, c, w, y0, yp0, tvals )

    % the particular solution is
        ypartic = partic_soln( a, b, c, w, tvals );

    % find roots of characteristic polynomial
        char_roots = roots( [a b c] ); 
        realr = real( char_roots(1) );
        imagr = real( char_roots(1) );

        if b > 0
            den = (c - a*w^2) + (b*w)^2;
            % the values of the particular solution and its derivative
            %    at zero are
            p0  = (c - a*w^2)*(-1*f0) / den;
            pp0 = (-1*b*w^2*f0) / den;
            
            if imagr ~= 0
            % solve for c1, c2
                tmpc = [1 0; realr imagr]^-1 * [y0 - p0; yp0 - pp0];
                out = exp(realr*tvals).*( tmpc(1)*cos(imagr*tvals) + ...
                                          tmpc(2)*sin(imagr*tvals) ) + ...
                      ypartic;
            else
                r1 = real( char_roots(1) );
                r2 = real( char_roots(2) );
                
                if r1 ~= r2
                    tmpc = [1 1; r1 r2]^-1 * [y0 - p0; yp0 - pp0];
                    out = tmpc(1)*exp(r1*tvals) + tmpc(2)*exp(r2* ...
                                                              tvals) + ...
                          ypartic;
                else
                    tmpc = [1 0; r1 1]^1 * [y0 - p0; yp0 - pp0];
                    out = tmpc(1)*exp(r1*tvals) + ...
                          tmpc(2)*tvals.*exp(r1*tvals) + ...
                          ypartic;
                end
            end
        else
            if c - a*w^2 ~= 0
                p0 = -1*f0/(c - a*w^2);
            else
                p0 = 0;
            end
            
            if det( [1 0; realr imagr] ) ~= 0
                tmpc = [1 0; realr imagr]^-1 * [y0 - p0; yp0];
            else
                tmpc = [ yp0; yp0/sqrt(c/a) ];
            end
            out = tmpc(1)*cos(imagr*tvals) + tmpc(2)*sin(imagr*tvals) + ...
                  ypartic;
        end
    end
end

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