% ### EXharmOSCode45.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] {0}
P.y0(2) = 1.0;   % initial velocity [m/s] {1}
P.b= 0.05;  % damping coefficient [kg/s] {0.05}
P.k= 250.0;   % stiffness [N m] {250}
P.m= 0.01;    % mass [kg] {0.01}

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

% Integration limits
P.t0 = 0.0;   % Start value {0}
P.tf = 3.0;   % Finish value {3}
P.dt = 0.0001;  % time step {0.0001}
% ----------------------------------------------------------------------
% +++
% 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]')

% ------------------------------------------------------
% compute Hilbert transform to extract inst. mag. and phase
if 1==0
    wH= hilbert(y(:,1));
    figure(3); clf;
    subplot(211); hold on; grid on; title('Analytic signal representation')
    plot(t,abs(wH)); xlabel('t [s]'); ylabel('Magnitude')
    subplot(212); hold on; grid on; title('Analytic signal representation')
    ang= cycs(wH)- fD*t;    % subtract accumulating phase settling into driving freq. (better way to do this?)
    plot(t,ang); xlabel('t [s]'); ylabel('Phase [cyc]')
    
    %theta= angle(wH)-angle(exp(-i*P.wr*t));  
    %plot(t,(theta)); xlabel('t [s]'); ylabel('Phase [cyc]')
end