% Creates "return map" for the rimless wheel with point mass at "hip".
% Katie Byl, 1/8/2009

gamma = 7*pi/180;  % slope downhill
g=9.8; 
L=1; 
nspokes = 8;
alpha = (1/2)*(2*pi/nspokes); % alpha is 1/2 the interleg angle

omega_up_crit = -sqrt(2*g/L*(1-cos(alpha+gamma)));
omega_dn_crit = sqrt(2*g/L*(1-cos(alpha-gamma)));
omega_uphill = omega_up_crit+[-20:.01:0];
omega_middle = omega_up_crit + [0:.01:1]*(omega_dn_crit-omega_up_crit);
omega_dnhill = omega_dn_crit+[0:.01:6];

% RETURN MAP, for each of 3 cases:
omega_dnhill_2 = cos(2*alpha)*sqrt(omega_dnhill.^2 + 4*g/L*sin(alpha)*sin(gamma));
fixed_rolling = interp1(omega_dnhill_2-omega_dnhill,omega_dnhill,0)  % rad/s
omega_uphill_2 = -cos(2*alpha)*sqrt(omega_uphill.^2 - 4*g/L*sin(alpha)*sin(gamma));
omega_middle_2 = -cos(2*alpha)*omega_middle;
fixed_stopped = interp1(omega_middle_2-omega_middle,omega_middle,0)  % rad/s

% Now, just plot everything:
figure(10); clf
fs = 12; % Font Size
plot(omega_uphill, omega_uphill_2,'b-'); hold on; grid on
plot(omega_dnhill, omega_dnhill_2,'b-');
plot(omega_middle, omega_middle_2,'b-');
plot([min(omega_uphill) max(omega_dnhill)],...
    [min(omega_uphill) max(omega_dnhill)],'k:')
plot([min(omega_uphill) max(omega_dnhill)],...
    -[min(omega_uphill) max(omega_dnhill)],'k:')
plot(omega_up_crit+[0 0],[-10 8],'r--');
plot(omega_dn_crit+[0 0],[-10 8],'r--');
plot(fixed_rolling,fixed_rolling,'ko','MarkerFaceColor',[ 0 .7 0]);
plot(fixed_stopped,fixed_stopped,'ko','MarkerFaceColor',[0 0 0]);
xlabel('\Omega_{n}^+','FontSize',fs)
ylabel('\Omega_{n+1}^+','FontSize',fs)
%title('Return map for the rimless wheel (8 degree slope)',...
%    'FontSize',fs)
axis equal
axis([-8 4 -8 4])
set(gca,'XTick',[-8:5],'YTick',[-8:5])

