%Bifurcation diagram for the Furuta pendulum
%   kja 070816
%Parameters
aminit;
m=0.098;g=9.8;r=0.215;l=0.148;Jp=2.62e-3;J0=m*r^2
%m1=0.098;l1=0.215;r1=0.148;J1=4.72e-4;Jp=J1+m1*l1^2;
w0=sqrt(m*g*l/Jp);dw=w0/100;
%Bifurcation diagram
w=0:dw:20;one=ones(size(w));
ws=0:dw:w0;onestable=ones(size(ws));
wu=w0:dw:20;oneunstable=ones(size(wu));
subplot(221)
pl=plot(w,0*one,'b--');
set(pl,'Linewidth', AM_data_linewidth);
hold on
pl=plot(ws,onestable,'b-',ws,-onestable,'b-');
set(pl,'Linewidth', AM_data_linewidth);
pl=plot(wu,oneunstable,'b--',wu,-oneunstable,'b--');
set(pl,'Linewidth', AM_data_linewidth);
pl=plot(wu,-1+acos(m*g*l/Jp./wu.^2)/pi,'b',wu,1-acos(m*g*l/Jp./wu.^2)/pi,'b');
set(pl,'Linewidth', AM_data_linewidth);
amaxis([0 20 -1.1 1.1]);
xlabel('w');ylabel('th/pi');

amprint('furuta-bifurcation.eps');
