% MITES08: Forced oscillations [forcedSHM.m]
% June 30, 08 [Mon.]
% Plots position of block subject to forced sinusoidal SHM.

clear;

% ***** USER: Change the parameters and explor *****
gamma = 0;  % gamma = b/(2m), where b = damping constant.
phi = 0; % phase constant.
C = 5; % Amplitude of SHM part of solution.
m =1; % mass of block
F0 = 300000; % Amplitude of force (Assuming sinusoidal force here).
omega0 = 200; % Natural angular frequency.
omega= 19;    % Angular frequency of force.
omegaenne = 0;    % REPLACE omega0 with omegaenne if gamma is not zero!



% ***** Change the following parameter for plotting ***** 
n = 3;  % For how many multiples of 2pi/omega do you want to plot the motion?

% After solving EOM for sinusoidally forced SHM (Pg. 44 of lecture notes),
% we obtain "forcedSHM" being x(t) = position of particle:
t = linspace(0, n*2*pi/omega, 100000); % Time: 0 to n*2pi/omega, divided into 1000 equal increments.
forcedSHO = C*exp(-gamma*t).*cos(omega0*t - phi) + (F0/(m*(omega0^2 - omega^2)))*cos(omega*t);

% Plot the position of block:
plot(t, forcedSHO, '--r');
hold on
xlabel('time');
ylabel('Position of forced SHO');
title('Motion of particle in sinusoidally forced SHO');
hold off



