line_cylinder_intersect Subroutine

subroutine line_cylinder_intersect(l0, l, p0, p, t)
    !+ Calculates the intersection of a line and a cylinder
    real(Float64), dimension(3), intent(in)  :: l0
        !+ Point on line
    real(Float64), dimension(3), intent(in)  :: l
        !+ Ray of line
    real(Float64), dimension(3), intent(in)  :: p0
        !+ Point on cylinder
    real(Float64), dimension(3), intent(out) :: p
        !+ Line-cylinder intersect point
    real(Float64), intent(out)               :: t
        !+ "time" to intersect

    real(Float64), dimension(2) :: times
    logical, dimension(2) :: mask
    real(Float64) :: r, vx, vy, x0, y0
    real(Float64) :: radicand, npos

    r = sqrt(p0(1) * p0(1) + p0(2) * p0(2))
    x0 = l0(1) ; y0 = l0(2)
    vx = l(1)  ; vy = l(2)

    if((vx.eq.0.d0).and.(vy.eq.0.d0)) then
        t = 0.d0        ! Parallel to a plane tangent to the cylinder
    else
        radicand = r**2 * (vx**2 + vy**2) - (vy * x0 - vx * y0)**2
        if(radicand.lt.0) then
            t = 0.d0    ! Parallel to a plane tangent to the cylinder
        else
            times(1) = (- vx * x0 - vy * y0 - sqrt(radicand)) / (vx**2 + vy**2)
            times(2) = (- vx * x0 - vy * y0 + sqrt(radicand)) / (vx**2 + vy**2)
            mask = times.gt.0
            npos = count(mask)
            if(npos.gt.0) then
                t = minval(times, mask=times.gt.0)
            else
                t = maxval(times, mask=times.le.0)
            endif
        endif
    endif

    p = l0 + l * t

end subroutine line_cylinder_intersect