% IAL n-point quadrature approximation to double integrals 
%     over the unit square with algebraic and algebraic/
%     logarithmic singularities along the diagonal
%     diagonal
%
function [Ialg,Ilog]=Ial(n)
global abj abl 
A=zeros(n); B=zeros(n); C=zeros(n);
xwj=gauss(n,abj); xwl=gauss(n,abl);
xj=xwj(:,1); wj=xwj(:,2);
xl=xwl(:,1); wl=xwl(:,2);
for i=1:n
  for j=1:n
    A(i,j)=wj(i)*xj(i)*wj(j)*(f((1-xj(j))*xj(i),xj(i)) ...
     +f(1-(1-xj(j))*xj(i),1-xj(i)));
    B(i,j)=wl(i)*xl(i)*wj(j)*(f((1-xj(j))*xl(i),xl(i)) ...
     +f(1-(1-xj(j))*xl(i),1-xl(i)));
    C(i,j)=wj(i)*xj(i)*wl(j)*(f((1-xl(j))*xj(i),xj(i)) ...
     +f(1-(1-xl(j))*xj(i),1-xj(i)));
  end
end
Ialg=sum(sum(A));
Ilog=-sum(sum(B+C));
%
% The routine for the function f, and other supporting routines,
% if any, must be appended as subfunctions to this routine.
%
