[NoNewPts,IsetClipPoly,NoInterpPts, ...

function [NoNewPts,IsetClipPoly,NoInterpPts, ...
  Xclip,Yclip]=clip(Xclip,Yclip,NoClipPts, ...
  Line,Epsilon)
%
% [NoNewPts,IsetClipPoly,NoInterpPts,Xclip, ...
%   Yclip]=clip(Xclip,Yclip,NoClipPts, ...
%   Line,Epsilon)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function clips a reference polygon
% using a line defined in vector Line.
%
% Xclip,Yclip  - coordinates of reference 
%                polygon on input
% NoClipPts    - length of Xclip/Yclip
% Line         - vector containing coefficients
%                of line equation in form:
%                Line(1)+Line(2)*X+Line(3)*Y=0
%                The normal vector to this 
%                line should be directed away
%                from the region to be kept.
% Epsilon      - for determining zero values
%
% NoNewPts     - length of Xclip/Yclip on
%                output
% IsetClipPoly - vector containing the 
%                indices for the clipped
%                polygon in correct order
% NoInterpPts  - number of new points
%                defined by interpolation
% Xclip,Yclip  - coordinates of clipped
%                polygon on output
%
% User m functions called:  inout, intrsect
%----------------------------------------------

%...Initialize
iPtConsider=1; iStartPoint=1; NoInterpPts=0; 
NoNewPts=0; IdxNewPt=NoClipPts+1;
[IsFptIn]=inout(Xclip(iPtConsider), ...
                Yclip(iPtConsider),Line);

%...Loop on each side
for j=1:NoClipPts
  iStartPoint=j; iPtConsider=j+1;
  if j == NoClipPts, iPtConsider=1; end;

  [IsNptIn]=inout(Xclip(iPtConsider), ...
                  Yclip(iPtConsider),Line);

  if IsFptIn == 0 & IsNptIn == 1
    %...Current point is outside and next point 
    %...is inside. Therefore, generate 
    %...intermediate points.
    [IsectErr,Xclip,Yclip]=intrsect ...
      (Xclip,Yclip,Line,iStartPoint, ...
      iPtConsider,IdxNewPt,Epsilon);
    if IsectErr == 1
      error('\n\nError #1 in clip');
    end
    NoNewPts=NoNewPts+1;
    IsetClipPoly(NoNewPts)=IdxNewPt;
    IdxNewPt=IdxNewPt+1;
    NoInterpPts=NoInterpPts+1;

  elseif IsFptIn == 1
    %...Current point inside, so save it
    NoNewPts=NoNewPts+1;
    IsetClipPoly(NoNewPts)=j;
    if IsNptIn == 0
      %...Next point is outside 
      %...Generate intermediate point
      [IsectErr,Xclip,Yclip]=intrsect ...
        (Xclip,Yclip,Line,iStartPoint, ...
        iPtConsider,IdxNewPt,Epsilon);
      if IsectErr == 1
        error('\n\nError #2 in clip');
      end
      NoNewPts=NoNewPts+1;
      IsetClipPoly(NoNewPts)=IdxNewPt;
      IdxNewPt=IdxNewPt+1;
      NoInterpPts=NoInterpPts+1;
    end
  end
  IsFptIn=IsNptIn;
end