function slide=splexmpl
%SPLEXMPL This is a slideshow file for use with playshow.m and makeshow.m
%   To see it run, type 'playshow splexmpl', 

%   Carl de Boor 01 dec 97
%   Copyright (c) 1987-2000 by C. de Boor and The MathWorks, Inc.
% $Revision: 1.6 $

if nargout<1,
   playshow splexmpl
else
  %========== Slide 1 ==========

  slide(1).code={
   'x = 0:4; y=-2:2;',
   's2 = 1/sqrt(2);',
   'clear v',
   'v(3,:,:) = [0 1 s2 0 -s2 -1 0].''*[1 1 1 1 1];',
   'v(2,:,:) = [1 0 s2 1 s2 0 -1].''*[0 1 0 -1 0];',
   'v(1,:,:) = [1 0 s2 1 s2 0 -1].''*[1 0 -1 0 1];',
   'sph = csape({x,y},v,{''clamped'',''periodic''});',
   'fnplt(sph), axis equal, axis off',
   '' };
  slide(1).text={
   '',
   '                     Here are some simple tasks for the spline toolbox.',
   '',
   '',
   ''};

  %========== Slide 2 ==========

  slide(2).code={
   'x = 2*pi*[0 1 .1:.2:.9]; y = cos(x); cs = csapi(x,y);',
   'fnplt(cs,2); xxyy =[-1 7 -1.2 1.2]; axis(xxyy)',
   'hold on, plot(x,y,''o'')' };
  slide(2).text={
   '                                                     Interpolation',
   '',
   '>> x = 2*pi*[ 0 1 .1:.2:.9 ]; y = cos(x);',
   '>> cs = csapi( x, y );',
   '',
   'constructs a cubic spline which matches cosine at the given points.',
   '',
   'To view it, use',
   '>> fnplt( cs )'};

  %========== Slide 3 ==========

  slide(3).code={
   'fnval(fnder(cs),[0 2*pi])*[-1;1];' };
  slide(3).text={
   'The cosine is 2pi-periodic. How well does our interpolant in CS do in that regard?',
   '',
   'We compute the difference in the first derivative at the two endpoints:',
   '',
   '>> fnval( fnder(cs), [0 2*pi] ) * [ -1; 1 ]',
   '',
   'ans =  -0.1375',
   ''};

  %========== Slide 4 ==========

  slide(4).code={
   'csp = csape(x,y,''periodic'');',
   'fnval(fnder(csp),[0 2*pi])*[-1;1];',
   'fnval(fnder(csp,2),[0 2*pi])*[-1;1];' };
  slide(4).text={
   'To enforce periodicity, use CSAPE instead of CSAPI:',
   '>> csp = csape( x, y, ''periodic'' );',
   '',
   ' Now, the check gives',
   '>>fnval( fnder(csp), [0 2*pi] ) * [ -1; 1 ]',
   ' ans =  3.3307e-016',
   '',
   'Even the second derivative now matches at the endpoints:',
   '>>fnval( fnder( csp, 2 ), [0 2*pi] ) * [ -1; 1 ]',
   ' ans =  1.3323e-015',
   ''};

  %========== Slide 5 ==========

  slide(5).code={
   ' pl = spapi(2,x,y);',
   ' hold on, fnplt(pl,''r'',2), hold off'};
  slide(5).text={
   'The   piecewise linear interpolant   to the same data is available via',
   '',
   '>> pl = spapi (2, x, y );',
   '',
   'Here we add it to the previous plot, in red:',
   '',
   '>> hold on; fnplt (pl, ''r'', 2 ), hold off'};

  %========== Slide 6 ==========

  slide(6).code={
   'x = linspace(0,2*pi,51); noisy_y = cos(x) + .2*(rand(size(x))-.5);',
   'fnplt( csapi(x,noisy_y) ), axis(xxyy), hold on, plot(x,noisy_y,''x'')',
   'tol = (.05)^2*(2*pi); fnplt( spaps(x,noisy_y,  tol ), ''r'', 2 ),  hold off' };
  slide(6).text={
   '                                                       Smoothing',
   'If the data are noisy, one would approximate rather than interpolate. For example,',
   'with',
   '>> x = linspace(0,2*pi,51); noisy_y = cos(x) + .2*(rand(size(x))-.5);',
   'interpolation would give the wiggly interpolant marked in blue',
   '>> fnplt( csapi( x, noisy_y ) )',
   'while',
   ' >> tol = (.05)^2*(2*pi); hold on,  fnplt( spaps( x, noisy_y,  tol ), ''r'', 2 ),  hold off',
   'gives the smoothed approximation, shown here in red.'};

  %========== Slide 7 ==========

  slide(7).code={
   'noisy_y([1 end]) = mean( noisy_y([1 end]) );',
   'lx = length(x); lx2 = round(lx/2); range = [lx2:lx 2:lx 2:lx2];',
   'sps = spaps([x(lx2:lx)-2*pi x(2:lx) x(2:lx2)+2*pi],noisy_y(range),2*tol);',
   'hold on, fnplt( sps, [0 2*pi], ''k'', 2), hold off' };
  slide(7).text={
   'The approximation is much worse near the ends of the interval, and is far',
   'from periodic. To enforce periodicity, approximate to periodically extended',
   'data, then restrict approximation to the original interval:',
   '>> noisy_y([1 end]) = mean(noisy_y([1 end]));',
   '>>lx = length(x);  lx2 = round(lx/2); range = [lx2:lx 2:lx 2:lx/2];',
   '>> sps = spaps( [x(lx2:lx)-2*pi x(2:lx) x(2:lx2)+2*pi], noisy_y(range), 2*tol);',
   '>> fnplt( sps, [0 2*pi], ''k'', 2)',
   'giving the more nearly periodical approximation, shown in black.'};

  %========== Slide 8 ==========

  slide(8).code={
   'spl2 = spap2(4, 4, x, noisy_y);',
   'fnplt(spl2,''b'',2); axis(xxyy), hold on, plot(x,noisy_y,''x''), hold off',
   '' };
  slide(8).text={
   '                                        Least-squares approximation',
   '',
   'Alternatively, one could use least-squares approximation to the noisy data by a',
   'spline with few degrees of freedom.',
   '',
   'For example, one might try a cubic spline with just four pieces:',
   '>> spl2 = spap2( 4, 4, x, noisy_y );',
   '>> fnplt( spl2, ''b'', 2 )'};

  %========== Slide 9 ==========

  slide(9).code={
   ' k = 5;   % i.e., we are working with quartic splines',
   ' x = 2*pi*sort([0 1 rand(1,10)]); y = cos(x);',
   'sp = spapi( optknt(x,k), x, y ); fnplt(sp,2,''g'');',
   ' hold on; plot(x,y,''o''); hold off, axis([-1 7 -1.1 1.1])' };
  slide(9).text={
   '                                                     Knot selection',
   'When using an SP... command to construct a spline, one usually has to specify a particular spline space. This is done by specifying a   k n o t  s e q u e n ce  and an  o r d e r , and this may be a bit of a problem.',
   'When doing spline interpolation, to data  X , Y by splines of order K , then OPTKNT will supply a good knot sequence, as in the following example:',
   '>> k = 5;   % i.e., we are working with quartic splines',
   '>> x = 2*pi*sort([0 1 rand(1,10)]); y = cos(x);',
   '>> sp = spapi( optknt(x,k), x, y );',
   '>> fnplt(sp,2,''g''), hold on; plot(x,y,''o''); hold off'};

  %========== Slide 10 ==========

  slide(10).code={
   'x = linspace(0,10,101); y = exp(x);',
   'sp0 = spap2( augknt( 0:2:10, 4), 4 ,x, y ); plot(x, y-fnval(sp0,x), ''r'' )',
   'sp1 = spap2( newknt(sp0), 4,x,y);',
   'hold on, plot(x,y-fnval(sp1,x),''k'',''linew'',2), hold off',
   '' };
  slide(10).text={
   '                                                     Knot selection (cont.)',
   'When doing least-squares approximation, one can use the current approximation',
   'to determine a possibly better knot selection with the aid of NEWKNT.',
   '',
   '>> x = linspace(0,10,101); y = exp(x);',
   '>> sp0 = spap2( augknt(0:2:10,4), 4, x, y );',
   'whose error is plotted above in red, and which isn''t all that good,',
   'compared to the following approximation with the  s a m e  number of knots,',
   ' but better distributed (whose error is plotted in black):',
   '>>sp1 = spap2( newknt(sp0), 4, x, y );',
   ''};

  %========== Slide 11 ==========

  slide(11).code={
   'x =.0001+[-4:.2:4]; y = -3:.2:3;',
   '[yy,xx] = meshgrid(y,x); r = pi*sqrt(xx.^2+yy.^2); z = sin(r)./r;',
   'bcs = csapi( {x,y}, z ); fnplt( bcs ), axis([-5 5 -5 5 -.5 1])' };
  slide(11).text={
   '                                                     Gridded data',
   '',
   'All the interpolation and approximation commands in the toolbox can also',
   'handle gridded data (in any number of variables).',
   '',
   'For example, here is a bicubic spline interpolant to the mexican hat function:',
   '>> x =.0001+[-4:.2:4]; y = -3:.2:3;',
   '>> [yy,xx] = meshgrid(y,x); r = pi*sqrt(xx.^2+yy.^2); z = sin(r)./r;',
   '>> bcs = csapi( {x,y}, z ); fnplt( bcs )'};

  %========== Slide 12 ==========

  slide(12).code={
   ' knotsx = augknt(linspace(x(1), x(end), 21), 4);',
   ' knotsy = augknt(linspace(y(1), y(end), 15), 4);',
   ' bsp2 =  spap2({knotsx,knotsy},[4 4], {x,y},z+.02*(rand(size(z))-.5));',
   ' fnplt(bsp2) , axis([-5 5 -5 5 -.5 1])' };
  slide(12).text={
   '... and here is the  l e a s t - s q u a r e s  approximation to noisy values',
   'of that function on the same grid:',
   '',
   '>> knotsx = augknt(linspace(x(1), x(end), 21), 4);',
   '>> knotsy = augknt(linspace(y(1), y(end), 15), 4);',
   '>> bsp2 = spap2( {knotsx,knotsy}, [4 4], {x,y} ,z+.02*(rand(size(z))-.5) );',
   ' fnplt(bsp2)',
   ''};

  %========== Slide 13 ==========

  slide(13).code={
   ' t = 0:8; xy = [0 0;1 1; 1.7 0;1 -1;0 0; -1 1; -1.7 0; -1 -1; 0 0].'';',
   ' infty = csape( t , xy, ''periodic'');  fnplt( infty , 2 )',
   ' axis([-2 2 -1.1 1.1]), hold on, plot(xy(1,:),xy(2,:),''o''),hold off' };
  slide(13).text={
   '                                                       Curves',
   '',
   'Gridded data can be handled easily because the toolbox can deal with',
   '  v e c t o r  -  v a l u e d   splines.  This also makes it easy to work with curves.',
   '',
   'Here, for example, is an approximation to infinity (obtained by putting',
   'a cubic spline curve through the marked points):',
   '>> t = 0:8; xy = [0 0;1 1; 1.7 0;1 -1;0 0; -1 1; -1.7 0; -1 -1; 0 0].'';',
   '>> infty = csape( t , xy, ''periodic'');  fnplt( infty , 2 ), axis([-2 2 -1.1 1.1])'};

  %========== Slide 14 ==========

  slide(14).code={
   'roller = csape( t , [ xy ;0 1/2 1 1/2 0 1/2 1 1/2 0], ''periodic'');',
   'fnplt( roller , 2, [0 4],''b'' ) , hold on, fnplt( roller, 2, [4 8], ''r''), plot3(0,0,0,''o'')  , hold off' };
  slide(14).text={
   '... and here is the same curve, but with motion in a third dimension:',
   '',
   '>>roller = csape( t , [ xy ;0 1/2 1 1/2 0 1/2 1 1/2 0], ''periodic'');',
   '>>fnplt( roller , 2, [0 4],''b'' ) , hold on, fnplt( roller, 2, [4 8], ''r''), plot3(0,0,0,''o'')',
   '>> hold off',
   '',
   'I have plotted the two halves of the curve in different colors and have marked the origin, as an aid to visualizing this two-winged space curve.'};

  %========== Slide 15 ==========

  slide(15).code={
   'x = 0:4; y=-2:2;  R = 4; r = 2; clear v',
   'v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].''*[1 1 1 1 1];',
   'v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[0 1 0 -1 0];',
   'v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[1 0 -1 0 1];',
   'dough0 = csape({x,y},v,''periodic'');',
   'fnplt(dough0)',
   'axis equal, axis off',
   '',
   '' };
  slide(15).text={
   '                                                        Surfaces',
   '',
   'Finally, bivariate tensor-product splines with values in R^3 give surfaces.',
   '',
   'For example, here is a good approximation to a doughnought:',
   '>> x = 0:4; y=-2:2;  R = 4; r = 2;',
   '>> v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].''*[1 1 1 1 1];',
   '>> v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[0 1 0 -1 0];',
   '>> v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[1 0 -1 0 1];',
   '>> dough0 = csape({x,y},v,''periodic''); fnplt(dough0)',
   ''};

  %========== Slide 16 ==========
  slide(16).code={
   'x = 0:4; y=-2:2;  R = 4; r = 2; clear v',
   'v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].''*[1 1 1 1 1];',
   'v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[0 1 0 -1 0];',
   'v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].''*[1 0 -1 0 1];',
   'dough0 = csape({x,y},v,''periodic'');',
   'fnplt(dough0)',
   'axis equal, axis off',
   'nx = 43; xy = [ones(1,nx);linspace(2,-2,nx)]; points = fnval(dough0,xy)'';',
   'ders = fnval(fndir(dough0,eye(2)),xy);',
   'normals = cross(ders(4:6,:),ders(1:3,:));',
   'normals = (normals./repmat(sqrt(sum(normals.*normals)),3,1))'';',
   'pn = [points;points+normals];'
   'hold on, for j=1:nx',
   '   plot3(pn([j,j+nx],1),pn([j,j+nx],2),pn([j,j+nx],3)),end,hold off',
   '' };
  slide(16).text={
   '                                                     Surfaces (cont.)',
   '',
   'And here is a crown of normals to that surface.',
'>> nx = 43; xy = [ones(1,nx);linspace(2,-2,nx)]; points = fnval(dough0,xy)'';',
'>> ders = fnval(fndir(dough0,eye(2)),xy);',
'>> normals = cross(ders(4:6,:),ders(1:3,:));',
'>> normals = (normals./repmat(sqrt(sum(normals.*normals)),3,1))'';',
'>> pn = [points;points+normals];'
'>> hold on, for j=1:nx',
'>>    plot3(pn([j,j+nx],1),pn([j,j+nx],2),pn([j,j+nx],3)) ,end, hold off',
   ''};

  %========== Slide 17 ==========

  slide(17).code={
   'fnplt(fncmb(dough0, [1 0 0; 0 1 0]))',
   'axis equal, axis off',
   '',
   '' };
  slide(17).text={
   '                                                     Surfaces (cont.)',
   '',
   'And here is its projection onto the (x,y)-plane.',
   '',
   '>> fnplt( fncmb( dough0, [1 0 0; 0 1 0]) )',
   ''};
end
