function Control1_1(model)
% Solution to problem 1 of control 1
% J.M.Soler, Oct.2018

% Set parameters
if model=='A'
  a = 1;        % radius of external circle (m)
  b = 2/5;      % radius of internal circle (m)
  w = pi/5;     % angular velocity (rad/s)
  tf = 20;      % final time (s)
elseif model=='B'
  a = 2;
  b = 1/5;
  w = pi/4;
  tf = 8;
elseif model=='C'
  a = 1;
  b = 2/7;
  w = pi/3;
  tf = 12;
elseif model=='D'
  a = 2;
  b = 3/4;
  w = pi/2;
  tf = 12;
else
  error('Control1_1 ERROR: bad model')
end
dt = 0.01;     % time interval (s)

% Find velocity
t = 0:dt:tf;
vx = (a-b)*w*( -cos(w*t) + cos((a-b)/b*w*t) );
vy = (a-b)*w*( -sin(w*t) - sin((a-b)/b*w*t) );
vt = sqrt(vx.^2+vy.^2);

% Find acceleration by numerical derivation
ax = derivative(t,vx);
ay = derivative(t,vy);
at = sqrt(ax.^2+ay.^2);

% Find trajectory by numerical integration
x = 0+integral(t,vx);
y = a+integral(t,vy);
l =   integral(t,vt);   % distance travelled

% Plot velocity
figure(1)
plot(t,vx, t,vy, t,vt)
legend('v_x','v_y','v_{tot}')
xlabel('t (s)')
ylabel('v (m/s)')

% Plot acceleration
figure(2)
plot(t,ax, t,ay, t,at)
legend('a_x','a_y','a_{tot}')
xlabel('t (s)')
ylabel('a (m/s^2)')

% Plot trajectory
figure(3)
plot(x,y)
xlabel('x (m)')
ylabel('y (m)')
axis equal
axis([-a,a,-a,a])

% Plot distance travelled
figure(4)
plot(t,l)
xlabel('t (s)')
ylabel('L (m)')

end
