clear; close all
%
%  Response of a spring-mass-damper system
%  to initial displacement
%  Different formulas for three damping cases
%  Save for MECH 141
%
m = 100;                %% mass, kg
k = 1000;               %% spring constant, N/m
omega_n = sqrt(k/m);    %% undamped frequency, Hz
x0 = 0.01;              %% Initial displacement, m
period = 2*pi/omega_n;  %% vibration period
t = (0:0.01:2)*period;
leg = zeros(0,8);
for b = [100 400 700 1200]  %% damping constant, N-sec/m
    zeta = (b/2)/sqrt(k*m);
    switch sign(zeta-1)
        case -1         %% zeta < 1, under-critically damped
            sq = sqrt(1-zeta^2);
            phase = atan(zeta/sq);
            omega_d = omega_n*sq; %% damped freq
            eterm = exp(-zeta*omega_n*t);
            cterm = cos(omega_d*t - phase);
            x = (x0/sq) * eterm .* cterm;
            dtype = 'under-damped';
            color = 'b';
        case 0          %% zeta = 1, critically damped
            x = x0*(1 + omega_n*t).*exp(-omega_n*t);
            dtype = 'critically damped';
            color = 'g';
        case 1          %% zeta > 1, over-critically damped
            sq = sqrt(zeta^2-1);
            p = omega_n * (zeta+sq);
            term1 = ((-zeta+sq)/(2*sq))*exp(-p*t);
            p = omega_n * (zeta-sq);
            term2 = ((zeta+sq)/(2*sq))*exp(-p*t);
            x = x0*(term1 + term2);
            dtype = 'over-damped';
            color = 'r';
    end
    plot(t,x*1000,color);   %% plot mm
    hold on
    leg = [leg; sprintf('b=%6.1f',b)];
end
legend(leg)
grid
fmt = 'Response to initial displacement, k=%.1f N/m, m=%.1f kg';
title(sprintf(fmt,k,m));
xlabel('Time, sec');
ylabel('Displacement, mm');