%-------------------------------------------------------------------------------
% Response_Curve_GUI.m
%
% Explores the effects of changing frequency, w, of the forcing term for the
% mechanical system given by
%
%           m y'' + c y' + k x = F0 cos(wt)
%
% Recommended variables to change for exploration:
%       (Please note that these can also be changed via the control panel 
%       that appears on the left after running the program )
%
%      wmin,wmax,defaultm,defaultc,defaultk,defaultF0,defaultx0,...
%      defaultxprime0,tmax,numw
% 
%
% Running this function displays a figure with two subplots.
%
% First, for the numw values between wmin and wmax
%   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.
%       Previous plots will vanish, but can be recovered by moving the
%       slider in the middle of the figure window
%
%   2) the bottom subplot shows the amplitude of the response as a dot 
%       placed at (w,amp(w)) on the amplitude graph.
%
%   The program shows all of the ws in sequence, and then the user
%   can see a previously shown set of graphs by moving the slider
% 
% 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_Curve_GUI()
%-------------------------------------------------------------------------------
% configuration variables
    
% these are the values that are initially displayed in the GUI
    defaultm  = 1;
    defaultc  = 0.1;
    defaultk  = 1;
    defaultF0 = 1;
       
    wmin = 0.5;
    wmax = 1.5;
    numw = 101;

    tmin = 0;
    tmax = 100; 

% this is y(0)
    defaultx0 = 1;
% and y'(0)
    defaultxprime0 = 0;
    
% colors
    solutionColor    = 'k';
    homogeneousColor = 'b';
    responseColor    = 'k';
    dotColor         = 'r';
    
    bggray=[.6 .6 .6];
    
%-------------------------------------------------------------------------------
% initialization   
    
% set variables so that they can have global scope    
    m  = defaultm;
    c  = defaultc;
    k  = defaultk;
    x0 = defaultx0;
    xprime0 = defaultxprime0;
    w  = wmin;
    F0 = defaultF0;    

    currentInd = 0;
    wrange     = 0;
    myPlots    = 0;
    myDots     = 0;
    myAmps     = 0;
    homPlot    = 0;
    ampPlot    = 0;
    
    pause on
% create the figure windows
    [figh,controls,subplotSol,subplotResp]=drawWindow();        
% create the GUI
    [slider,cwmin,cwmax,wbox]=setupUI();
% run the demonstration
    mainProcess();

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

    function mainProcess()
        
        % hide the slider's value of w (as it will be changing a lot)
        % set(wbox,'Visible','off');
        % establish a vector for w
        wrange = linspace(wmin,wmax,numw);
        % establish a vector for t
        myt = linspace(tmin,tmax,1000);
                    
        % initialize containers for exact solutions and associated amplitudes
        parExact = zeros(numw,length(myt));
        responseAmpMax = zeros(numw,1);
        myPlots = zeros(numw,1);
        myDots = zeros(numw,1);
        myAmps = zeros(numw,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:numw
            w = wrange(j); 
            [parExact(j,:),responseAmpMax(j)] = getSolAndAmp(m,c,k,w,myt);
        end
        
        % set the axis limits now that we know them
        axis( subplotSol, [tmin tmax -max(responseAmpMax)-1 ...
                           max(responseAmpMax)+1] );
        axis( subplotResp, [wmin wmax min(responseAmpMax)-1 ...
                            max(responseAmpMax)] );
      
        % plot the homogeneous solution    
        homPlot = plot( myt, homsolExact, [homogeneousColor '-'], ...
                        'Parent', subplotSol, 'LineWidth', 2 );
        
        % plot the first solution, and amplitude information    
        myPlots(1) = plot( myt, parExact(1,:), [solutionColor '-'], ...
                           'Parent', subplotSol, 'LineWidth', 2 );
        myDots(1) = plot( wrange(1), responseAmpMax(1), [dotColor '.'], ...
                          'MarkerSize', 30, 'Parent', subplotResp );
        myAmps(1) = plot( myt, responseAmpMax(1)*ones(1,length(myt)), ...
                          [dotColor ':'], 'LineWidth', 2, ...
                          'Parent', subplotSol );
        
        % put in some legends
        solnLegendH = legend( [ homPlot, myPlots(1), myAmps(1) ], ...
                              { 'homg sol', ...
                                'part sol', ...
                                'resp amp' } );
        % respLegendH(1) = legend( myDots(1), ['w = ' num2str(wrange(1))] );
        respLegendH = legend( ['w = ' num2str(wrange(1))] );
        
        % plot the full response amplitude graph   
        ampPlot = plot( wrange, responseAmpMax, [responseColor '-'], ...
                        'Parent', subplotResp, 'LineWidth', 2 );
        
        getframe(); % for animation
           
        % for the remaining w, plot solutions
        for j=2:numw
            % hide the previous plots
            set( myPlots(j-1), 'Visible', 'off' );
            set( myDots(j-1), 'Visible', 'off' );
            set( myAmps(j-1), 'Visible', 'off' );

            % plot the current information    
            myPlots(j) = plot( myt, parExact(j,:), [solutionColor '-'], ...
                               'Parent', subplotSol, 'LineWidth', 2);
            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 );
            % update the slider
            set( slider, 'Value', min(numw,j) );

            % set( respLegendH(j-1), 'Visibility', 'off' );
            % respLegendH(j) = legend( myDots(j), ['w = ' num2str(wrange(j))] );
            respLegendH = legend( ['w = ' num2str(wrange(j))] );
            
            getframe(); %animate
        end
        
        % display the slider's value of w   
        % set( wbox, 'String', ['w=' num2str(wrange(numw))] );
        % set( wbox,'Visible','on' );
        
        % tell the slider that the currently displayed w is the last one    
        currentInd = numw;
    end
    % end mainProcess
    
%---------------------------------------
% set up control panel
%    variables are returned so that they can have global scope

    function [slider,cwmin,cwmax,wbox]=setupUI()    
    
        bgroup = uibuttongroup( 'visible', 'off', 'Position', [0 0 1 1], ...
                                'Parent', controls );
        
        %--------------------
        % the control panel buttons and dialogs
        
        % put the equation at the top
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 500 180 60 ], ...
                   'String', 'm y'''' + c y'' + k y = f0 cos(w t)', ...
                   'HandleVisibility', 'off', 'FontSize', 15 );
        
        % add some instructions
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 440 180 60 ], ...
                   'String', 'Set parameters, then click Draw', ...
                   'HandleVisibility', 'off', 'FontSize', 15 );
        
        % the draw and exit buttons
        uicontrol( 'Style', 'PushButton', 'Parent', bgroup, ...
                   'pos', [ 20 400 60 30 ], ...
                   'String', 'Draw', ...
                   'Callback', @startButtonCall );
        uicontrol( 'Style', 'PushButton', 'Parent', bgroup, ...
                   'pos', [ 100 400 60 30 ], ...
                   'String', 'Exit', ...
                   'Callback', @exitButtonCall );
       
        % parameter values:
        %    m
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 360 60 25 ], ...
                   'String', 'm', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 360 100 25 ], ...
                   'String', num2str(m), ...
                   'Callback', @t1cb );
        %    c
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 320 60 25 ], ...
                   'String', 'c', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 320 100 25 ], ...
                   'String', num2str(c), ...
                   'Callback', @t2cb );
        %    k
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 280 60 25 ], ...
                   'String', 'k', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 280 100 25 ], ...
                   'String', num2str(k), ...
                   'Callback', @t3cb );
        %    f0
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 240 60 25 ], ...
                   'String', 'f0', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 240 100 25 ], ...
                   'String', num2str(F0), ...
                   'Callback', @tf0cb );
        %    y0
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 200 60 25 ], ...
                   'String', 'y(0)', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 200 100 25 ], ...
                   'String', num2str(x0), ...
                   'Callback', @t4cb );
        %    yp0
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 160 60 25 ], ...
                   'String', 'y''(0)', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 160 100 25 ], ...
                   'String', num2str(xprime0), ...
                   'Callback', @t5cb );
        % range of t values
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 40 60 25 ], ...
                   'String', 'max t', ...
                   'HandleVisibility','off', 'FontSize',13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 80 40 100 25 ], ...
                   'String', num2str(tmax), ...
                   'Callback', @t8cb );
        
        % range of w values
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 20 120 160 25 ], ...
                   'String', 'range for w', ...
                   'HandleVisibility', 'off', 'FontSize', 13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                   'pos', [ 20 80 80 25 ], ...
                   'String', num2str(wmin), ...
                   'Callback', @wminCallback );
        uicontrol( 'Style', 'text', 'Parent', bgroup, ...
                   'pos', [ 80 80 30 25 ], ...
                   'String', '--', ...
                   'HandleVisibility', 'off', 'FontSize', 13 );
        uicontrol( 'Style', 'edit', 'Parent', bgroup, ...
                    'pos', [ 110 80 80 25 ], ...
                    'String', num2str(wmax), ...
                    'Callback', @wmaxCallback );

        %--------------------
        % the slider
        %    we place the slider relative to the other two subplots
        
        lowerplotPos  = get( subplotResp, 'Position' );
        lowerplotYMax = lowerplotPos(2) + lowerplotPos(4);
        upperplotPos  = get( subplotSol, 'Position' );
        upperplotYMin = upperplotPos(2);
        lowerplotXMin = lowerplotPos(1);
        
        lowerplotWidth  = lowerplotPos(3);
        lowerplotHeight = lowerplotPos(4);
        sliderHeight    = lowerplotHeight/4;
        sliderY         = lowerplotYMax + ...
            (upperplotYMin-lowerplotYMax)/2 - sliderHeight/2;
        posVector       = [ lowerplotXMin, sliderY, lowerplotWidth, ...
                            sliderHeight];
    
        % place the slider
        slider = uicontrol( figh, 'Style', 'Slider', ...
                            'Units', get(subplotResp, 'Units'), ...
                            'Pos', posVector,...
                            'Min', 1, 'Max', numw, ...
                            'SliderStep', [1,1/(numw)], ...
                            'Callback',@wSliderCall, 'Value', 1 );
    
        set( slider, 'Units', 'pixels' );
        sliderPos = get(slider, 'Pos');
        set( slider, 'Units', 'normalized' );

        % place the text boxes flanking the slider
        cwmin = uicontrol( 'Parent', figh, 'Style', 'text', ...
                           'String', num2str(wmin), ...
                           'pos', [ sliderPos(1)-60, sliderPos(2), 50, 30 ],...
                           'BackgroundColor', get(figh,'Color') );
        cwmax = uicontrol( 'Parent', figh, 'Style', 'text', ...
                           'String', num2str(wmax), ...
                           'pos', [ sliderPos(1)+sliderPos(3)+10, ...
                                    sliderPos(2), 50, 30 ], ...
                           'BackgroundColor', get(figh,'Color') );
        
        uicontrol( 'Parent', figh, 'Style', 'text', ...
                   'String', '<- w ->', ...
                   'pos', [ sliderPos(1)+sliderPos(3)/2 - 40, ...
                            sliderPos(2)-30, 80, 30 ], ...
                   'HandleVisibility', 'on', 'FontSize',15, ...
                   'BackgroundColor', get(figh,'Color') );
        
        % wbox = uicontrol( 'Style', 'text', 'Parent', figh, ...
        %                   'pos', [sliderPos(1)+sliderPos(3)/2-40, ...
        %                           sliderPos(2)-60, 120, 30], ...
        %                   'String', ['current w: ' num2str(0)], ...
        %                   'BackgroundColor', get(figh,'Color'), ...
        %                   'FontSize', 15, 'Visible', 'off' );
        wbox = 0;
        
        set(bgroup,'Visible','on');   
    end

%---------------------------------------
% start and exit buttons

    function startButtonCall(hObject,eventdata)
        figure(figh);
        % hide all of the graphs
        set(homPlot,'Visible','off');
        set(ampPlot,'Visible','off');
        set(myPlots(currentInd),'Visible','off');
        set(myDots(currentInd),'Visible','off');
        set(myAmps(currentInd),'Visible','off');
    
        % reset the slider
        set(slider,'Value',1);
        set(slider,'Min',1,'Max',numw);
        set(cwmin,'String',num2str(wmin));
        set(cwmax,'String',num2str(wmax));

        getframe(); %animate to update
        mainProcess();
    end
    
    function exitButtonCall(hObject,eventdata)
        exitFlag=1;
        close all;
    end

% callbacks to set parameters
    function t1cb(hObject,eventdata)
        m=str2double(get(hObject,'String'));
    end
    
    function t2cb(hObject,eventdata)
        c=str2double(get(hObject,'String'));
    end

    function t3cb(hObject,eventdata)
        k=str2double(get(hObject,'String'));
    end

    function t4cb(hObject,eventdata)
        x0=str2double(get(hObject,'String'));
    end
      
    function t5cb(hObject,eventdata)
        xprime0=str2double(get(hObject,'String'));
    end
    
    function t6cb(hObject,eventdata)
        func=get(hObject,'String');
    end

    % function t7cb(hObject,eventdata)
    %    numw=str2double(get(hObject,'String'));
    % end
    
    function t8cb(hObject,eventdata)
        tmax=str2double(get(hObject,'String'));
    end

    function tf0cb(hObject,eventdata)
        F0=str2double(get(hObject,'String'));
    end

    function wminCallback(hObject,eventdata)
        wmin=str2double(get(hObject,'String'));
    end

    function wmaxCallback(hObject,eventdata)
        wmax=str2double(get(hObject,'String'));
    end

%---------------------------------------
% callback function for when slider is moved and released

    function wSliderCall(hObject,eventdata)
        % hide the old plots
        set(myPlots(currentInd),'Visible','off');
        set(myDots(currentInd),'Visible','off');
        set(myAmps(currentInd),'Visible','off');
        % get the current index of the slider    
        currentInd=round(get(hObject,'Value'));
        % show the desired plots and value of w
        % set(wbox,'String',['current w: ' num2str(wrange(currentInd))]);
        set(myPlots(currentInd),'Visible','on');
        set(myDots(currentInd),'Visible','on');
        set(myAmps(currentInd),'Visible','on');
        
        % update the legend in the w graph
        respLegendH = legend( [ 'w = ', num2str(wrange(currentInd)) ] );

    end

%---------------------------------------
% set up the controls and figure windows

    function [figh,controls,subplotSol,subplotResp]=drawWindow()

        scrsz    = get(0,'ScreenSize');
        scrsz(3) = min(1500,scrsz(3));
        scrsz(4) = min(900,scrsz(4));

        % the controls window
        % controls=figure('OuterPosition',[1 scrsz(4)/1.5 240 scrsz(4)/1.25]);
        controls = figure( 'OuterPosition', [1 200 240 600] );
        set( controls, 'MenuBar', 'none', 'Name', 'Controls' );
        
        % the figure window
        figh = figure( 'OuterPosition', ...
                       [1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/1.5],...
                       'Name', 'Response Curves' );
        set( figh, 'Toolbar', 'figure' );
        axes( 'Parent', figh, 'Color', bggray);

        subplotSol  = subplot( 3,1,1, 'parent', figh );
        subplotResp = subplot( 3,1,3, 'parent', figh );

        xlabel( subplotSol,  't', 'FontSize', 15 );
        xlabel( subplotResp, 'w', 'FontSize', 15 );
        ylabel( subplotSol,  'x(t)','FontSize',15);
        ylabel( subplotResp, 'response amplitude', 'FontSize', 15 );

        hold( subplotSol );
        hold( subplotResp );   
    end

%---------------------------------------
% find 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

%---------------------------------------
% find 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 
            % as described in the text
            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
end

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