function [tout,xout] = bead_on_hoop(X0, show_animation)

% Run this code by using the m-file "run_bead_on_hoop.m"
%
% NONLINEAR PHENOMENA, Spring 2013, UCSB (ECE 183 / ME 169 / PHYS 106)
% Katie Byl
% Lecture 8 (April 25, 2013)
%
% Code to demonstrate the "bead on a hoop", presented in Lecture 6.

global bfac

if ~exist('show_animation','var')
    show_animation = false;
end
if show_animation
    S = odeset('AbsTol',1e-12,'RelTol',1e-12);
else
    S = odeset('AbsTol',1e-6,'RelTol',1e-6); % to speed up ode solver
end
[tout,xout] = ode45(@bead_ode,[0 100],X0,S);

g = 9.8;
r = 1.0; % ring radius
m = 0.1; % mass of bead
w = 4; % angular velocity of hoop spinning
b = sqrt(bfac*m*m*g*r);  % we want "b^2" >> "m^2*g*r", for overdamped case
theta_crit = acos(g/(r*w*w));

if show_animation
    figure(2); clf
    plot(tout,xout(:,1)); hold on
    plot(tout,0*tout+theta_crit,'r--')
    plot(tout,0*tout-theta_crit,'r--')
    figure(4); clf
    %ti = interp1(tout,1:length(tout),10);
    %ti = ceil(ti);
    %plot(tout(ti:end),xout(ti:end,1),'b');
    plot(tout,xout(:,1)-xout(end,1)); hold on
    axis([0 50 xout(end,1)*0+1e-6*[-1 1]]);
    grid on
    title('Zooming in to look for Oscillations...')
    xlabel('Time (s)')
    ylabel('x(t) - x_{final}');
end

figure(3);
plot(xout(:,1),xout(:,2),'k-','LineWidth',1); hold on
plot(xout(end,1),xout(end,2),'r.')


if show_animation
    figure(1); clf
    av = 2*pi*[0:.005:1];
    figure(1); clf
    
    plot(r*cos(av),r*sin(av),'k-','LineWidth',2); hold on; axis image;
    axis(1.2*[-1 1 -1 1]);
    theta = X0(1)-pi/2;
    p1 = plot(r*cos(theta),r*sin(theta),'ko','MarkerFaceColor','r',...
        'MarkerSize',12);
    dt = .05;
    max_t = 5;
    for t = dt:dt:max_t %max(tout);
        theta = interp1(tout,xout(:,1),t)-pi/2;
        %dtheta = interp1(tout,xout(:,2),t);
        set(p1,'XData',r*cos(theta),'YData',r*sin(theta));
        drawnow;
        pause(dt)
    end
    
    theta = xout(end,1)-pi/2;
    set(p1,'XData',r*cos(theta),'YData',r*sin(theta));
end

end

function dX = bead_ode(t,X)

global bfac

g = 9.8;
r = 1.0; % ring radius
m = 0.1; % mass of bead
w = 4; % angular velocity of hoop spinning
%bfac = 10;
%b = sqrt(bfac*m*m*g*r);  % we want "b^2" >> "m^2*g*r", for overdamped case
b = bfac;

%roots([m*r, b, m*g-m*r*w*w])

theta = X(1);
dtheta = X(2);

d2theta = (-b*dtheta - (m*g-m*r*w*w*cos(theta))*sin(theta))/(m*r);
dX = [dtheta; d2theta];
end


