function c = CylinderIntercept(u,v,radius,BoxLength1, BoxLength2)
%Compute the coordinates of the intercept of a box surface and the line conecting two points
%intercept point is c, k is vector connecting points, n is distance along k that c is located

%Check if sphere is between vectors u and v
swaped = 0;
if (inCylinder(v,radius,BoxLength1, BoxLength2) & ~inCylinder(u,radius,BoxLength1, BoxLength2))
   %moving into the cylinder as oposed to leaving it.
   temp = u;
   u = v;
   v = temp;
   swapped = 1;
end
if or(swaped,(inCylinder(u,radius,BoxLength1, BoxLength2) & ~inCylinder(v,radius,BoxLength1, BoxLength2)))
   %OK to compute intercept and u is the inner vector
      %cuved wall computation
      k = v-u;
      a=k(1)^2+k(2)^2;
      b=2*(k(1)*u(1)+k(2)*u(2));
      c=u(1)^2+u(2)^2-radius^2;
      n=(-b-sqrt(b^2-4*a*c))/(2*a);
      if n<0
         n=(-b+sqrt(b^2-4*a*c))/(2*a);
      end
      DistanceToCurvedWall = norm(n*k);
      
   if v(3)<BoxLength1 & u(3)>BoxLength1% lower endcap could be first
      nend = abs((BoxLength1-u(3))/(u(3)-v(3)));
      DistanceToWall = norm(nend*(v-u));
      
      if DistanceToWall<DistanceToCurvedWall
         c = u + nend*(v-u);
         c(3) = BoxLength1; %above calulation not perfict due to roundings?
      else %curve wall first
         c=u+n*k;
      end
            
   elseif v(3)>BoxLength2 & u(3)<BoxLength2 % upper endcap
      nend = abs((BoxLength2-u(3))/(u(3)-v(3)));
      DistanceToWall = norm(nend*(v-u));
      
      if DistanceToWall<DistanceToCurvedWall
         c = u + nend*(v-u);
         c(3) = BoxLength2; %above calulation not perfict due to roundings?
      else %curve wall first
         c=u+n*k;
      end
     
   else %hit a curved wall and an edge is not involved
      c=u+n*k;
   end
         
else
   c = [0;0;0];
end
