% Examination of the homogeneous and particular solutions to a second-
% order constant-coefficient linear problem with forcing terms (that
% can be solved with the Method of Undetermined Coefficients).  We
% show two graphs: in the top, the solution to the IVP, y, along with
% the homogeneous and particular solution parts of this; in the bottom,
% the forcing function and the particular solution.
%
% by Justin Dunmyre
% rewritten by Gavin LaRose, 10/5/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 MUC_GUI()
%-------------------------------------------------------------------------------
% configuration variables

% we solve a y'' + b y' + c y = f
    defaulta = 1;
    defaultb = 1;
    defaultc = 1;
    
% the function f.  all of the following are good definitions for f
%    defaultf = '2*cos(0.5*t)';
%    defaultf = '2*t*exp(-t/10)';
%    defaultf = '2*exp(-t/10)*cos(2*t)';
    defaultf = '2*cos(0.5*t)';

% initial conditions to use: y(0) = y0, y'(0) = yp0
    defaulty0  = 1;
    defaultyp0 = 0;
    
% the t values on which to graph the solutions
    t_range = [ 0 100 ];
    
% color choices
    homog_color  = 'b';
    partic_color = 'g';
    full_color   = 'k';
    force_color  = 'c';
    
    background_grey = [ 0.6 0.6 0.6 ];
    
%-------------------------------------------------------------------------------
% initialization

% define initial values for the variables
    a = defaulta;
    b = defaultb;
    c = defaultc;
    
    y0 = defaulty0;
    yp0 = defaultyp0;

    func = defaultf;
    funcVec = regexprep( func, '\*', '.*' );
    funcVec = regexprep( funcVec, '\^', '.^' );
    funcVec = regexprep( funcVec, '\/', './' );
    
    tstart = t_range(1);
    tend   = t_range(2);
    
% set up the figures
    scrsz = get(0,'ScreenSize');
    figh = figure( 'Position', [1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2] );
    fighax = axes( 'Parent', figh, 'Color', background_grey );

    controls = figure( 'OuterPosition', [1 200 205 600] );
    set( controls, 'MenuBar', 'none', 'Name', 'Controls' );
    
    subplotSol = subplot(2,1,1, 'Parent', figh, 'Color', background_grey);
    subplotF   = subplot(2,1,2, 'Parent', figh, 'Color', background_grey);
    hold( subplotSol );
    hold( subplotF );
    
    xlabel( 't', 'FontSize', 16, 'Parent', subplotSol );
    ylabel( 'y', 'FontSize', 16, 'Parent', subplotSol );
    xlabel( 't', 'FontSize', 16, 'Parent', subplotF );
    ylabel( 'y', 'FontSize', 16, 'Parent', subplotF );
    title( 'Full, homogeneous and particular solutions', ...
           'FontSize', 16, 'Parent', subplotSol )
    title( 'Particular solution and forcing term', ...
           'FontSize', 16, 'Parent', subplotF )
    
% flags
    topPlot=0;
    botPlot=0;

% the equation we're solving
    eqn = [ num2str(a) '*D2y + ' num2str(b) '*Dy + ' num2str(c) '*y ' ...
            '= ' func ];
% define a function for the right-hand side
    function out = rhs(t)
        evalc( [ 'out=' funcVec ] );
    end
    
% set up the GUI and run the process
    setupUI();
    mainProcess();
    
%-------------------------------------------------------------------------------
% main process

    function mainProcess()

        % first we get a function from which we can extricate the particular
        %   solution
        % turn off warnings, insofar as we can; we'll test for this later
        warnstate = warning('off','symbolic:mupadmex:MuPADTextWarning');
        ys = dsolve( eqn );
        warnstate = warning('off','symbolic:mupadmex:MuPADTextWarning');

        % the parameters and variables are in this are
        var_list = symvar( ys );
    
        if length(var_list) > 0
            ygeneral = matlabFunction( ys, 'vars', var_list );

        % for the solution with initial condition, we just use ode45
            [ tvals, ysoln ] = ode45( ...
                @(t,y) [ y(2); rhs(t)-(b*y(2)+c*y(1))/a ], ...
                tstart:.1:tend, [ y0, yp0 ] );
        
        % the homogeneous solution, which we define to be the full
        %    solution with the particular solution subtracted off
            homogeneoussol = ysoln(:,1) - ygeneral(0,0,tvals);
        % the full solution
            sol = ysoln(:,1);
            
        % the particular solution
            particularsol = ygeneral(0,0,tvals);
            
        % plot all of these
            topPlot = plot( tvals, sol, [ '-' full_color ], ...
                            tvals, homogeneoussol, [ '-' homog_color ], ...
                            tvals, particularsol, [ '-', partic_color ], ...
                            'LineWidth', 2, 'Parent', subplotSol );
            legend( subplotSol, 'full solution y', ...
                    'homogeneous solution y_c', 'particular solution y_p' );
            botPlot = plot( tvals, rhs(tvals), [ '-', force_color ], ...
                            tvals, particularsol, [ '-', partic_color ], ...
                            'LineWidth', 2, 'Parent', subplotF );
            legend( subplotF, 'forcing f(t)', 'particular solution y_p' );
        else
            warning(warnstate);
            disp( [ 'Ouch!  Matlab failed to solve the differential ' ...
                    'equation.  Is your forcing term one for which we ' ...
                    'can solve this?' ] );
        end
    end

    function setupUI()
	figure( controls );
        
        bgroup = uibuttongroup( 'visible', 'off', 'Position', [0 0 1 1], ...
                                'Parent', controls );

    % controls, from the bottom of the control window
    %    y'(0)
        c6 = uicontrol( 'Style', 'Text', 'String', 'y''(0)', ...
                        'FontSize', 15, 'pos', [20 20 60 30], ...
                        'parent', bgroup, 'HandleVisibility','off' );
        t6 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 20 30 30], 'Callback', @t6cb, ...
                        'String', num2str( defaultyp0) );
        
    %    y(0)
        c5 = uicontrol( 'Style', 'Text', 'String', 'y(0)', 'FontSize', 15, ...
                        'pos', [20 70 60 30], 'parent', bgroup, ...
                        'HandleVisibility', 'off' );
        t5 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 70 30 30], 'Callback', @t5cb, ...
                        'String', num2str(defaulty0) );

    %    c
        c4 = uicontrol( 'Style', 'Text', 'String', 'c', 'FontSize', 15, ...
                        'pos', [20 120 60 30], 'parent', bgroup, ...
                        'HandleVisibility', 'off' );
        t4 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 120 30 30], 'Callback', @t4cb, ...
                        'String', num2str(defaultc) );

    %    b
        c3 = uicontrol( 'Style', 'Text', 'String', 'b', 'FontSize',15, ...
                        'pos', [20 170 60 30], 'parent', bgroup, ...
                        'HandleVisibility', 'off' );
        t3 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 170 30 30], 'Callback', @t3cb, ...
                        'String', num2str(defaultb) );

    %    a
        c2 = uicontrol( 'Style', 'text', 'String', 'a', 'FontSize',15, ...
                        'pos', [20 220 60 30], 'parent', bgroup, ...
                        'HandleVisibility', 'off' );
        t2 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 220 30 30], 'Callback', @t2cb, ...
                        'String', num2str(defaulta) );
        
    %    f(t)
        c1 = uicontrol( 'Style', 'Text', 'String', 'f(t)', 'FontSize', 15,...
                        'pos', [20 270 60 30], 'parent', bgroup, ...
                        'HandleVisibility', 'off' );
        t1 = uicontrol( 'Style', 'edit', 'parent', bgroup, ...
                        'pos', [80 270 100 30], 'Callback', @t1cb, ...
                        'String', defaultf );
            
        startbutton = uicontrol( 'Parent', bgroup, 'Style', 'PushButton',...
                                 'Position', [20 320 60 30], ...
                                 'String', 'Draw', ...
                                 'Callback',@startButtonCall );
        exitbutton = uicontrol( 'Parent', bgroup, 'Style', 'PushButton',...
                                'Position', [100 320 60 30], ...
                                'String', 'Exit', ...
                                'Callback', @exitButtonCall );
   
        c0b = uicontrol( 'Style', 'text', 'parent', bgroup, ...
                         'HandleVisibility', 'off', 'pos', [20 370 180 60], ...
                         'String', ...
                         'Set parameters, then click Draw', ...
                         'FontSize', 15 );

        c0a = uicontrol( 'Style', 'text', 'parent',bgroup, ...
                        'HandleVisibility', 'off',...
                        'pos', [20 450 180 60], ...
                        'String', 'ay''''+by''+cy =f(t)', 'FontSize', 15);

        
        set( bgroup, 'SelectionChangeFcn', @uibuttongroup1_SelectionChangeFcn );
        set( bgroup, 'Visible', 'on' );   
        
        figure(figh);
    
    end

    function exitButtonCall(hObject,eventdata)
        exitFlag=1;
        close all;
    end

    function startButtonCall(hObject,eventdata)
        figure(figh);
        set( topPlot,'Visible','off' );
        set( botPlot,'Visible','off' );
   
        mainProcess();
    end

    function t1cb(hObject,eventdata)
        func = get( hObject,'String' );
        funcVec = regexprep( func, '\*', '.*' );
        funcVec = regexprep( funcVec, '\^', '.^' );
        funcVec = regexprep( funcVec, '\/', './' );

        % redefine a function for the right-hand side
        function out = rhs(t)
            evalc( [ 'out=' funcVec ] );
        end

        % redefine the equation we're solving
        eqn = [ num2str(a) '*D2y + ' num2str(b) '*Dy + ' num2str(c) '*y ' ...
                '= ' func ];
    end
    
    function t2cb(hObject,eventdata)
        a = str2double( get( hObject, 'String' ) );
    end

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

    function t4cb(hObject,eventdata)
        c = str2double( get( hObject, 'String' ) );
    end
      
    function t5cb(hObject,eventdata)
        y0 = str2double( get( hObject, 'String' ) );
    end
    
    function t6cb( hObject,eventdata )
        yp0 = str2double( get( hObject, 'String' ) );
    end

end
