%Molecular flywheel:
function [t1, t2,ttt] = wheels(r1,r2,alpha,tspan)% System parameters
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);
th = (r2*ca*s20-sa)/(r2*c20);
y10 = atan(th);                   % 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);
t1 = y(:,1)'; t2 = y(:,2)';ttt = t';
return

