SUBROUTINE Poly_Intercept (a, b, x, y, n, u, v, m, num, ierr) !----------------------------------------------------------------------- ! INTERSECTION OF A STRAIGHT LINE ! AND POLYGONAL PATH !----------------------------------------------------------------------- ! The polygon is defined by the set of points (xi, yi), i = 1, 2, ..., n. ! The straight line is from (a1,a2) to (b1,b2). ! On exit, the arrays U and V contain the num points at which the line ! crosses the polygon in order, provided that num <= m. ! Error indicator: ! ierr = 0 no error detected ! = 1 if a = b ! = 2 U and V require more storage, i.e. num > m. ! = -i if the ith segment of the polygon is coincident with part of the ! line. ! Based upon routine PFIND from the NSWC Mathematics Library. ! Code converted using TO_F90 by Alan Miller ! Date: 2000-07-04 Time: 12:24:01 IMPLICIT NONE REAL, INTENT(IN) :: a(2) REAL, INTENT(IN) :: b(2) REAL, INTENT(IN) :: x(:) REAL, INTENT(IN) :: y(:) INTEGER, INTENT(IN) :: n REAL, INTENT(OUT) :: u(:) REAL, INTENT(OUT) :: v(:) INTEGER, INTENT(IN) :: m INTEGER, INTENT(OUT) :: num INTEGER, INTENT(OUT) :: ierr ! Local variables INTEGER :: i, ind, nm1 REAL :: d, diff, diff1, eps, h, hi, k, ki, onem, onep, p, q, s, & t, tmax, tmin, tol, tol0 !---------------------- ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE ! SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 . eps = EPSILON(1.0) !---------------------- num = 0 IF (n < 2) GO TO 200 h = b(1) - a(1) k = b(2) - a(2) IF (h == 0.0 .AND. k == 0.0) GO TO 200 ierr = 0 nm1 = n - 1 tol = 4.0*eps tol0 = 2.0*eps onep = 1.0 + tol onem = 0.5 + (0.5 - tol0) ind = 0 DO i = 1, nm1 hi = x(i + 1) - x(i) ki = y(i + 1) - y(i) IF (hi == 0.0 .AND. ki == 0.0) CYCLE ind = 1 ! CHECK IF THE LINE FROM A TO B AND THE I-TH ! LINE IN THE PATH ARE PARALLEL s = hi*k t = h*ki d = s - t IF (ABS(d) <= tol*MAX(ABS(s), ABS(t))) GO TO 40 !----------------------------------------------------------------------- ! THE LINES ARE NOT PARALLEL !----------------------------------------------------------------------- p = x(i) - a(1) q = y(i) - a(2) s = hi*q t = ki*p diff = s - t IF (ABS(diff) <= tol*MAX(ABS(s),ABS(t))) diff = 0.0 s = h*q t = k*p diff1 = s - t IF (ABS(diff1) <= tol*MAX(ABS(s),ABS(t))) diff1 = 0.0 s = diff/d t = diff1/d IF (s < 0.0 .OR. s > onep) CYCLE IF (t < 0.0 .OR. t > onep) CYCLE IF (num > 0 .AND. t == 0.0) CYCLE IF (s > 0.0) GO TO 20 ! POINT A IS ON THE I-TH LINE 10 num = num + 1 IF (num > m) GO TO 210 u(num) = a(1) v(num) = a(2) CYCLE ! POINT B IS ON THE I-TH LINE 20 IF (s < onem) GO TO 30 21 num = num + 1 IF (num > m) GO TO 210 u(num) = b(1) v(num) = b(2) CYCLE ! THE INTERIOR OF THE LINE FROM A TO B ! INTERSECTS WITH THE I-TH LINE 30 num = num + 1 IF (num > m) GO TO 210 u(num) = a(1) + s*h v(num) = a(2) + s*k CYCLE !----------------------------------------------------------------------- ! THE LINES ARE PARALLEL !----------------------------------------------------------------------- 40 IF (ABS(hi) > ABS(ki)) GO TO 50 d = a(2) - y(i) IF (ABS(d) <= tol0*MAX(ABS(a(2)),ABS(y(i)))) d = 0.0 s = d/ki p = x(i) + s*hi IF (ABS(a(1) - p) > tol*MAX(ABS(a(1)),ABS(p))) CYCLE d = b(2) - y(i) IF (ABS(d) <= tol0*MAX(ABS(b(2)),ABS(y(i)))) d = 0.0 t = d/ki GO TO 60 50 d = a(1) - x(i) IF (ABS(d) <= tol0*MAX(ABS(a(1)),ABS(x(i)))) d = 0.0 s = d/hi p = y(i) + s*ki IF (ABS(p - a(2)) > tol*MAX(ABS(p),ABS(a(2)))) CYCLE d = b(1) - x(i) IF (ABS(d) <= tol0*MAX(ABS(b(1)),ABS(x(i)))) d = 0.0 t = d/hi ! THE 2 LINES ARE PORTIONS OF THE SAME ! STRAIGHT INFINITE LINE 60 IF (s > 0.0 .AND. s < onem) GO TO 220 IF (t > 0.0 .AND. t < onem) GO TO 220 tmin = MIN(s,t) tmax = MAX(s,t) IF (tmax <= 0.0) GO TO 70 IF (tmin >= onem) GO TO 80 GO TO 220 70 IF (tmax < 0.0) CYCLE IF (num > 0) CYCLE IF (tmax == s) GO TO 10 GO TO 21 80 IF (tmin > 1.0) CYCLE IF (tmin == s) GO TO 10 GO TO 21 END DO IF (ind == 0) GO TO 200 IF (num < 2) RETURN IF (u(num) == x(1) .AND. v(num) == y(1)) num = num - 1 RETURN ! ERROR RETURN 200 ierr = 1 RETURN 210 ierr = 2 num = num - 1 RETURN 220 ierr = -i RETURN END SUBROUTINE Poly_Intercept