function [ n, m ] = circularscanline( fractureShp, circleCentre, radius)
% 
% [ n, m ] = circularscanline( fractureShp, circleCentre, radius)
% Applies the Circular Scanline fracture sampling method of Mauldon et al, 2001 to a
% vectorized fracture map. 
% 
% Input:
% The fracture map is specified by fractureShp,
% which must be shapefile structure as produced by the 'shaperead'
% function. Circlecentre is a 2x1 vector specifying the centre coordinates
% of the counting circle, radius is a scalar specifying the radius of the
% counting circle. 
% 
% Output:
% n : the number of fracture intersections with the circular scanline.
% m : the number of fracture end points inside the circle.
% 
% M. Mauldon, W.M. Dunne, and M.B. Rohrbaugh, “Circular scanlines and circular windows: new tools for characterizing the geometry of fracture traces,” Journal of Structural Geology,  vol. 23, Feb. 2001, pp. 247-258.


% circle equation: (x-circleX)^2 + (y-circleY)^2 = radius^2
%                   x = sqrt( radius^2 -( y - circleY)^2 + circleX; pwèèèèèp
% 
% Line equation  y = mx +b
%                y = dy/dx 
xend = [];
yend=[];
xtwoint = [];
ytwoint = [];
xoneint = [];
yoneint = [];


n = 0; % Scanline intersections
m = 0; % TRace endpoints in Circle
for fractureIdx = 1 : length( fractureShp),
   l = length( fractureShp( fractureIdx ).X)-1;
   
   if l>1,
     x = [fractureShp( fractureIdx ).X];
     y = [fractureShp( fractureIdx ).Y];
     intersectingscanline = false;
     
     for idx = 1 : l-1,
        a1 = [x(idx), y(idx)];
        a2 = [x(idx + 1), y( idx + 1 )];
        
        intersectingscanline = or( intersectingscanline, linecircleintersect2( a1, a2, circleCentre, radius, false) );               
     end
     
     firstPointInCircle = sqrt( sum( ([x(1),y(1)]-circleCentre).^2) ) < radius;
     lastPointInCircle = sqrt( sum( ([x(l),y(l)]-circleCentre).^2) ) < radius;
     
     %    if firstPointInCircle && lastPointInCircle, n = n+0; end;
     if xor( firstPointInCircle, lastPointInCircle) && intersectingscanline, 
         n = n+1; 
         xoneint = [xoneint, x(1:end)];
         yoneint = [yoneint, y(1:end)];
%          BUG in here, if the line intersects the circle, but the line
%          segment doed not intersect with the normal, it is not counted.
     end; % ONE endpoint in circle, and intersecting, so one intersection
     if ~( firstPointInCircle || lastPointInCircle) && intersectingscanline, 
         n = n+2; 
         xtwoint = [xtwoint, x(1:end)];
         ytwoint = [ytwoint, y(1:end)];
     end; % no endpoints in circle and intersecting, so two intersections 
     if firstPointInCircle,   
         m = m+1; 
         xend = [xend, x(1)]; 
         yend = [yend, y(1)]; 
     end
     if lastPointInCircle,    
         m = m+1; 
         xend = [xend, x(l)]; 
         yend = [yend, y(l)]; 
     end
     
   end

       
    
end