% damposcseries.m
% 
% Function to plot x(t) for a damped oscillator
%    for a series of damping values
%    Takes initial conditions so that initial velocity is (nearly) zero 
% See also damposc.m
%
% RP 16 Oct. 2006; slight modifications Oct. 19, 2007

clear all
close all

t = 0.0:0.1:20;   % array of time values
A = 1.0;  % "amplitude," in zero damping limit

% gamma over 2 omega_0 -- make an array with logarithmically spaced
% increments, using the "logspace" command.  Could also do:
% logg = -1.0:0.05:1.0;  gw = 10.^logg;
gw = logspace(-1,1,40);

% Assign omega_0 = 1.0;

N = length(gw);

x = zeros(length(t),N);  % initialize

% Create the figure window
figure;
axis([min(t(:)) max(t(:)) -1 1]); hold on
xlabel('time', 'FontWeight', 'bold');
ylabel('x', 'FontWeight', 'bold');
title('Damped Oscillator -- varying \gamma', 'FontWeight', 'bold');
text(10,-0.8, '\gamma / 2 \omega_0 = ', 'FontWeight', 'bold');
text(4,0.9, 'Initial amplitude A, initial velocity 0; \omega_0 = 1', ...
    'FontWeight', 'bold');

            
for j=1:N,
    if (gw(j) < 1.0)
        % underdamped:  g < 2*w0
        w = sqrt(1.0 - gw(j)*gw(j));  % damped angular frequency
        alpha = atan(-1.0*gw(j));
        Ampl = A / cos(alpha);
        x(:,j) = Ampl*exp(-gw(j)*t).*cos(w*t + alpha);
        % make  previous plot gray
        if (j > 1)
          plot(t, x(:,j-1), 'color', [0.8 0.8 1.0], 'linewidth', 2.0);
          text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]);
        end
        plot(t, x(:,j), 'color', [0 0 0.5], 'linewidth', 2.0);
        text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold');
        if (j==1)
            % pause, to explain...
           text(14, -0.4, 'paused...', 'FontWeight', 'bold')
           pause
           text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0])
        end
        lastj = j;
    elseif (abs(gw(j)-1.0) < 0.001)
        % critically damped:  g = 2*w0
        if (j == lastj + 1)
            text(14, -0.4, 'paused...', 'FontWeight', 'bold')
            pause
            text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0])
            text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]);
        end
        % for zero inital velocity, B = A*omega_0
        x(:,j) = (A + A*t).*exp(-1.0*t);
        % make  previous plot gray
          plot(t, x(:,lastj), 'color', [0.8 0.8 1.0], 'linewidth', 2.0);
          text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]);
        plot(t, x(:,j), 'k', 'linewidth', 3.0);
        text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold');
        lastj = j;
    else
        % overdamped:  g > 2*w0
        if (j == lastj + 1)
            text(14, -0.4, 'paused...', 'FontWeight', 'bold')
            pause
            text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0])
            text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]);
        end
        beta = sqrt(gw(j)*gw(j) - 1.0);
        % for zero inital velocity, A1 = -A2*(gamma/2 - beta)/(gamma/2 + beta)'
        A2 = A / (1.0 - ((gw(j) - beta)/(gw(j) + beta)));
        A1 = -A2*(gw(j) - beta)/(gw(j) + beta);
        x(:,j) = A1*exp(-(gw(j) + beta)*t) + A2*exp(-(gw(j) - beta)*t);
        % make  previous plot gray
        if (j > (lastj+1))
          plot(t, x(:,j-1), 'color', [1.0 0.8 0.8], 'linewidth', 2.0);
          text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]);
        end
        plot(t, x(:,j), 'r', 'linewidth', 2.0);
        text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold');
    end
    pause(0.5);
end