
r(t)=(2cos(t),sin(t))
a = 2; b = 1; t1 = 0; t2 = 2*pi;
r = [a*cos(t), b*sin(t)]; % ellipse
r1 = diff(r, t); r2 = diff(r1, t); d = sqrt(r1*r1');
T = simplify(r1/d); N = [-T(2), T(1)] % unit tangent and normal vectors
k = simplify((r1(1)*r2(2) - r1(2)*r2(1))/d^2) % the curvature k = 23sin2 t+1
Circ = r + N/k+ (N*cos(2*pi*s) + T*sin(2*pi*s))/k% the osculating circle
m = 18; 
for n = 0 : m; % prepare animation
ezplot(r(1),r(2),[t1,t2]); hold on; axis equal;
tn = t1 + n*(t2 - t1)/m; 
C = subs(Circ, t, tn);
rT = subs(r + s*T, t, tn); 
rN = subs(r + s*N, t, tn);
hn = ezplot(rN(1), rN(2), [0, 1]); set(hn, 'Color', 'r');
ht = ezplot(rT(1), rT(2), [0, 1]); set(ht, 'Color', 'r');
hc = ezplot(C(1), C(2), [0, 1]); set(hc, 'Color', 'g');
pause(0.1); hold off;
end;
