% demonstrates curvature of a parabola
% Author: Fernando Guevara Vasquez, Aug 29 2007

% the parabola
xx=linspace(-4,4,100);
yy=xx.^2/2;

theta=linspace(0,2*pi,100);

% loop over several points of the curve
for t=linspace(-2,2,50),
 px=t; py=t^2/2; % current point (px,py)
 kappa = (1+t^2)^(-3/2); % curvature at point 
 % unit normal vector
 nx = -t/sqrt(1+t^2); ny =  1/sqrt(1+t^2);
 % unit tangent vector
 tx = 1/sqrt(1+t^2);  ty = t/sqrt(1+t^2);
 % center of 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([-4,4,0,8])
 axis equal;
 pause(0.1);
end;
