function Control1_2_JMS(version)

% Control1_2.m
% Solution to problem 2 of control 1
% JMS, Oct.2019

% Set parameters
m = 0.3;        % particle mass (Kg)
r0 = 0.12;      % radial parameter (m)
if strcmp(version,'A')
  w = 3;        % angular velocity (rad/s)
  gamma = pi/9; % angular phase (rad)
  tf = 10;      % final time (s)
elseif strcmp(version,'B')
  w = 2;
  gamma = -pi/3;
  tf = 12;
elseif strcmp(version,'C')
  w = 8;
  gamma = pi/6;
  tf = 3;
elseif strcmp(version,'D')
  w = 5;
  gamma = -pi/12;
  tf = 6;
elseif strcmp(version,'Z')
  w = 6;
  gamma = pi/8;
  tf = 5;
end
dt = 1e-3;      % time interval (s)

% Find and plot trajectory
t = (0:dt:tf);
alpha = 1-exp(-w*t/5);
x = (4*r0/3)*alpha.*cos(w*t+gamma);
y = (3*r0/4)*alpha.*sin(w*t);
r = [x;y];
figure(11)
plot(x,y)
axis equal
axis([-0.2,0.2,-0.15,0.15])
grid on
xlabel('x (m)')
ylabel('y (m)')
title('Position')

% Find and plot velocity and acceleration
vx = derivative(t,x);
vy = derivative(t,y);
ax = derivative(t,vx);
ay = derivative(t,vy);
v = [vx;vy];
a = [ax;ay];
f = m*a;
[ft,fn] = forces_tan_nor_JMS(v,f);
vmod = sqrt(sum(v.^2));
fmod = sqrt(sum(f.^2));
ftmod = sqrt(sum(ft.^2));
fnmod = sqrt(sum(fn.^2));
figure(12)
plot(t,ftmod, t,fnmod)
legend('|F_t|','|F_n|','Location','Northwest')
grid on
xlabel('t (s)')
ylabel('forces (N)')
title('Force')

% Find kinetic energy and work of force
Ekin = 0.5*m*vmod.^2;
W = integral(t,sum(f.*v));
figure(13)
plot(t,Ekin(1)+W,'r-', t,Ekin,'k--')
legend('Work','E_{kin}','Location','Northwest')
grid on
xlabel('t (s)')
ylabel('Energy (J)')
title('Energy')

% Find fraction of time with y<r0
fraction_r_LT_r0 = numel(find(y<=r0/2))/numel(t)

end
