function Control1_3_JMS(model)
% Solution to problem 3 of control 1. 
% JMS, Oct.2022

% Set data
if model=='A'
  R = 0.8;     % trajectory radius (m)
  T = 5;       % trajectory period (s)
elseif model=='B'
  R = 0.5;
  T = 10;
elseif model=='C'
  R = 0.75;
  T = 7;
elseif model=='D'
  R = 1.2;
  T = 3;
else
  error('unknown model')
end
tmax = 2*T;  % max. time (s)
dt = 0.001;  % time interval (s)

% Find velocity
w = 2*pi/T;
t = 0:dt:tmax;
vx = -R*w*sin(2*w*t).*(5-8*sin(w*t).^2);
vy = -R*w*(1-4*sin(2*w*t).^2+2*sin(w*t).^2);
v = sqrt(vx.^2+vy.^2);

% Plot velocity
figure(11)
plot(t,vx,t,vy,'LineWidth',1.5)
grid on
legend('v_x','v_y')
xlabel('time (s)')
ylabel('velocity (m/s)')
set(gca,'FontSize',15)

% Find position and travelled distance
x = R+integral_JMS(t,vx);
y = 0+integral_JMS(t,vy);
s = integral_JMS(t,v);
fprintf('total distance = %f m\n',s(end))

% Plot trajectory
figure(12)
plot(x,y,'LineWidth',1.5)
axis equal
grid on
xlabel('x (m)')
ylabel('y (m)')
set(gca,'FontSize',15)

% Find acceleration
ax = derivative_JMS(t,vx);
ay = derivative_JMS(t,vy);
a = sqrt(ax.^2+ay.^2);

% Plot acceleration
figure(13)
plot(t,ax,t,ay,t,a,'LineWidth',1.5)
grid on
legend('a_x','a_y','|a|')
xlabel('time (s)')
ylabel('acceleration (m/s^2)')
set(gca,'FontSize',15)

end
