%SDSU: MATH 636
%By:   Greg Behm
%Date: Nov 24, 2015

function [] = HWK636BesselJ2()

    %period
    p = pi/2:2*pi:21*pi/2;

    %Define the grid
    [th,r] = meshgrid((0:5:360)*pi/180,0:.05:1);
    [x,y] = pol2cart(th,r);

    %speed/frequency 
    c = 1;

    %The Modes You Want
    Modes = [1,3,4]; %the FB-Modes present in our analysis
    Roots = [7.01559 6.38016 7.58834]; %roots of besselj phi_12,31,41 respectively
    u=zeros(size(x));
    
    %Evolve in time and calls the subroutine GBplot to conserve on space
    %and show how to utilize helper functions inside the main function.
    for t = p(1):(p(end)-p(1))/1000:p(end)
        u=zeros(size(x));
        if (t<=p(2))
            %Marginal Stability Curves indicate the order of modes as followed
            %phi_31 --> phi_12,31 --> phi_12,31,41 --> phi_12,41 --> phi_41  
            for j=2:2
                u = u + besselj(Modes(j),Roots(j)*r).*cos(Modes(j)*th)*cos(c*Roots(j)*t);
            end
            titleString = ['Mode(s)  \Phi_{41} Present   t = ' num2str(t,'%.4f')];
            GBplot(x,y,u,titleString)
        end
        if (p(2)<=t && t<=p(3))
            for j=1:2
                u = u + besselj(Modes(j),Roots(j)*r).*cos(Modes(j)*th)*cos(c*Roots(j)*t);
            end
            titleString = ['Mode(s)  \Phi_{12}   \Phi_{31} Present   t = ' num2str(t,'%.4f')];
            GBplot(x,y,u,titleString)
        end
        if (p(3)<=t && t<=p(4))
            for j=1:3
                u = u + besselj(Modes(j),Roots(j)*r).*cos(Modes(j)*th)*cos(c*Roots(j)*t);
            end
            titleString = ['Mode(s)  \Phi_{12}  \Phi_{31}  \Phi_{41} Present    t = ' num2str(t,'%.4f')];
            GBplot(x,y,u,titleString)
        end
        if (p(4)<=t && t<=p(5))
            for j=[1 3]
                u = u + besselj(Modes(j),Roots(j)*r).*cos(Modes(j)*th)*cos(c*Roots(j)*t);
            end
            titleString = ['Mode(s)  \Phi_{12}  \Phi_{41} Present    t = ' num2str(t,'%.4f')];
            GBplot(x,y,u,titleString)
        end
        if (p(5)<=t)
            for j=3:3
                u = u + besselj(Modes(j),Roots(j)*r).*cos(Modes(j)*th)*cos(c*Roots(j)*t);
            end
            titleString = ['Mode(s)  \Phi_{41} Present    t = ' num2str(t,'%.4f')];
            GBplot(x,y,u,titleString)
        end
        
    end
end

%Subfunction to save on copy paste and to show anyone who hasn't seen subfunctions before        
function [] = GBplot(x,y,u,titleString)
    %One way to plot a global title in subplots
    clf %clf easiest way to remove annotation it shouldn't be visually noticed
    annotation('textbox', [0 0.9 1 0.1], ...
                   'String', titleString, ...
                   'EdgeColor', 'none', ...
                   'HorizontalAlignment', 'center');
               
    %side profile
    subplot(1,2,1)
    mesh(x,y,u);
    zlim([-2 2])
    axis off
    view(-5,30)
    
    %Top down view
    subplot(1,2,2)
    surf(x,y,u);
    zlim([-3 3])
    grid off
    axis square
    axis off
    view(-90,90)
    pause(0.05)
end
