%-------------------------------------------------------------------------------
% Response_Simple2.m
%
% Explores the effects of changing frequency, w, of the forcing term for the
% mechanical system given by
%
%           m y'' + c y' + k y = F0 cos(wt)
%
% Recommended variables to change for exploration:
%      wrange,wmin,wmax,m,c,k,F0,x0,xprime0,tmax,widthScaler
%
% Running this function displays a figure with two subplots.
%
% First, for each value included in wrange
%   1) the top subplot shows the analytical solution to the ODE x vs t, as 
%       well as the amplitude of the solution particular to the current w
%       Previously displayed solutions will become harder to see (see
%       widthScaler variable)
%   2) the bottom subplot shows the amplitude of the response as a dot 
%       placed at (w,amp(w))
%   The program pauses after each w, until the user strikes any key
% 
% Finally, the full graph for the amplitude of the response is displayed in 
% the lower subplot.
% 
% by Justin Dunmyre 8/22/12
% lightly edited 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_Simple2()
%-------------------------------------------------------------------------------
% configuration variables

% individual values of w at which we will plot solutions
    wrange = [.5 .75 1 1.25 1.5];
    
% range of w on which to plot the full amplitude response graph
    wmin = 0;
    wmax = 1.5;
    
% values of m, c, k, F0, y(0) and y'(0) in the equation
%    m y'' + c y' + k y = F0 cos(w t)
    m = 1;    
    c = 0.1;    
    k = 1;    
    F0 = 1;
% this is y(0)
    x0 = 1;
% and y'(0)
    xprime0 = 0;
    
% start and end times to plot the solution
    tmin = 0;
    tmax = 100;
    
% previously viewed plots will have their visibility reduced by this
% factor.  Set to 0 if you wish to see only the current plot 
    widthScaler = 0.3;
    
% colors
    solutionColor    = 'k';
    homogeneousColor = 'b';
    responseColor    = 'k';
% color to draw the bar for the amplitude of the response
    dotColor         = 'r';
    
% deep blue color for the previously viewed plots
    offColor = [0 0 0.4];
% gray background to use for the figure (heightens visibility of bright
%    colors)
    bggray = [.6 .6 .6];
    
%-------------------------------------------------------------------------------
% initialization

% turn on pausing, draw the figures
    pause on;
    [figh,subplotSol,subplotResp,wbox] = drawWindow();
% run the demonstration
    mainProcess();

%-------------------------------------------------------------------------------
% functions

    function mainProcess()
        
        % initialize containers for plot handles
        myPlots = zeros(length(wrange),1);
        myDots  = zeros(length(wrange),1);
        myAmps  = zeros(length(wrange),1);
        % initialize absissca vectors for solution and amplitude graphs
        myt  = linspace(tmin,tmax,1000);
        ampw = linspace(wmin,wmax,100);
        
        % initialize ordinate vectors for solution and amplitude graphs for
        %    the individual frequencies in wrange
        parExact = zeros(length(wrange),length(myt));
        responseAmpMax = zeros(length(wrange),1);
            
        % get the homogeneous solution
        homsolExact = getHomogeneous(m,c,k,myt);
        
        % for each w, get the particular soltuions along with their amplitudes
        for j=1:length(wrange)
            w=wrange(j); 
            [parExact(j,:),responseAmpMax(j)]=getSolAndAmp(m,c,k,w,myt);        
        end
        
        % get the full amplitude graph
        fullresponseAmpMax = getFullResponseAmp(ampw);
 
        % Setup the figure axis, now that we know the bounds
        axis( subplotSol, [ tmin tmax -max(responseAmpMax)-1 ...
                            max(responseAmpMax)+1] );
        axis( subplotResp, [ wmin wmax min(fullresponseAmpMax)-1 ...
                            max(fullresponseAmpMax)] );
        
        % plot the homogeneous solution
        homPlot = plot( myt,homsolExact, [homogeneousColor '-'], ...
                        'Parent', subplotSol, 'LineWidth', 2);
        
        % plot the first particular solution and amplitude bar in the top
        %    subplot, and the corresponding dot in the lower subplot
        myPlots(1) = plot( myt, parExact(1,:), [solutionColor '-'], ...
                           'Parent', subplotSol, 'LineWidth', 2 );
        myAmps(1)  = plot( myt, responseAmpMax(1)*ones(1,length(myt)), ...
                           [dotColor ':'], 'LineWidth', 2, ...
                           'Parent', subplotSol );
        myDots(1)  = plot( wrange(1), responseAmpMax(1), [dotColor '.'], ...
                           'MarkerSize', 30, 'Parent', subplotResp );
        
        % put in some legends
        solnLegendH = legend( [ homPlot, myPlots(1), myAmps(1) ], ...
                              { 'homogeneous solution', ...
                                'particular solution', ...
                                'response amplitude' } );
        respLegendH = legend( myDots(1), ['w = ' num2str(wrange(1))] );
        
        % plot the rest of the particular solutions and amplitudes as above
        for j=2:length(wrange)
            pause();     

            % make the previous plots less visible
            for ell=1:j-1
                widths = get(myPlots(ell),'LineWidth');
                set( myPlots(ell), 'LineStyle',':', ...
                     'LineWidth', max(widths*widthScaler,1) );
                if(widthScaler==0)
                    set( myPlots(ell), 'Visible', 'off' );
                end
            end
            % change the color of the previous set of plots
            set( myDots(j-1), 'Color', offColor);
            set( myPlots(j-1), 'Color', offColor);
            set( myAmps(j-1), 'Visible', 'off');
            
            % reset the legends
            delete( solnLegendH );
            delete( respLegendH );

            % plot the particular solution and amplitudes as above
            myPlots(j) = plot( myt, parExact(j,:), [solutionColor '-'], ...
                               'Parent', subplotSol, 'LineWidth',3 );
            myDots(j)  = plot( wrange(j), responseAmpMax(j), ...
                               [dotColor '.'], 'MarkerSize', 30, ...
                               'Parent', subplotResp );
            myAmps(j)  = plot( myt, responseAmpMax(j)*ones(1,length(myt)), ...
                               [dotColor ':'], 'LineWidth', 2, ...
                               'Parent', subplotSol );

            % add the legends back in
            solnLegendH = legend( [ homPlot, myPlots(j), myAmps(j) ], ...
                                  { 'homogeneous solution', ...
                                    'particular solution', ...
                                    'response amplitude' } );
            respLegendH = legend( myDots(j), ['w = ' num2str(wrange(j))] );
            
            % Update the value of w in the middle of the figure
            % set( wbox, 'String', ...
            %      ['w=' num2str(wrange(j)) ' ... press enter to continue'] );
            % set( wbox, 'Visible', 'on');
        end
        pause();
        
        % update the value of w in the middle of the figure    
        set( wbox, 'Visible', 'off' );
        
        % draw the full amplitude plot
        ampPlot = plot( ampw, fullresponseAmpMax, [responseColor '-'], ...
                        'Parent', subplotResp, 'LineWidth', 2 );
    end
    
%----------------------------------------
% set up figure

    function [figh,subplotSol,subplotResp,wbox]=drawWindow()
        scrsz = get(0,'ScreenSize'); %Get the screen size
        scrsz(3)=min(1500,scrsz(3));
        scrsz(4)=min(900,scrsz(4));

        figh = figure( 'OuterPosition', ...
                       [200 scrsz(4)/1.5 scrsz(3)/1.5 scrsz(4)/1.5], ...
                       'Name', 'Response Curve Simulation' );
        axes( 'Parent', figh, 'Color', bggray );

        subplotSol = subplot( 2,1,1, 'Parent', figh );
        set( subplotSol, 'Color', bggray );
        subplotResp = subplot( 2,1,2, 'Parent', figh );
        set( subplotResp, 'Color', bggray );
        hold(subplotSol);
        hold(subplotResp);
        
        title(subplotSol,[num2str(m) 'y'''' + ' num2str(c) 'y'' + ' ...
                          num2str(k) 'y= ' num2str(F0) 'cos(wt)'], ...
              'FontSize',15);

        xlabel( subplotSol, 't', 'FontSize', 15 );
        xlabel( subplotResp, 'w', 'FontSize', 15 );
        ylabel( subplotSol, 'y(t)', 'FontSize', 15);
        ylabel( subplotResp, 'response amplitude', 'FontSize', 15 ); 
        
        wbox = uicontrol( 'Style', 'text', 'Parent', figh, ...
                          'FontSize', 15, ...
                          'BackgroundColor', get(figh,'Color'), ...
                          'Units', 'normalized', ...
                          'String', '(press enter to continue)', ...
                          'Visible', 'On', ...
                          'pos', [.3 .45 .4 .05] );
        set( figh, 'Toolbar', 'figure' );
        %                   'String', ['w= ' num2str(wrange(1)) ...
        %                     ' ... press enter to continue'], 'Visible', 'On');
    end

%----------------------------------------
% the amplitude of the graph

    function fullresponseAmpMax=getFullResponseAmp(ampw)
        
        for j=1:length(ampw)
            w=ampw(j); 
            if( c>0 ) % damping case
                
                A = (k-m*w^2)*F0/((k-m*w^2)^2+(c*w)^2);
                B = c*w*F0/((k-m*w^2)^2+(c*w)^2);
                fullresponseAmpMax(j)=sqrt(A^2+B^2);                
                
            else   % undamped case
                if (k-m*w^2~=0) %natural and w frequencies do not match   
                    A=F0/(k-m*w^2);
                    fullresponseAmpMax(j)=abs(A);
                    
                else %frequencies match -> resonance
                   w0=sqrt(k/m);
                   A=0;
                   B=F0/(2*m*w0);
                   fullresponseAmpMax(j)=B*tmax;
                end
            end
         end 
    end

%----------------------------------------
% the particular solution and associated amplitude graph

    function [parExact,AmpMax]=getSolAndAmp(m,c,k,w,myt)
        % find roots of characteristic polynomial
        % a = real part of root, b = imaginary part of root
        myRoots = roots([m c k]); 
        a = real(myRoots(1));
        b = imag(myRoots(1));
            
        % Determine what kind of solutions we will have
        if ( c > 0 ) %damping 
            A = (k-m*w^2)*F0/((k-m*w^2)^2+(c*w)^2);
            B = c*w*F0/((k-m*w^2)^2+(c*w)^2);

            % form particular solution
            xp = A*cos(w*myt) + B*sin(w*myt);
                
            AmpMax = sqrt(A^2+B^2);
            % what kind of damping?
            if (b~=0) % underdamped
                % solve for c1 c2
                tempc = [1 0; a b]^-1*[x0-A;xprime0-B*w];
                c1 = tempc(1);
                c2 = tempc(2);
                % form solution
                parExact = exp(a*myt).*(c1*cos(b*myt)+c2*sin(b*myt)) + xp;
            else
                r1 = real(myRoots(1));
                r2 = real(myRoots(2));
                if( r1~=r2 ) % overdamped
                    % solve for c1 c2
                    tempc = [1 1; r1 r2]^-1*[x0-A;xprime0-B*w];
                    c1 = tempc(1);
                    c2 = tempc(2);
                    % form solution
                    parExact = c1*exp(r1*myt) + c2*exp(r2*myt) + xp;
                else % critically damped
                    % solve for c1 c2
                    tempc = [1 0; r1 1]^-1*[x0-A;xprime0-B*w];
                    c1 = tempc(1);
                    c2 = tempc(2);
                    % form solution
                    parExact = c1*exp(r1*myt) + c2*myt.*exp(r1*myt) + xp;
                end
            end                
        else % no damping
            if ( k-m*w^2 ~= 0 ) % frequencies do not match
                A = F0/(k-m*w^2);
                B = 0;
                AmpMax = abs(A);
                xp = A*cos(w*myt);
                % solve for c1 c2
                tempc = [1 0; a b]^-1*[x0-A;xprime0];
                c1 = tempc(1);
                c2 = tempc(2);
                % form solution
                parExact = exp(a*myt).*(c1*cos(b*myt) + c2*sin(b*myt)) + xp;
            
            else % frequencies do match -> resonance
                w0 = sqrt(k/m);
                A = 0;
                B = F0/(2*m*w0);
                AmpMax = B*tmax;                    
                xp = B*myt.*sin(w*myt);
                
                % solve for c1 c2
                tempc = [1 0; a b]^-1*[x0;xprime0];
                c1 = tempc(1);
                c2 = tempc(2);
                % form solution
                parExact = exp(a*myt).*(c1*cos(b*myt) + c2*sin(b*myt)) + xp;
            end
        end
    end

%----------------------------------------
% the homogeneous solution

    function homsolExact=getHomogeneous(m,c,k,myt)
        % find roots of characteristic polynomial
        % a = real part of root, b = imaginary part of root        
        myRoots = roots([m c k]);
        a = real(myRoots(1));
        b = imag(myRoots(1));
        
        if ( b~=0 ) % underdamped
            % solve for c1 c2
            tempc = [1 0;a b]^-1*[x0;xprime0];
            homc1 = tempc(1);
            homc2 = tempc(2);
            % form solution
            homsolExact = exp(a*myt).*(homc1*cos(b*myt) + homc2*sin(b*myt));
        else
           r1 = real(myRoots(1));
           r2 = real(myRoots(2));
           if ( r1~=r2 ) % overdamped
               % solve for c1 c2
               tempc = [1 1; r1 r2]^-1*[x0;xprime0];
               homc1 = tempc(1);
               homc2 = tempc(2);
               % form solution
               homsolExact = homc1*exp(r1*myt) + homc2*exp(r2*myt);
           else % critically damped
               % solve for c1 c2
               tempc = [1 0; r1 1]^-1*[x0;xprime0]; 
               homc1 = tempc(1);
               homc2 = tempc(2);
               % form solution
               homsolExact = homc1*exp(r1*myt) + homc2*myt.*exp(r1*myt);
           end
        end
    end
end

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