
  clear; close
  disp(' ')
  disp('       Curvature and Osculating circle for the curve ')
  disp(' ')
  disp('          r(t)  =  (cos(t)^3, sin(t)^3, t) ')
  disp('             for t in the interval [0,pi] ')
  disp(' At t = 0 and at pi, the osculating circle lies in the yz plane.')
  disp(' At t = pi/2, the osculating circle lies in the xz plane. ') 

  t = linspace(0, pi, 101);
  x = inline(' cos(t).^3');
  y = inline('sin(t).^3');
  z = inline('t');

   plot3(x(t),y(t),z(t))
   hold on
   or = [0 0 0];
   arrow3([-1.5,0,0], [3,0,0],'k');
   arrow3(or, [0,1.5,0], 'k');
   arrow3(or, [0,0,4], 'k');
   view(130, 30)
   axis image
   text(1,.2,0, '\fontsize{14}x')
   text(0,1.2,0, '\fontsize{14}y')
   text(0,0,3.2,'\fontsize{14}z')
   axis off
   rotate3d on
   pause
   tau = linspace(0, 2*pi, 101);
   ss= linspace(0, pi, 21);
   format compact
   for j = 1:21
      s = ss(j);
      point = [x(s),y(s),z(s)]; 
      hpoint = plot3(x(s),y(s),z(s), '*');
      v = [-3*cos(s).^2*sin(s), 3*sin(s)^2*cos(s),1];
      T = v/norm(v);
      a = [ 6*cos(s)*sin(s)^2-3*cos(s)^3, 6*sin(s)*cos(s)^2-3*sin(s)^3,0];
      num = cross(a,v);
      N = num/norm(num);
      hT = arrow3(point, T, 'r');
      hN = arrow3(point, N, 'k');
      kappa = norm(cross(a,v))/norm(v)^3;
      rho = 1/kappa;
      disp(' ')
      disp('   t      kappa      rho   ')
      [s,kappa, rho ]
      disp(' ')
      disp('      T        N  ')
      [T', N']
      center = point + rho*N;
      xcircle = center(1) +rho*cos(tau)*T(1) + rho*sin(tau)*N(1);
      ycircle = center(2) +rho*cos(tau)*T(2) + rho*sin(tau)*N(2);
      zcircle = center(3) +rho*cos(tau)*T(3) + rho*sin(tau)*N(3);
      hcircle = plot3(xcircle, ycircle, zcircle, 'r');
      rotate3d on
      pause
      delete(hpoint); delete(hT); delete(hN); delete(hcircle);
      
   end   

   hold off
   format loose
   close
