function Control1_1(model)

% Exercise 1 of control 1, 30/10/2015
% Find and plot the trajectory and energy of a damped harmonic oscillator
% Written as a function to include all models of control 1.
% Input:
%   model : ('A'|'B'|'C'|'D')
% Sample usage:
%   Control1_1('B')
% J.M.Soler, Oct.2015

% Set system parameters
if model=='A'
    m = 1.4;  % system mass (kg)
    k = 6.5;  % spring constant (N/m)
    c = 0.8;  % damping constant (kg/s)
    x0 = 2.8; % initial position (m)
    tf = 18;  % final time (s)
elseif model=='B'
    m = 3.6;
    k = 10.1;
    c = 1.1;
    x0 = -6.2;
    tf = 11;
elseif model=='C'
    m = 5.1;
    k = 4.8;
    c = 0.3;
    x0 = -4.5;
    tf = 38;
elseif model=='D'
    m = 7.9;
    k = 3.4;
    c = 2.5;
    x0 = 3.9;
    tf = 24;
else
    error('Wrong model')
end
w0 = sqrt(k/m);       % oscillator frequency (rad/s)
dt = (2*pi/w0)/100;   % time interval (s) (100 times per period)

% Find oscillator positions
t = (0:dt:tf); % time mesh
x = dho_position( m, c, k, x0, t );

% Plot trajectory (exercise 1B)
figure(11)
subplot(3,1,1)
plot(t,x)
xlabel('time (s)')
ylabel('position (m)')

% Find and plot velocity (exercise 1C)
v = derivative( t, x );
subplot(3,1,2)
plot(t,v)
xlabel('time (s)')
ylabel('velocity (m/s)')

% Find and plot acceleration (exercise 1D)
a = derivative( t, v );
subplot(3,1,3)
plot(t,a)
xlabel('time (s)')
ylabel('acceleration (m/s^2)')

% Find work of damping force (exercise 1E)
W = integral( t, -c*v.^2 );

% Find and plot energies
T = 0.5*m*v.^2;
U = 0.5*k*x.^2;
E = T+U;
figure(12)
plot( t,T, t,U, t,E, t(1:10:end),E(1)+W(1:10:end),'ok' )
legend( 'Kinetic','Potential','Total','Initial-dampingWork')
xlabel('time (s)')
ylabel('energy (J)')
