SUBROUTINE envelope(x, y, n, vertex, nvert, iwk)
! Find the vertices (in clockwise order) of a polygon enclosing
! the points (x(i), y(i), i=1, ..., n.
! On output, vertex(i), i=1, ..., nvert contains the numbers of the vertices.
! iwk() is an integer work array which must have dimension at least n
! in the calling program.
! There is a limit of 100 vertices imposed by the dimension of array next.
! Programmer: Alan Miller
! Latest revision - 12 September 1987
! Fortran 90 version - 8 August 1996
IMPLICIT NONE
INTEGER :: n, vertex(n), nvert, iwk(n)
REAL :: x(n), y(n)
! Local variables
INTEGER :: next(100), i, i1, i2, j, jp1, jp2, i2save, i3, i2next
REAL :: xmax, xmin, ymax, ymin, dist, dmax, dmin, x1, y1, dx, dy, x2, y2,
dx1, dx2, dmax1, dmax2, dy1, dy2, temp, zero = 0.0
IF (n < 2) RETURN
! Choose the points with smallest & largest x- values as the
! first two vertices of the polygon.
IF (x(1) > x(n)) THEN
vertex(1) = n
vertex(2) = 1
xmin = x(n)
xmax = x(1)
ELSE
vertex(1) = 1
vertex(2) = n
xmin = x(1)
xmax = x(n)
END IF
DO i = 2, n-1
temp = x(i)
IF (temp < xmin) THEN
vertex(1) = i
xmin = temp
ELSE IF (temp > xmax) THEN
vertex(2) = i
xmax = temp
END IF
END DO
! Special case, xmax = xmin.
IF (xmax == xmin) THEN
IF (y(1) > y(n)) THEN
vertex(1) = n
vertex(2) = 1
ymin = y(n)
ymax = y(1)
ELSE
vertex(1) = 1
vertex(2) = n
ymin = y(1)
ymax = y(n)
END IF
DO i = 2, n-1
temp = y(i)
IF (temp < ymin) THEN
vertex(1) = i
ymin = temp
ELSE IF (temp > ymax) THEN
vertex(2) = i
ymax = temp
END IF
END DO
nvert = 2
IF (ymax == ymin) nvert = 1
RETURN
END IF
! Set up two initial lists of points; those points above & those below the
! line joining the first two vertices. next(i) will hold the pointer to the
! point furthest from the line joining vertex(i) to vertex(i+1) on the left
! hand side.
i1 = vertex(1)
i2 = vertex(2)
iwk(i1) = -1
iwk(i2) = -1
dx = xmax - xmin
y1 = y(i1)
dy = y(i2) - y1
dmax = zero
dmin = zero
next(1) = -1
next(2) = -1
DO i = 1, n
IF (i == vertex(1) .OR. i == vertex(2)) CYCLE
dist = (y(i) - y1)*dx - (x(i) - xmin)*dy
IF (dist > zero) THEN
iwk(i1) = i
i1 = i
IF (dist > dmax) THEN
next(1) = i
dmax = dist
END IF
ELSE IF (dist < zero) THEN
iwk(i2) = i
i2 = i
IF (dist < dmin) THEN
next(2) = i
dmin = dist
END IF
END IF
END DO
! Ends of lists are indicated by pointers to -ve positions.
iwk(i1) = -1
iwk(i2) = -1
nvert = 2
j = 1
! Start of main process.
! Introduce new vertex between vertices j & j+1, if one has been found.
! Otherwise increase j. Exit if no more vertices.
40 IF (next(j) < 0) THEN
IF (j == nvert) RETURN
j = j + 1
GO TO 40
END IF
jp1 = j + 1
DO i = nvert, jp1, -1
vertex(i+1) = vertex(i)
next(i+1) = next(i)
END DO
jp2 = jp1 + 1
nvert = nvert + 1
IF (jp2 > nvert) jp2 = 1
i1 = vertex(j)
i2 = next(j)
i3 = vertex(jp2)
vertex(jp1) = i2
! Process the list of points associated with vertex j. New list at vertex j
! consists of those points to the left of the line joining it to the new
! vertex (j+1). Similarly for the list at the new vertex.
! Points on or to the right of these lines are dropped.
x1 = x(i1)
x2 = x(i2)
y1 = y(i1)
y2 = y(i2)
dx1 = x2 - x1
dx2 = x(i3) - x2
dy1 = y2 - y1
dy2 = y(i3) - y2
DMAX1 = zero
dmax2 = zero
next(j) = -1
next(jp1) = -1
i2save = i2
i2next = iwk(i2)
i = iwk(i1)
iwk(i1) = -1
iwk(i2) = -1
60 IF (i /= i2save) THEN
dist = (y(i) - y1)*dx1 - (x(i) - x1)*dy1
IF (dist > zero) THEN
iwk(i1) = i
i1 = i
IF (dist > DMAX1) THEN
next(j) = i
DMAX1 = dist
END IF
ELSE
dist = (y(i) - y2)*dx2 - (x(i) - x2)*dy2
IF (dist > zero) THEN
iwk(i2) = i
i2 = i
IF (dist > dmax2) THEN
next(jp1) = i
dmax2 = dist
END IF
END IF
END IF
i = iwk(i)
ELSE
i = i2next
END IF
! Get next point from old list at vertex j.
IF (i > 0) GO TO 60
! End lists with -ve values.
iwk(i1) = -1
iwk(i2) = -1
GO TO 40
END SUBROUTINE envelope