% ISTARAL n-point quadrature approximation to double 
%     integrals over the unit square with algebraic 
%     singularities along the sides and logarithmic 
%     singularities along the diagonal
%
function [Istara,Istarl]=Istaral(n)
global alpha beta abj abj1 abj2 abj3 abj4 abl1 abl2 abl3 abl4 
A=zeros(n); B=zeros(n); C=zeros(n); 
D=zeros(n); E=zeros(n);
xwj=gauss(n,abj);
xwj1=gauss(n,abj1); xwj2=gauss(n,abj2); 
xwj3=gauss(n,abj3); xwj4=gauss(n,abj4);
xwl1=gauss(n,abl1); xwl2=gauss(n,abl2);
xwl3=gauss(n,abl3); xwl4=gauss(n,abl4);
xj=xwj(:,1); wj=xwj(:,2);
xj1=xwj1(:,1); xj2=xwj2(:,1); xj3=xwj3(:,1); xj4=xwj4(:,1);
wj1=xwj1(:,2); wj2=xwj2(:,2); wj3=xwj3(:,2); wj4=xwj4(:,2);
xl1=xwl1(:,1); xl2=xwl2(:,1); xl3=xwl3(:,1); xl4=xwl4(:,1);
wl1=xwl1(:,2); wl2=xwl2(:,2); wl3=xwl3(:,2); wl4=xwl4(:,2);
for i=1:n
  for j=1:n
    A(i,j)=wj(i)*wj(j)*f(xj(j),xj(i));
    B(i,j)=wl3(i)*wj2(j)*f((1-xj2(j))*xl3(i),xl3(i)) ... 
     *(1-xl3(i)+xj2(j)*xl3(i))^alpha;
    C(i,j)=wl4(i)*wj1(j)*f(1-xl4(i)+xj1(j)*xl4(i),1-xl4(i)) ...
     *(1-xl4(i)+xj1(j)*xl4(i))^beta;
    D(i,j)=wj3(i)*wl2(j)*f((1-xl2(j))*xj3(i),xj3(i)) ...
    *(1-xj3(i)+xl2(j)*xj3(i))^alpha;
    E(i,j)=wj4(i)*wl1(j)*f(1-xj4(i)+xl1(j)*xj4(i),1-xj4(i)) ...
    *(1-xj4(i)+xl1(j)*xj4(i))^beta;
  end
end
Istara=sum(sum(A));
Istarl=-sum(sum(B+C+D+E));

% 
% The routine for the function f, and other supporting routines,
% if any, must be appended as subfunctions to this routine.
%
