%Molecular flywheel:
clear all
close all
% System parameters
r1=.4; r2=.8;alpha=.7*pi/2; 
sa = sin(alpha); ca = cos(alpha);
r12 = r1*r2;
% Initial conditions: put (2) in extreme position, (1) in proximal
y20 = asin(r2*ca/sa); c20 = cos(y20); s20 = sin(y20);
t = (r2*ca*s20-sa)/(r2*c20);
y10 = atan(t);
tspan = [0 12];                % solve for 0 <= t <= 50
yzero = [y10;y20];               % Initial conditions; defines phi0
 options = odeset('RelTol',1e-8,'AbsTol',[1e-12 1e-12]);
 [t,y] = ode45('flywheel',tspan,yzero,[],alpha,r1,r2);
 %[t,y] = ode45('flywheel',tspan,yzero,options);
c1 = cos(y(:,1)); s1 = sin(y(:,1));
c2 = cos(y(:,2)); s2 = sin(y(:,2));
k = sa*(r2*s2 - r1*s1) + r1*r2*(c1.*c2+s1.*s2*ca);
%E = .5*(y(:,2).^2 + y(:,4).^2) -k;
l1 = zeros(length(t)); l1(:) = -pi/2;
l2 = zeros(length(t)); l2(:) =  pi/2;
plot(t,y(:,1),t,y(:,2),t,k,t,l1,t,l2)
figure
plot(y(:,1),y(:,2))
xlabel('y_1','FontSize',14)
ylabel('y_2','FontSize',14)%,'Rotation',0,'Horizontal Alignment','right')
title('flywheel equations','FontSize',16)