% INTSIN01 Integration relative to the weight function 
%    w(t)=sin(1/t+t)
%
f1='%5.0f %21.13e\n';
fprintf('\n')
disp('    n     n-point Gauss')
load -ascii absin0; abs=absin0;
load -ascii abcos0; abc=abcos0;
abl=r_jacobi01(40);
for n=1:7
  xwl=gauss(n,abl); xws=gauss(n,abs); xwc=gauss(n,abc);
  intl=sum(xwl(:,2).*exp(-xwl(:,1)).*(cos(xwl(:,1))...
    +sin(xwl(:,1))));
  ints=sum(xws(:,2).*exp(-xws(:,1)).*cos(xws(:,1)));
  intc=sum(xwc(:,2).*exp(-xwc(:,1)).*sin(xwc(:,1)));
  int=ints+intc-intl;
  fprintf(f1,n,int)
end
fprintf('\n')
abL=r_jacobi01(1000);
for N=200:200:1000
  xwL=gauss(N,abL);
  intL=sum(xwL(:,2).*exp(-xwL(:,1)).*sin(1./xwL(:,1)...
    +xwL(:,1)));
  fprintf(f2,N,intL)
end