% demonstrates curvature of an ellipse
% Author: Fernando Guevara Vasquez, Aug 29 2007

tt=linspace(0,2*pi,100);
% equation of the ellipse
xx=2*cos(tt); yy=5*sin(tt);

% loop over several points of the curve
for t=linspace(0,2*pi,41),
 px=2*cos(t); py=5*sin(t); % current point (px,py)
 kappa = 10/(4+21*cos(t)^2)^(3/2);
 % unit tangent vector
 tx = -2*sin(t)/sqrt(4+21*cos(t)^2); ty= 5*cos(t)/sqrt(4+21*cos(t)^2);
 % unit normal vector
 nx = -5*cos(t)/sqrt(4+21*cos(t)^2); ny=-2*sin(t)/sqrt(4+21*cos(t)^2);
 % center of the curvature circle
 xc = px + nx/kappa;
 yc = py + ny/kappa;
 % parametric equation of the circle centerd at (xc,yc) of rad. 1/kappa
 xxc=xc+(1/kappa)*cos(theta);
 yyc=yc+(1/kappa)*sin(theta);

 % plotting
 plot(xx,yy,'g',xxc,yyc,'r',...
      px,py,'+',xc,yc,'+',[px,xc],[py,yc],'b',[px,px+tx],[py,py+ty],'b');
 title(sprintf('curvature circle for P=(%g,%g)',px,py));
 axis([-5,5,-5,5])
 axis equal;
 pause(0.1);
end;