% ### HOode45EX.m ###       09.19.14
% Numerically integrate the damped/driven harmonic oscillator
%   m*x''+ b*x' + k*x = A*sin(wt)

clear
% -----------------------------------------------------
% User input (Note: All paramters are stored in a structure)
P.y0(1) = 0.0;   % initial position [m]
P.y0(2) = 1.0;   % initial velocity [m/s]
P.b= 0.1;  % damping coefficient [kg/s]
P.k= 250.0;   % stiffness [N m]
P.m= 0.01;    % mass [kg]

% sinusoidal driving term
P.A= 0.0;   % amplitude [N] (set to zero to turn off)
fD= 1.05*sqrt(P.k/P.m)/(2*pi);  % freq. (Hz) [expressed as fraction of resonant freq.]

% Integration limits
P.t0 = 0.0;   % Start value
P.tf = 3.0;   % Finish value
P.dt = 0.0001;  % time step
% ----------------------------------------------------------------------
% +++
% spit back out some basic derived quantities 
P.wr= 2*pi*fD;  % convert to angular freq.
disp(sprintf('Resonant frequency ~%g [Hz]', sqrt(P.k/P.m)/(2*pi)));
Q = (sqrt(P.k/P.m))/(P.b/P.m);  % quality factor
disp(sprintf('Q-value = %g', Q));
% +++
% use built-in ode45 to solve
[t y] = ode45('HOfunction', [P.t0:P.dt:P.tf],P.y0,[],P);

% ------------------------------------------------------
% visualize
figure(1); clf;
plot(t,y(:,1)); hold on; grid on; 
xlabel('t [s]');    ylabel('x(t) [m]')
% Phase plane 
figure(2); clf;
plot(y(:,1), y(:,2)); hold on; grid on; 
xlabel('x [m]');    ylabel('dx/dt [m/s]')