%
%  gibbs.m  -- djm, 20 june 05
%

xp = 0:0.02:8;
yp = zeros(size(xp));

for jj = 2:size(xp,2)
    yp(jj) = quad(@sinc,0,xp(jj));
end

err = 2*(0.5-yp);

figure(1);  clf;  hold on
plot(xp,err,'b')
ylabel('\bf 2(f(x)-F_N(x))')
xlabel('\bf (N + 1/2) x/\pi')

figure(2);  clf;  hold on
th = 0.02:0.02:2;

g = 1./sin(th/2) - (2./th);
plot([0 th],[0 g],'b')
ylabel('\bf g(\theta)')
ylabel('\bf g(\theta)')
xlabel('\bf \theta')