
clear 
close all
disp(' ')
disp(' graph of curve t -> (t^2 cos(t), t^2 sin(t), (2-t))  ')
disp(' on the interval 0 to 2. ')

  x = inline('t.^2.*cos(t)');
  y = inline('t.^2.*sin(t)');
  z = inline('2-t');
t = linspace(0, 2, 101);
plot3(x(t),y(t),z(t))
hold on
axis([-3,1, -.5, 4, 0, 3])
axis off
origin = [0 0 0];
arrow3(origin, [1,0 0],'k');
arrow3(origin, [0,1,0], 'k');
arrow3(origin, [0,0,3],'k');
text(1.2,0,0, '\fontsize{14}x')
text(0,1,0, '\fontsize{14} y')
text(.4,0,3, '\fontsize{14}z')
text(0,0,2, '(0,0,2)')
view(123, 30) 
rotate3d on
disp(' ')
disp(' We add the osculating circle at several points . ') 
disp(' the curvature is greatest at t = 0 when kappa = 2 ')
disp(' and then decreases. The osculating circle is in the xz')
disp(' plane when t = 0')
pause
tau = linspace(0, 2*pi, 101);

ss = [0,.2 .4 .6 .8 1 1.2 1.4 1.6];
disp('   ')
    format compact
   for j = 1:9

   s = ss(j);
   p = [x(s), y(s), z(s)];
   v = [2*s*cos(s)-s^2*sin(s), 2*s*sin(s)+s^2*cos(s), -1];
   a = [2*cos(s)-4*s*sin(s)-s^2*cos(s), 2*sin(s)+4*s*cos(s)-s^2*sin(s), 0];
   T = v/norm(v);
   num = dot(v,v)*a - dot(a,v)*v;
   N = num/norm(num);
   hN = arrow3(p,N, 'k')
   hT = arrow3(p,T, 'r')
   kappa = norm(cross(a,v))/norm(v)^3;
   rho = 1/kappa;
   disp(' ')
   disp(' ')
   disp('     kappa    rho')
   [kappa, rho]
   disp('      T        N')
   [T', N']
    center = p+ rho*N;
   xcircle = rho*cos(tau)*T(1) + rho*sin(tau)*N(1) + center(1);
   ycircle = rho*cos(tau)*T(2) + rho*sin(tau)*N(2) + center(2);
   zcircle = rho*cos(tau)*T(3) + rho*sin(tau)*N(3) + center(3);
   hcircle = plot3(xcircle,ycircle,zcircle,'r');
   hpoint = plot3(x(s), y(s),z(s), '*');
   rotate3d on
   pause
   delete(hT); delete(hN); delete(hcircle); delete(hpoint);
end
format loose
hold off
pause
close
