euclidean_geo_m.f90 Source File

module euclidean_geo_m
    !! Collection of routines related with Euclidean geometry
    use precision_m, only: FP, FP_EPS, FP_NAN
    implicit none
    private

    public :: intersection_lines

contains

    subroutine intersection_lines(pa, qa, pb, qb, xi, yi, info, ta, tb)
        !! Computes intersection point of two lines 
        !! going through points [pa, qa] and [pb, qb]
        real(FP), intent(in), dimension(2) :: pa
        !! x/y-coordinate of first point of first line
        real(FP), intent(in), dimension(2) :: qa
        !! x/y-coordinate of second point of first line
        real(FP), intent(in), dimension(2) :: pb
        !! x/y-coordinate of first point of second line
        real(FP), intent(in), dimension(2) :: qb
        !! x/y-coordinate of second point of second line
        real(FP), intent(out) :: xi
        !! If info = 0: x-coordinate of intersection point
        !! If info /= 0: NaN
        real(FP), intent(out) :: yi
        !! If info = 0: y-coordinate of intersection point
        !! If info /= 0: NaN
        integer, intent(out) :: info
        !! Returns info on intersection point
        !! info = 0: Unique intersection point found
        !! info = 1: Both lines are identical
        !! info = -1: Both lines are parallel and not identical
        !! info = -2: Some line is degenerate, i.e. pa = qa or pb = qb
        real(FP), intent(out), optional :: ta
        !! If info = 0: Parameter where intersection occurs w.r.t. first line,
        !! i.e. pa + ta * (pa - qa)
        !! If info /= 0: NaN
        real(FP), intent(out), optional :: tb
        !! If info = 0: Parameter where intersection occurs w.r.t. second line,
        !! i.e. pb + tb * (pb - qb)
        !! If info /= 0: NaN

        real(FP) :: mx, my, nx, ny, dx, dy, det, ovl, ra, rb
        logical :: gax, gay, gbx, gby

        ! Check if any line is degenerate (both points of any line are equal)
        gax = ( abs(qa(1) - pa(1)) < FP_EPS )
        gay = ( abs(qa(2) - pa(2)) < FP_EPS )
        gbx = ( abs(qb(1) - pb(1)) < FP_EPS )
        gby = ( abs(qb(2) - pb(2)) < FP_EPS )

        if ( (gax .and. gay) .or. (gbx .and. gby) ) then
            info = -2
            xi = FP_NAN
            yi = FP_NAN
            if (present(ta)) then
                ta = FP_NAN
            endif
            if (present(tb)) then
                tb = FP_NAN
            endif
            return
        endif

        ! Determine slope of lines
        mx = qa(1) - pa(1)
        my = qa(2) - pa(2)

        nx = qb(1) - pb(1)
        ny = qb(2) - pb(2)

        dx = pb(1) - pa(1)
        dy = pb(2) - pa(2)

        det = mx * ny - my * nx
        
        ! Check if lines are parallel or identical
        if (abs(det) < FP_EPS) then
            xi = FP_NAN
            yi = FP_NAN
            if (present(ta)) then
                ta = FP_NAN
            endif
            if (present(tb)) then
                tb = FP_NAN
            endif
            ! Check if point qa is on second line
            ovl = mx * dy - dx * my
            if (abs(ovl) < FP_EPS) then
                info = 1
            else
                info = -1
            endif
            return
        endif

        ! Determine parameters of intersection of lines
        ra = (ny * dx - nx * dy) / det
        rb = (my * dx - mx * dy) / det

        if (present(ta)) then
            ta = ra
        endif
        if (present(tb)) then
            tb = rb
        endif

        xi = pb(1) + nx * rb
        yi = pb(2) + ny * rb

        info = 0

    end subroutine

end module