% **********************************************************************
%    Subroutine OkQuadrature.m by www.boundary-element-method.com      *
% **********************************************************************

% OkQuadrature.m Version 1, Sept 2015
% Stephen Kirkup, School of Engineering
% University of Central Lancashire
% http://www.researchgate.net/profile/Stephen_Kirkup

% A test on the quadrature rule input to *2lc (and *alc) routines in 
% the 2D and axisymmetric boundary element method.
% Source of this code:
%  http://www.boundary-element-method.com/mfiles/OkQuadrature.m
% Source of the user-guide:
%  http://www.boundary-element-method.com/mfiles/OkQuadrature_m.pdf
% Source of the test code:
%  http://www.boundary-element-method.com/mfiles/OkQuadrature_tests.m

function [passed]=OkQuadrature(n, w, x, xStar, tolQuad)
    passed=true;
    
% Test 1: Verify tolQuad
%  Verify that tolQuad is positive
%   Tests 1(a)
    if (tolQuad <= 0) 
        passed = false;
        'Error(*lc)- Test 1(a) : tolerance (tolQuad) must be positive'
    end
%  Verify that tolQuad is not too large
%   Tests 1(b)
    if (tolQuad > 0.001) 
        passed = false;
        'Error(*lc) - Test 1(b) : tolerance (tolQuad) is too large'
    end
    
% Test 2: Verify n is positive
    if (n <= 0) 
        passed = false;
        'Error(*lc)- Test 2 : number of quadrature points (n) must be positive'
    end 

% Test 3: Verify that all of the abscissae are in [0,1]
    LocalError = false;
    for k = 1:n
        if (x(k) < 0 | x(k) > 1) 
            LocalError = True;
        end
    end
   if (LocalError) 
       'Error(*lc) - Test 3 : quadrature abscissae must be on [0,1]'
       passed = false;
   end
   
% Test 4: Verify that all of the weights sum to one
   Sum = 0.0;
   for k = 1:n
       Sum = Sum + w(k);
   end
   if (abs(Sum - 1) > n * tolQuad) 
      'Error(*lc) - Test 4 : quadrature weights must sum to one'
      passed = false;
   end

% Test 5: Verify that the integral of f(x)=x is 0.5
   Sum = 0.0;
   for k = 1:n
      Sum = Sum + w(k) * x(k);
   end
   if (abs(Sum - 0.5) > n * tolQuad) 
      'Error(*lc) - Test 5 : integration of f(x)=x failed'
      passed = false;
   end
  
% The following tests are particularly for the case in which there is
%  a possible discontinuity at xStar

   if (xStar > 0 & xStar < 1) 

% Test D1: Verify x-discontinuity for f(x)=1 (x<xStar),
%  f(x)=-1 (x>xStar)
       Sum = 0.0;
       for k = 1:n
           if (x(k) < xStar) 
               Sum = Sum + w(k);
           else
               Sum = Sum - w(k);
           end 
       end
       if (abs(Sum - 1 + 2 * xStar) > n * tolQuad) 
           'Error(*lc) - Test D1 : split quadrature rule must integrate function with discontinuity at x_discontinutiy accurately'
           passed = false;
       end
  
% Test D2: Verify x-discontinuity for f(x)=1 (x<xStar),
%  f(x)=1-x (x>xStar)
       Sum = 0.0;
       for k = 1:n
               if (x(k) < xStar) 
               Sum = Sum + w(k);
           else
               Sum = Sum - w(k) * 0.5;
           end
       end
       if (abs(Sum + 0.5 - 1.5 * xStar) > n * tolQuad) 
           'Error(*lc) - Test D2: split quadrature rule must integrate function with discontinuity at x_discontinutiy accurately'
           passed = false;
       end 

% Test D3: Verify x-discontinuity for f(x)=x (x<xStar),
%  f(x)=1 (x>xStar)
       Sum = 0.0;
       for k = 1:n
           if (x(k) < xStar)
               Sum = Sum + w(k) * x(k);
           else
               Sum = Sum + w(k);
           end
       end
       if (abs(Sum - 1 + xStar - 0.5 * xStar ^ 2) > n * tolQuad) 
           'Error(*lc) - Test D3: split quadrature rule must integrate function with discontinuity at x_discontinutiy accurately'
           passed = false
       end 
   end 


    