function [pass,xint,yint,jint,kint,kfrac] = seg(xtrack,ytrack,xa,ya,xb,yb)
%
%     This function determines whether a set of n tracks (given by 
%     coordinates xtrack,ytrack) intersects a series of m-1 line segments
%     whose end points are given by xa(m), ya(m), xb(m), yb(m). Pass=0  
%     means no intersection; pass=1 means intersection. 
%
%     The quantities xint and yint constitute the intersection points,
%     along the line segment numbered jint and the kint'th point along the
%     track. kfrac is the fraction of the distance between track points i
%     and i+1 that the intersection occurs. 
%
%     Note that pass is a vector of length n, where n is the number of
%     tracks, and the other quantities are matrices of dimension
%     (n,maxtimes), where maxtimes is the specified upper limit on the
%     number of times a particular track may intersect the set of line
%     segments. 
%
%     Copyright 2012 WindRiskTech L.L.C.
%     Last modified September, 2012.  
%-------------------------------------------------------------------------
maxtimes=10;    % Maximum number of times a track may intersect the lines segments
%-------------------------------------------------------------------------
[n,~]=size(xtrack);
m=max(size(xa));
for i=1:m
    if xa(i) < 0
        xa(i)=xa(i)+360;
    end
    if xb(i) < 0
        xb(i)=xb(i)+360;
    end
end    
%
% default values
%
pass=zeros(1,n);
jint=zeros(n,maxtimes);
kint=zeros(n,maxtimes);
xint=zeros(n,maxtimes);
yint=zeros(n,maxtimes);
kfrac=zeros(n,maxtimes);
%
[~, jmax]=min(abs(ytrack),[],2);
jmax=jmax-1;
jmax=max(jmax,2);
%
for k=1:n
    %
    nj=0;
    for i=1:jmax(k)-1
        x1=xtrack(k,i);
        x2=xtrack(k,i+1);
        y1=ytrack(k,i);
        y2=ytrack(k,i+1);
        for j=1:m
              x3=xa(j);
              y3=ya(j);
              x4=xb(j);
              y4=yb(j);
    %
              denom=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1);
              if abs(denom) >= 1e-8
                   ua=((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/denom;
                   ub=((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/denom;
                   if (ua >= 0 && ua <= 1) && (ub >= 0 && ub <= 1)
                        nj=nj+1;   
                        pass(k)=1;
                        xint(k,nj)=x1+ua*(x2-x1);
                        yint(k,nj)=y1+ua*(y2-y1);
                        jint(k,nj)=j;
                        kint(k,nj)=i;
                        kfrac(k,nj)=ua;
                   end 
              end 
              if nj == maxtimes
                  break
              end    
        end 
    end 
%
end
%
end
        