% This code shows the position of a damped oscillator as a function of
% time.

% First Piece

set(0,'DefaultAxesFontSize',24);
set(0,'DefaultLineLineWidth',3);

A=30; % (Amplitude)
omega=0.5;
beta1=0; % No damping case.
beta2=0.01; %
%beta3=0.02; % .
beta3=0.25; 
omega1_1=(omega.^2-beta1.^2).^(1/2);
omega1_2=(omega.^2-beta2.^2).^(1/2);
omega1_3=(omega.^2-beta3.^2).^(1/2);

t=[0:.1:300];
x1=A.*exp(-beta1.*t).*(cos(omega1_1.*t));
x2=A.*exp(-beta2.*t).*(cos(omega1_2.*t));
x3=A.*exp(-beta3.*t).*(cos(omega1_3.*t));

figure
subplot(3,1,1)
hold on
plot(t,A.*exp(-beta1.*t),'r')
plot(t,-A.*exp(-beta1.*t),'r')
plot(t,x1,'k')
ylabel('Amplitude')
title(strcat('\beta =',sprintf('%0.2f',beta1),'      \omega_\circ =',sprintf('%0.2f',omega)))
subplot(3,1,2)
plot(t,x2,'k')
hold on
plot(t,A.*exp(-beta2.*t),'r')
plot(t,-A.*exp(-beta2.*t),'r')
ylabel('Amplitude')
title(strcat('\beta =',sprintf('%0.2f',beta2),'      \omega_\circ =',sprintf('%0.2f',omega)))
subplot(3,1,3)
plot(t,x3,'k')
hold on
plot(t,A.*exp(-beta3.*t),'r')
plot(t,-A.*exp(-beta3.*t),'r')
ylabel('Amplitude')
title(strcat('\beta =',sprintf('%0.2f',beta3),'      \omega_\circ =',sprintf('%0.2f',omega)))
xlabel('Time (a.u.)')

% Second Piece
figure
plot(t(1:1000),x1(1:1000),'k')
hold on
plot(t(1:1000),x2(1:1000),'r')
plot(t(1:1000),x3(1:1000),'m')
xlabel('Time (a.u.)')
legend('No Damping','Weak Damping','Less Weak Damping')