C NCLFORTSTART
      SUBROUTINE hulla(xp,yp,n,xhull,yhull,nhmax,nh)
      integer n
      REAL xp(n), yp(n)
      integer nhmax, nh
      REAL xhull(nhmax), yhull(nhmax)
C NCLEND

C Driver for the  "convex" algoritm
C References:
C Eddy, W. F., A new convex hull algorithm for planar sets,
C ACM Transactions on mathematical software, 3, #4, 398---403.
      PARAMETER (mxn=5000)
      INTEGER in(mxn), ia(mxn), ib(mxn), ih(mxn), il(mxn)
      REAL xx(2,mxn)
      IF (mxn .LT. n) THEN
          WRITE (6,*) ' error in convex hull mxn, n:', mxn, n
      END IF
      m = n
C copy data:
      DO 10 i = 1, n
          in(i) = i
          xx(1,i) = xp(i)
          xx(2,i) = yp(i)
   10 CONTINUE

C get convex hull
      CALL convex(n,xx,m,in,ia,ib,ih,nh,il)

C get linked list
      ik = il(1)
      DO 20 i = 1, nh
          j = ih(ik)
          xhull(i) = xp(j)
          yhull(i) = yp(j)
          ik = il(ik)
   20 CONTINUE
      RETURN
      END

      SUBROUTINE convex(n,x,m,in,ia,ib,ih,nh,il)
C Reference:
C Eddy, W. F., A new convex hull algorithm for planar sets,
C ACM Transactions on mathematical software, 3, #4, 398---403.
C
C THIS SUBROUTINE DETERMINES WHICH OF THE M POINTS OF ARRAY
C X WHOSE SUBSCRIPTS ARE IN ARRAY IN ARE VERTICES OF THE
C MINIMUM AREA CONVEX POLYGON CONTAINING THE M POINTS. THE
C SUBSCRIPTS OF THE VERTICES ARE PLACED IN ARRAY IH IN THE
C ORDER THEY ARE FOUND. NH IS THE NUMBER OF ELEMENTS IN
C ARRAY IH AND ARRAY IL. ARRAY IL IS A LINKED LIST GIVING
C THE ORDER OF THE ELEMENTS OF ARRAY IH IN A COUNTER
C CLOCKWISE DIRECTION. THIS ALGORITHM CORRESPONDS TO A
C PREORDER TRAVERSAL OF A CERTAIN BINARY TREE. EACH VERTEX
C OF THE BINARY TREE REPRESENTS A SUBSET OF THE M POINTS.
C AT EACH STEP THE SUBSET OF POINTS CORRESPONDING TO THE
C CURRENT VERTEX OF THE TREE IS PARTITIONED BY A LINE
C JOINING TWO VERTICES OF THE CONVEX POLYGON. THE LEFT SON
C VERTEX IN THE BINARY TREE REPRESENTS THE SUBSET OF POINTS
C ABOVE THE PARTITIONING LINE AND THE RIGHT SON VERTEX, THE
C SUBSET BELOW THE LINE. THE LEAVES OF THE TREE REPRESENT
C EITHER NULL SUBSETS OR SUBSETS INSIDE A TRIANGLE WHOSE
C VERTICES COINCIDE WITH VERTICES OF THE CONVEX POLYGON.
C FORMAL PARAMETERS
C INPUT
C N  INTEGER           TOTAL NUMBER OF DATA POINTS
C X  REAL ARRAY (2,N)  (X,Y) CO-ORDINATES OF THE DATA
C M  INTEGER           NUMBER OF POINTS IN THE INPUT SUBSET
C IN INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE POINTS
C                      IN THE INPUT SUBSET
C WORK AREA
C IA INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF LEFT SON
C                      SUBSETS. SEE COMMENTS AFTER DIMENSION
C                      STATEMENTS
C IB INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF RIGHT SON
C                      SUBSETS
C OUTPUT
C IH INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE
C                      VERTICES OF THE CONVEX HULL
C NH INTEGER           NUMBER OF ELEMENTS IN ARRAY IH AND
C                      ARRAY IL. SAME AS NUMBER OF VERTICES
C                      OF THE CONVEX POLYGON
C IL INTEGER ARRAY (M) A LINKED LIST GIVING IN ORDER IN A
C                      COUNTER-CLOCKWISE DIRECTION THE
C                      ELEMENTS OF ARRAY IH
      DIMENSION x(2,n)
      DIMENSION in(m), ia(m), ib(m), ih(m), il(m)
C THE UPPER END OF ARRAY IA IS USED TO STORE TEMPORARILY
C THE SIZES OF THE SUBSETS WHICH CORRESPOND TO RIGHT SON
C VERTICES, WHILE TRAVERSING DOWN THE LEFT SONS WHEN ON THE
C LEFT HALF OF THE TREE, AND TO STORE THE SIZES OF THE LEFT
C SONS WHILE TRAVERSING THE RIGHT SONS(DOWN THE RIGHT HALF)
      LOGICAL maxe, mine
      IF (m .EQ. 1) GO TO 220
      il(1) = 2
      il(2) = 1
      kn = in(1)
      kx = in(2)
      IF (m .EQ. 2) GO TO 210
      mp1 = m + 1
      min = 1
      mx = 1
      kx = in(1)
      maxe = .FALSE.
      mine = .FALSE.
C FIND TWO VERTICES OF THE CONVEX HULL FOR THE INITIAL
C PARTITION
      DO 60 i = 2, m
          j = in(i)
          IF (x(1,j)-x(1,kx)) 30, 10, 20
   10     CONTINUE
          maxe = .TRUE.
          GO TO 30
   20     CONTINUE
          maxe = .FALSE.
          mx = i
          kx = j
   30     CONTINUE
          IF (x(1,j)-x(1,kn)) 50, 40, 60
   40     CONTINUE
          mine = .TRUE.
          GO TO 60
   50     CONTINUE
          mine = .FALSE.
          min = i
          kn = j
   60 CONTINUE
C IF THE MAX AND MIN ARE EQUAL, ALL M POINTS LIE ON A
C VERTICAL LINE
      IF (kx .EQ. kn) GO TO 180
C IF MAXE (OR MINE) HAS THE VALUE TRUE THERE ARE SEVERAL
C MAXIMA (OR MINIMA) WITH EQUAL FIRST COORDINATES
      IF (maxe.OR.mine) GO TO 230
   70 CONTINUE
      ih(1) = kx
      ih(2) = kn
      nh = 3
      inh = 1
      nib = 1
      ma = m
      in(mx) = in(m)
      in(m) = kx
      mm = m - 2
      IF (min .EQ. m) min = mx
      in(min) = in(m-1)
      in(m-1) = kn
C BEGIN BY PARTITIONING THE ROOT OF THE TREE
      CALL split(n,x,mm,in,ih(1),ih(2),0,ia,mb,mxa,ib,ia(ma),mxbb)
C FIRST TRAVERSE THE LEFT HALF OF THE TREE
C START WITH THE LEFT SON
   80 CONTINUE
      nib = nib + ia(ma)
      ma = ma - 1
   90 CONTINUE
      IF (mxa .EQ. 0) GO TO 110
      il(nh) = il(inh)
      il(inh) = nh
      ih(nh) = ia(mxa)
      ia(mxa) = ia(mb)
      mb = mb - 1
      nh = nh + 1
      IF (mb .EQ. 0) GO TO 100
      ilinh = il(inh)
      CALL split(n,x,mb,ia,ih(inh),ih(ilinh),1,ia,mbb,mxa,ib(nib),
     $           ia(ma),mxb)
      mb = mbb
      GO TO 80
C THEN THE RIGHT SON
  100 CONTINUE
      inh = il(inh)
  110 CONTINUE
      inh = il(inh)
      ma = ma + 1
      nib = nib - ia(ma)
      IF (ma .GE. m) GO TO 120
      IF (ia(ma) .EQ. 0) GO TO 110
      ilinh = il(inh)
C ON THE LEFT SIDE OF THE TREE, THE RIGHT SON OF A RIGHT SON
C MUST REPRESENT A SUBSET OF POINTS WHICH IS INSIDE A
C TRIANGLE WITH VERTICES WHICH ARE ALSO VERTICES OF THE
C CONVEX POLYGON AND HENCE THE SUBSET MAY BE NEGLECTED.
      CALL split(n,x,ia(ma),ib(nib),ih(inh),ih(ilinh),2,ia,mb,mxa,
     $           ib(nib),mbb,mxb)
      ia(ma) = mbb
      GO TO 90
C NOW TRAVERSE THE RIGHT HALF OF THE TREE
  120 CONTINUE
      mxb = mxbb
      ma = m
      mb = ia(ma)
      nia = 1
      ia(ma) = 0
C START WITH THE RIGHT SON
  130 CONTINUE
      nia = nia + ia(ma)
      ma = ma - 1
  140 CONTINUE
      IF (mxb .EQ. 0) GO TO 160
      il(nh) = il(inh)
      il(inh) = nh
      ih(nh) = ib(mxb)
      ib(mxb) = ib(mb)
      mb = mb - 1
      nh = nh + 1
      IF (mb .EQ. 0) GO TO 150
      ilinh = il(inh)
      CALL split(n,x,mb,ib(nib),ih(inh),ih(ilinh),-1,ia(nia),ia(ma),mxa,
     $           ib(nib),mbb,mxb)
      mb = mbb
      GO TO 130
C THEN THE LEFT SON
  150 CONTINUE
      inh = il(inh)
  160 CONTINUE
      inh = il(inh)
      ma = ma + 1
      nia = nia - ia(ma)
      IF (ma .EQ. mp1) GO TO 170
      IF (ia(ma) .EQ. 0) GO TO 160
      ilinh = il(inh)
C ON THE RIGHT SIDE OF THE TREE, THE LEFT SON OF A LEFT SON
C MUST REPRESENT A SUBSET OF POINTS WHICH IS INSIDE A
C TRIANGLE WITH VERTICES WHICH ARE ALSO VERTICES OF THE
C CONVEX POLYGON AND HENCE THE SUBSET MAY BE NEGLECTED.
      CALL split(n,x,ia(ma),ia(nia),ih(inh),ih(ilinh),-2,ia(nia),mbb,
     $           mxa,ib(nib),mb,mxb)
      GO TO 140
  170 CONTINUE
      nh = nh - 1
      RETURN
C ALL THE SPECIAL CASES ARE HANDLED DOWN HERE
C IF ALL THE POINTS LIE ON A VERTICAL LINE
  180 CONTINUE
      kx = in(1)
      kn = in(1)
      DO 200 i = 1, m
          j = in(i)
          IF (x(2,j) .LE. x(2,kx)) GO TO 190
          mx = i
          kx = j
  190     CONTINUE
          IF (x(2,j) .GE. x(2,kn)) GO TO 200
          min = i
          kn = j
  200 CONTINUE
      IF (kx .EQ. kn) GO TO 220
C IF THERE ARE ONLY TWO POINTS
  210 CONTINUE
      ih(1) = kx
      ih(2) = kn
      nh = 3
      IF ((x(1,kn).EQ.x(1,kx)).AND. (x(2,kn).EQ.x(2,kx))) nh = 2
      GO TO 170
C IF THERE IS ONLY ONE POINT
  220 CONTINUE
      nh = 2
      ih(1) = in(1)
      il(1) = 1
      GO TO 170
C MULTIPLE EXTREMES ARE HANDLED HERE
C IF THERE ARE SEVERAL POINTS WITH THE (SAME) LARGEST
C FIRST COORDINATE
  230 CONTINUE
      IF (.NOT.maxe) GO TO 250
      DO 240 i = 1, m
          j = in(i)
          IF (x(1,j) .NE. x(1,kx)) GO TO 240
          IF (x(2,j) .LE. x(2,kx)) GO TO 240
          mx = i
          kx = j
  240 CONTINUE
C IF THERE ARE SEVERAL POINTS WITH THE (SAME) SMALLEST
C FIRST COORDINATE
  250 CONTINUE
      IF (.NOT.mine) GO TO 70
      DO 260 i = 1, m
          j = in(i)
          IF (x(1,j) .NE. x(1,kn)) GO TO 260
          IF (x(2,j) .GE. x(2,kn)) GO TO 260
          min = i
          kn = j
  260 CONTINUE
      GO TO 70
      END
      SUBROUTINE split(n,x,m,in,ii,jj,s,iabv,na,maxa,ibel,nb,maxb)
C Reference:
C Eddy, W. F., A new convex hull algorithm for planar sets,
C ACM Transactions on mathematical software, 3, #4, 398---403.
C
C THIS SUBROUTINE TAKES THE M POINTS OF ARRAY X WHOSE
C SUBSCRIPTS ARE IN ARRAY IN AND PARTITIONS THEM BY THE
C LINE JOINING THE TWO POINTS IN ARRAY X WHOSE SUBSCRIPTS
C ARE II AND JJ. THE SUBSCRIPTS OF THE POINTS ABOVE THE
C LINE ARE PUT INTO ARRAY IABV, AND THE SUBSCRIPTS OF THE
C POINTS BELOW ARE PUT INTO ARRAY IBEL. NA AND NB ARE,
C RESPECTIVELY, THE NUMBER OF POINTS ABOVE THE LINE AND THE
C NUMBER BELOW. MAXA AND MAXB ARE THE SUBSCRIPTS FOR ARRAY
C X OF THE POINT FURTHEST ABOVE THE LINE AND THE POINT
C FURTHEST BELOW, RESPECTIVELY. IF EITHER SUBSET IS NULL
C THE CORRESPONDING SUBSCRIPT (MAXA OR MAXB) IS SET TO ZERO
C FORMAL PARAMETERS
C INPUT
C N    INTEGER           TOTAL NUMBER OF DATA POINTS
C X    REAL ARRAY (2,N)  (X,Y) CO-ORDINATES OF THE DATA
C M    INTEGER           NUMBER OF POINTS IN INPUT SUBSET
C IN   INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE
C                        POINTS IN THE INPUT SUBSET
C II   INTEGER           SUBSCRIPT FOR ARRAY X OF ONE POINT
C                        ON THE PARTITIONING LINE
C JJ   INTEGER           SUBSCRIPT FOR ARRAY X OF ANOTHER
C                        POINT ON THE PARTITIONING LINE
C S    INTEGER           SWITCH TO DETERMINE OUTPUT. REFER
C                        TO COMMENTS BELOW
C OUTPUT
C IABV INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE
C                        POINTS ABOVE THE PARTITIONING LINE
C NA   INTEGER           NUMBER OF ELEMENTS IN IABV
C MAXA INTEGER           SUBSCRIPT FOR ARRAY X OF POINT
C                        FURTHEST ABOVE THE LINE. SET TO
C                        ZERO IF NA IS ZERO
C IBEL INTEGER ARRAY (M) SUBSCRIPTS FOR ARRAY X OF THE
C                        POINTS BELOW THE PARTITIONING LINE
C NB   INTEGER           NUMBER OF ELEMENTS IN IBEL
C MAXB INTEGER           SUBSCRIPT FOR ARRAY X OF POINT
C                        FURTHEST BELOW THE LINE. SET TO
C                        ZERO IF NB IS ZERO
      DIMENSION x(2,n)
      DIMENSION in(m), iabv(m), ibel(m)
      INTEGER s
C IF S = 2 DONT SAVE IBEL,NB,MAXB.
C IF S =-2 DONT SAVE IABV,NA,MAXA.
C OTHERWISE SAVE EVERYTHING
C IF S IS POSITIVE THE ARRAY BEING PARTITIONED IS ABOVE
C THE INITIAL PARTITIONING LINE. IF IT IS NEGATIVE, THEN
C THE SET OF POINTS IS BELOW.
      LOGICAL t
      t = .FALSE.
C CHECK TO SEE IF THE LINE IS VERTICAL
      IF (x(1,jj) .NE. x(1,ii)) GO TO 10
      xt = x(1,ii)
      dir = sign(1.,x(2,jj)-x(2,ii))*sign(1.,float(s))
      t = .TRUE.
      GO TO 20
   10 CONTINUE
      a = (x(2,jj)-x(2,ii))/ (x(1,jj)-x(1,ii))
      b = x(2,ii) - a*x(1,ii)
   20 CONTINUE
      up = 0.
      na = 0
      maxa = 0
      down = 0.
      nb = 0
      maxb = 0
      DO 60 i = 1, m
          is = in(i)
          IF (t) GO TO 30
          z = x(2,is) - a*x(1,is) - b
          GO TO 40
   30     CONTINUE
          z = dir* (x(1,is)-xt)
   40     CONTINUE
          IF (z .LE. 0.) GO TO 50
C THE POINT IS ABOVE THE LINE
          IF (s .EQ. -2) GO TO 60
          na = na + 1
          iabv(na) = is
          IF (z .LT. up) GO TO 60
          up = z
          maxa = na
          GO TO 60
   50     CONTINUE
          IF (s .EQ. 2) GO TO 60
          IF (z .GE. 0.) GO TO 60
C THE POINT IS BELOW THE LINE
          nb = nb + 1
          ibel(nb) = is
          IF (z .GT. down) GO TO 60
          down = z
          maxb = nb
   60 CONTINUE
      RETURN
      END
      SUBROUTINE locpt(x0,y0,x,y,n,l,m)
      REAL x(n), y(n)
C History
C
C Morris, A. H., 1990, NSWC library of mathematics subroutines,
C Naval Surface Warfare Center, Dahlgren, Vir. 22448
C
C (PJF)  4-4-1992 calculations of eps simplified
C
C-----------------------------------------------------------------------
C     GIVEN A POLYGONAL LINE CONNECTING THE VERTICES (X(I),Y(I))
C     (I = 1,...,N) TAKEN IN THIS ORDER. IT IS ASSUMED THAT THE
C     POLYGONAL PATH IS A LOOP, WHERE (X(N),Y(N)) = (X(1),Y(1))
C     OR THERE IS AN ARC FROM (X(N),Y(N)) TO (X(1),Y(1)).
C
C     (X0,Y0) IS AN ARBITRARY POINT AND L AND M ARE VARIABLES.
C     L AND M ARE ASSIGNED THE FOLLOWING VALUES ...
C
C        L = -1   IF (X0,Y0) IS OUTSIDE THE POLYGONAL PATH
C        L =  0   IF (X0,Y0) LIES ON THE POLYGONAL PATH
C        L =  1   IF (X0,Y0) IS INSIDE THE POLYGONAL PATH
C
C     M = 0 IF (X0,Y0) IS ON OR OUTSIDE THE PATH. IF (X0,Y0)
C     IS INSIDE THE PATH THEN M IS THE WINDING NUMBER OF THE
C     PATH AROUND THE POINT (X0,Y0).
C
C-----------------------
C
C     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
C            SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0
C
C for most of 32-bit unix workstation:
      DATA base, mant/2., 23/
      eps = base** (1-mant)
      n0 = n
      IF (x(1) .EQ. x(n).AND.y(1) .EQ. y(n)) n0 = n - 1
      pi = atan2(0.0,-1.0)
      pi2 = 2.0*pi
      tol = 4.0*eps*pi
      l = -1
      m = 0
C
      u = x(1) - x0
      v = y(1) - y0
      IF (u .EQ. 0.0.AND.v .EQ. 0.0) GO TO 20
      IF (n0 .LT. 2) RETURN
      theta1 = atan2(v,u)
C
      sum = 0.0
      theta = theta1
      DO 10 i = 2, n0
          u = x(i) - x0
          v = y(i) - y0
          IF (u .EQ. 0.0.AND.v .EQ. 0.0) GO TO 20
          thetai = atan2(v,u)
C
          angle = abs(thetai-theta)
          IF (abs(angle-pi) .LT. tol) GO TO 20
          IF (angle .GT. pi) angle = angle - pi2
          IF (theta .GT. thetai) angle = -angle
          sum = sum + angle
          theta = thetai
   10 CONTINUE
C
      angle = abs(theta1-theta)
      IF (abs(angle-pi) .LT. tol) GO TO 20
      IF (angle .GT. pi) angle = angle - pi2
      IF (theta .GT. theta1) angle = -angle
      sum = sum + angle
C
C     SUM = 2*PI*M WHERE M IS THE WINDING NUMBER
C
      m = abs(sum)/pi2 + 0.2
      IF (m .EQ. 0) RETURN
      l = 1
      IF (sum .LT. 0.0) m = -m
      RETURN
C
C     (X0, Y0) IS ON THE BOUNDARY OF THE PATH
C
   20 CONTINUE
      l = 0
      RETURN
      END

