%  Numerical integration of equations of
%  motion for particle in circular cone...

%  x(1)=r, x(2)=r_dot, x(3)=theta

global r1
r1 = 0.5;

t0 = 0;
tf = 6;
x0 = [r1 0 0]';

options = odeset('RelTol',1e-5,'AbsTol',1e-6);
[t,x] = ode45('cone_eqns', [t0,tf], x0);

%  Calculate auxiliary results and plot...

[rows,cols] = size(x);
for i=1:rows,
  vv = 4*x(i,2)*x(i,2) + 0.04/(x(i,1)^2);
  v(i) = sqrt(vv);
  h(i) = (r1-x(i,1))/tan(pi/6);
  energy(i) = 0.05*(vv) + 1.69914*(x(i,1)-r1);
end

figure(1)
plot(h,v)
xlabel ('h (m)')
ylabel ('v (m/s)')
title ('Speed versus Height Fallen')
%print -dps plot1

figure(2)
plot(t,energy)
xlabel ('t (s)')
ylabel ('energy (J)')
title ('Total Mechanical Energy versus Time')
%print -dps plot2

%  Plot of 3-d trajectory of particle on cone's surface.

figure(3)
xx = diag(x(:,1))*cos(x(:,3));
yy = diag(x(:,1))*sin(x(:,3));
zz = (x(:,1)-0.5)/tan(pi/6);
plot3(xx,yy,zz);
xlabel ('x (m)')
ylabel ('y (m)')
zlabel ('z (m)')
%print -dps plot3

%  For an animated trajectory, use:
%figure(4)
%comet3(xx,yy,zz);
%xlabel ('x (m)')
%ylabel ('y (m)')
%zlabel ('z (m)')