SUBROUTINE locpt (x0, y0, x, y, n, l, m) !----------------------------------------------------------------------- ! GIVEN A POLYGONAL LINE CONNECTING THE VERTICES (X(I),Y(I)) (I = 1,...,N) ! TAKEN IN THIS ORDER. IT IS ASSUMED THAT THE POLYGONAL PATH IS A LOOP, ! WHERE (X(N),Y(N)) = (X(1),Y(1)) OR THERE IS AN ARC FROM (X(N),Y(N)) TO ! (X(1),Y(1)). N.B. The polygon may cross itself any number of times. ! (X0,Y0) IS AN ARBITRARY POINT AND L AND M ARE VARIABLES. ! On output, L AND M ARE ASSIGNED THE FOLLOWING VALUES ... ! L = -1 IF (X0,Y0) IS OUTSIDE THE POLYGONAL PATH ! L = 0 IF (X0,Y0) LIES ON THE POLYGONAL PATH ! L = 1 IF (X0,Y0) IS INSIDE THE POLYGONAL PATH ! M = 0 IF (X0,Y0) IS ON OR OUTSIDE THE PATH. IF (X0,Y0) IS INSIDE THE ! PATH THEN M IS THE WINDING NUMBER OF THE PATH AROUND THE POINT (X0,Y0). ! Fortran 66 version by A.H. Morris ! Converted to ELF90 compatibility by Alan Miller, 15 February 1997 !----------------------- IMPLICIT NONE REAL, INTENT(IN) :: x0, y0, x(:), y(:) INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: l, m ! Local variables INTEGER :: i, n0 REAL :: angle, eps, pi, pi2, sum, theta, theta1, thetai, tol, u, v ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE ! SMALLEST NUMBER SUCH THAT 1.0 + EPS > 1.0 eps = EPSILON(1.0) !----------------------------------------------------------------------- n0 = n IF (x(1) == x(n) .AND. y(1) == y(n)) n0 = n - 1 pi = ATAN2(0.0, -1.0) pi2 = 2.0*pi tol = 4.0*eps*pi l = -1 m = 0 u = x(1) - x0 v = y(1) - y0 IF (u == 0.0 .AND. v == 0.0) GO TO 20 IF (n0 < 2) RETURN theta1 = ATAN2(v, u) sum = 0.0 theta = theta1 DO i = 2, n0 u = x(i) - x0 v = y(i) - y0 IF (u == 0.0 .AND. v == 0.0) GO TO 20 thetai = ATAN2(v, u) angle = ABS(thetai - theta) IF (ABS(angle - pi) < tol) GO TO 20 IF (angle > pi) angle = angle - pi2 IF (theta > thetai) angle = -angle sum = sum + angle theta = thetai END DO angle = ABS(theta1 - theta) IF (ABS(angle - pi) < tol) GO TO 20 IF (angle > pi) angle = angle - pi2 IF (theta > theta1) angle = -angle sum = sum + angle ! SUM = 2*PI*M WHERE M IS THE WINDING NUMBER m = ABS(sum)/pi2 + 0.2 IF (m == 0) RETURN l = 1 IF (sum < 0.0) m = -m RETURN ! (X0, Y0) IS ON THE BOUNDARY OF THE PATH 20 l = 0 RETURN END SUBROUTINE locpt