% **********************************************************************
%    Subroutine OkGeometry2lc.m by www.boundary-element-method.com     *
% **********************************************************************
% OkGeometry2lc.m Version 1, Sept 2015

% Matlab/Freemat/Octave/Scilab open source codes
% Stephen Kirkup, School of Engineering
% University of Central Lancashire
% http://www.researchgate.net/profile/Stephen_Kirkup

%  Checks that the given geometry is suitable for the boundary element
%  method.
% Source of thiS code:
%  http://www.boundary-element-method.com/mfiles/OkGeometry2lc.m
% Source of the user-guide:
%  http://www.boundary-element-method.com/mfiles/OkGeometry2lc_m.pdf
% Source of the test code:
%  http://www.boundary-element-method.com/mfiles/OkGeometry2lc_tests.m


function [passed]=OkGeometry2lc(p, vecp, qa, qb, lponel, tolGeom)
qbma=qb-qa;
pmqa=p-qa;
pmqb=p-qb;
qlen=vecsize(qbma);
pqalen=vecsize(pmqa);
pqblen=vecsize(pmqb);
passed = true;

% Verify tolGeom
%  Verify that tolGeom is positive
if (tolGeom <= 0) 
  passed = false;
  'Error(*2lc): tolGeom must be positive'
end
% Verify that tolGeom is not too large
if (tolGeom > (qlen + pqalen + pqblen) * 0.001) 
   passed = false;
  'Error(*2lc) - tolGeom is too large'
end

% Verify that qa and qb are not the same point
if (qlen < tolGeom * (vecsize(qb) + vecsize(qa) + 2))
   passed = false;
  'Error(*2lc) - The points qa and qb must be distinct'
end

% Verify vecp
%  Verify that vecp has unit magnitude
if (abs(vecsize(vecp) - 1.0) > tolGeom) 
   passed = false;
  'Error(*2lc) - vecp must be a unit vector'
end
%  Verify that vecp is equal to the unit normal to the panel when
%   the point p lies on the panel (lponel=True)
if (lponel)
  normq=norm2(qa,qb); 
  diff=normq-vecp;  
  if (vecsize(diff) > 2 * tolGeom)
     passed = false;
    'Error(*2lc) - If lponel is true then vecp must be the unit leftward normal to panel qa-qb'
  end
end

% Verify lponel
%  Set qStar; the point on the panel qa-qb that is nearest to p
if (passed)
    xStar = (qbma*pmqa')/ (qlen ^ 2);
  if (xStar <= 0) 
    qStar = qa;
  elseif (xStar >= 1) 
    qStar = qb;
  else
    qStar = qa + xStar * qbma;
  end
%  If lponel=true then qStar must be the same point as p.
  if (lponel) 
    if (dist(p, qStar) > tolGeom * (vecsize(qStar) + vecsize(p) + 2)) 
       passed = false;
      'Error(*2lc) - If lponel=true then the point p must lie on qa-qb'
    end
%  If lponel=false then qStar must not be the same point as p.
  else
    if (dist(p, qStar) < tolGeom * (vecsize(p) + vecsize(qa) + 2))
       passed = false;
      'Error(*2lc) - If lponel=false then the point p must not lie on qa-qb'
    end
  end
end

% Verify that p is not coincident with qa
if (dist(p, qa) < tolGeom * (vecsize(p) + vecsize(qa) + 2))
   passed = false
  'Error(*2lc) - The points p and qa must be distinct'
end

% Verify that p is not coincident with qb
if (dist(p, qb) < tolGeom * (vecsize(p) + vecsize(qb) + 2))
   passed = false;
  'Error(*2lc) - The points p and qb must be distinct'
end

% Warning if lponel=false and the point p is too close to qa-qb
if (~lponel & passed) 
  if (dist(p, qStar) < qlen / 10 & passed)
  'Warning(*2lc) - If lponel=false then if the observaton point p is close to' 
  ' qa-qb the expected error in solution likely to be magnified'
  end
end

% Warning if lponel=true and the point p is too close to qa or qb
if (lponel & passed)
  if (dist(p, qa) < qlen / 10 | dist(p, qb) < qlen / 10) 
  'Warning(*2lc) - If lponel=true then if the observaton point p is close to qa or qb'
  ' the expected error in solution likely to be magnified'
  end
end
