Simulación sistemas rotacionales (Rueda)

function rotacional
R=1;

% declara variables

time=0:1/10:30;
condicionesiniciales=[0 10];

% desarrolla las ecuaciones diff rk 4 orden
[t,theta]=ode45(@fk,time,condicionesiniciales)
x=R*sin(theta(:,1));
y=-R*cos(theta(:,1));

% Animación graficas

for i=1:length(time)

% plot velocidad angular
pause(1/100)
subplot(2,2,1)
plot(time(i),theta(i,2),'-')
xlabel('tiempo (s)')
ylabel('w (rad/s)')
hold on

% plot rotacion angular
subplot(2,2,2)
plot(time(i),theta(i,1),'-')
xlabel('tiempo (s)')
ylabel('Angulo (rad)')
hold on

% grafico circulo
xc=0; yc=0; r=1; % Centro y radio
n = 50; k=0:n; fi=2*pi*k/n;
xx=xc+r*cos(fi); yy = yc+r*sin(fi);
subplot(2,2,3)
plot(xc,yc,'x',xx,yy);
axis([-2 2 -2 2]), axis equal;
hold on

%plot linea para visualizar la velocidad angular
subplot(2,2,3)
plot([0 x(i)],[0 y(i)])
axis([-2 2 -2 2])
xlabel('x(m)')
ylabel('y(m)')
hold off

end

end

%Ecuacion diferencial

function rk=fk(t,y)
m=2;
b=0.1;
R=1;
rk=[y(2);-(b/(m*R))*y(2)];
end