%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% solves the harmonic oscillator with damping   %
%                                               %
%   x' = y...                                   %
%   y' = -(k/m)*x - (c/m)*y                     %
%                                               % 
% requires the use of file:  harmosc.m          % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

global c;
global m;
global k;


  c = .3;
  m = 1;
  k = 1;

% time span of integration
  tspan = [0 10];
% initial conditions
x0 = 1;
y0 = 0;
v0 = [x0,y0];

% solves ode using ode45 routine

[t,vout] = ode45(@harmosc,tspan,v0);

% 
%   plots y(t) for xstar=0.0 in an upper subplot
%

% note:
%%paramstring = sprintf( ', {\\itb}=%d, {\\itc}=%d,{\\itx_0}=%d', b, c, xstar);

figure(1)
clf
subplot(3,1,1)

plot( t, vout(:,1) );
xlabel( '{\itt}' );
ylabel( '{\itx}({\itt})' );
%title('Prey')
axis tight;
grid on;

subplot(3,1,2)

plot( t, vout(:,2) );
xlabel( '{\itt}' );
ylabel( '{\ity}({\itt})' );
%title( 'Predator')
axis tight;
grid on;

subplot(3,1,3)

plot( vout(:,1), vout(:,2) );
xlabel( '{\itx}' );
ylabel( '{\ity}' );
title( 'Phase plane')
%axis tight;
axis equal;
grid on;

%figure(2)

hbead=line(vout(1,1),vout(1,2),'marker','o','markersize',8,'erase','xor');
%axis([-1 1 -1 1]);
%axis('square');
for k=2:length(t)
    set(hbead,'xdata',vout(k,1),'ydata',vout(k,2));
    drawnow
    pause(.1)
end