%Uses the files from Matlab Central File Exchange "plot2svg"
clear all
t = 0:0.01:11;
damping = [0.5, 1, 1.5, 0];  % damping values
colors = {'b', [0, 0.7, 0], 'r', 'm'};
h_fig = figure;
set(h_fig, 'Position',[360 673 354 249])
h_axes = axes;
axis([0, t(end), -1.05, 1.05]);
xlabel('time(s)')
ylabel('x/x0')
hold on
omega_0 = 1;  % natural frequency
for k = 1:length(damping)
  gamma = damping(k) * omega_0;
  if damping(k) > 1  % over-damped
    gamma1 = gamma - sqrt(gamma^2 - omega_0^2);
    gamma2 = gamma + sqrt(gamma^2 - omega_0^2);
    x = (gamma2 * exp(-gamma1 * t) - gamma1 * exp(-gamma2 * t)) / (gamma2 - gamma1);
  else
    if damping(k) == 1  % critically damped
      x = (1 + gamma * t) .* exp(-gamma * t);
    else  % under-damped
      omega_prime = sqrt(omega_0^2 - gamma^2);
      alpha = atan(gamma / omega_prime);
      x = exp(-gamma*t) .* cos(omega_prime * t - alpha) / cos(alpha);
    end
  end
  plot(t, x, 'Color', colors{k}, 'LineWidth', 1)
end
plot([0; t(end)], [0; 0], 'k')
legend('\zeta < 1', '\zeta = 1', '\zeta > 1', '\zeta = 0')
hold off
plot2svg('damping.svg',h_fig)