clear 
close all
disp(' ')
disp(' ARC LENGTH ')
disp(' ')
disp(' First we show how to make polygonal approximations to ')
disp(' a curve and make an approximation to the arc length  ')
disp(' We consider the curve in R^2, t -> (t*cos(t), t*sin(t))')
disp(' on the interval 0 to 2*pi  ')

x = inline(' t.*cos(t)');
y = inline(' t.*sin(t)');
clear t
t = linspace(0, 2*pi, 101);
plot(x(t),y(t))
hold on
disp(' ')
disp(' Three different approximations are made: first with 5 segments,')
disp(' then with 10 segments and finally with 20 segments. ')
disp(' To make the 5 segment approximation, enter return 5 times ')
disp(' The other approximations are made automatically ')


l = zeros(3,1);

clear s
s = linspace(0, 2*pi, 6);
length = 0;

for n = 1:5
    pause
    p = [x(s(n)), y(s(n))];
    q = [x(s(n+1)), y(s(n+1))];
    length = length + norm(p-q)
    plot([x(s(n)), x(s(n+1))], [y(s(n)), y(s(n+1))], 'r')
end
l(1) = length;


disp('Now we do the 10 segment approximation  ')
pause
clf
plot(x(t),y(t))
hold on
clear s
s = linspace(0, 2*pi, 11);
length = 0;
for n = 1:10
    p = [x(s(n)), y(s(n))];
    q = [x(s(n+1)), y(s(n+1))];
    length = length +norm(p-q);
    plot([x(s(n)), x(s(n+1))], [y(s(n)), y(s(n+1))], 'r')
end
length
l(2) = length;


disp(' ')
disp(' Next we do the 20 segment approximation  ')
pause

clf
plot(x(t),y(t))
hold on
clear s
s = linspace( 0, 2*pi, 21);
length = 0;
for  n = 1:20
   p = [x(s(n)), y(s(n))];
   q = [x(s(n+1)), y(s(n+1))];
   length = length +norm(p-q);
   plot([x(s(n)), x(s(n+1))], [y(s(n)), y(s(n+1))], 'r')
end
length
l(3) = length;
disp(' ')
disp(' The results of the polygonal line approximation are ')
disp(' ')
disp('  number of segments     length of polygonal line ')
N = [5;10;20];
  [N, l]


pause
disp(' the approximate length using a numerical integrator  ')
disp(' on the speed function , speed = sqrt(1+t^2), is  ')
speed = inline('sqrt(1+t.^2)');
quadl(speed, 0, 2*pi)


pause
disp(' ')
disp(' In this case, the integral can be done analytically ')
disp(' with an exact answer of  ')
disp(' ')
disp(' pi*(1+4*pi^2)^(1/2)-1/2*log(-2*pi+(1+4*pi^2)^(1/2)) ')
disp( '  =  21.2563  ')

pause
disp(' ')
disp(' Another example, this time in 3D. Take the curve ')
disp(' t -> (t^2*cos(t), t^2*sin(t), t)  on the interval ')
disp('t from 0 to 2  ')
clear s t; clf
x = inline('t.^2 .*cos(t)');
y = inline('t.^2.*sin(t)');
z = inline('t');
t = linspace(0, 2, 101); 
plot3(x(t), y(t), z(t))
hold on
axis equal
disp(' ')
disp('Again we make a polygonal line approximation with 5 segments')
disp(' You must hit return 5 times ')

N = 5;
ss = linspace(0, 2, N+1);
length = 0;
for n = 1:N
    pause
   s1 = ss(n)
   s2 = ss(n+1)
   p = [x(s1), y(s1), z(s1)];
   q = [x(s2), y(s2), z(s2)];
   length = length + norm(p-q)
   plot3([x(s1), x(s2)], [y(s1), y(s2)], [z(s1),  z(s2)], 'r')
end
pause
disp(' ')
disp(' Now we make a 20 segment approximation ')
clf
plot3(x(t),y(t),z(t))
hold on
N = 20;
ss = linspace(0, 2, N+1);
length = 0;
for n = 1:N
   s1 = ss(n);
   s2 = ss(n+1);
   p = [x(s1), y(s1), z(s1)];
   q = [x(s2), y(s2), z(s2)];
   length = length + norm(p-q);
   plot3([x(s1), x(s2)], [y(s1), y(s2)], [z(s1),  z(s2)], 'r')
end
length 

disp(' Now we use a numerical integrator to integrate the speed ')
disp(' function, speed = sqrt( 1 +4t^2 + t^4)  ')
speed = inline('sqrt(t.^4 + 4*t.^2 + 1)');
quadl(speed, 0,2)

