% Matlab program, Jorg Peters
% apply generic subdivison algorithm
% in      s=1 dimensions
%         n   times
% with    d+1 averaging steps

clear;
clf;

    d = 3;
    n = 5;

    prt = 1; % 0 = do not print figures
    colr = ['r', 'g', 'c', 'm', 'k'];  % color scheme for subsequent iterates

    values = [2 -1 1 2 -1 2 0];
    seq = values;   % seq = sequence of numbers
    N = length(seq); N0 = N;
    left = 0;  % leftmost index

    plot([-1,N], [min(seq)-1, max(seq)+1],'w.'); hold on; 
    plot([0:N-1], seq, '+'); axis off; plot([0:N-1], seq, ':');
	if prt==1, print -dpsc uni0.eps; end;
    pause;

    for sub=1:n,
	spacezero = [seq; zeros(1,N)];
	seq = (2*spacezero(:))';
	for avg=0:d,
	    seq = ( seq(1:2*N-1-avg) + seq(2:2*N-avg) )/2;
	end;
	N = length(seq);
	left = 2*left + (d+1)/2;
	plot((left+[0:N-1])/2^sub, seq, [colr(sub) '+']);
	plot((left+[0:N-1])/2^sub, seq, [colr(sub) ':']);
	sub
	    if prt==1,
		eval(['print -depsc uni', sprintf('%d',sub), '.eps']);
	    end;
	pause;
    end;
    clf;
    plot([-1,N0], [min(values)-1, max(values)+1],'w.'); hold on;  axis off;
    plot((left+[0:N-1])/2^n, seq,colr(n));
    plot([0:N0-1], values,'+'); plot([0:N0-1], values,':');
    tit = sprintf('spline of degree %d', d);
    title (tit);
    if prt==1, print -dpsc unifin.eps; end;