subroutine catalog(ncat,cat,x,ec,e,cc,ne,mc)
      implicit none
      integer ncat,ne,mc
      double precision x(3,ne)
      integer ec(ne)
      integer e(2,mc)
      integer cc(2,mc)
      integer cat(4,mc)
      integer i
      integer kc, kc1
      integer n, n0, n1, n2, nc, nc0, nc1
      integer mcat
      double precision x0, y0, x1, y1, x2, y2
      double precision xp, slope, slope0
      double precision u, v, u0, v0
c
c     set connection numbers negative until catalog is computed
c
      nc = 0
      do i = 1, ne
         if(ec(i) .gt. nc) nc = ec(i)
         ec(i) = -ec(i)
      enddo
      do i = 1, nc
         cc(1,i) = -cc(1,i)
         cc(2,i) = -cc(2,i)
      enddo
c
c     compute catalog
c
      ncat = 0
      do 20 n0 =1, ne
         nc0 = ec(n0)
 10      if(nc0 .lt. 0) then
            nc0 = -nc0
            ec(n0) = nc0
         else if(nc0 .ne. 0) then
            nc1 = nc0
            i = (n0 - e(1,nc0))/(e(2,nc0) - e(1,nc0))
            nc0 = cc(i+1,nc0)
            if(nc0 .gt. 0) goto 10
            nc0 = -nc0
            cc(i+1,nc1) = nc0
         endif
         if(nc0 .eq. 0) goto 20
         n1 = n0
         i = (n1 - e(1,nc0))/(e(2,nc0) - e(1,nc0))
         mcat = ncat + 1
         cat(1,mcat) = n1
         x1 = x(1,n1)
         y1 = x(2,n1)
         n2 = e(2-i,nc0)
         n = 1
         do while(n2 .ne. n0)
            if(n .ge. 4) then
               cat(1,mcat+1) = cat(1,mcat)
               cat(2,mcat+1) = cat(n,mcat)
               mcat = mcat + 1
               n = 2
            endif
            x0 = x1
            y0 = y1
            n1 = n2
            x1 = x(1,n1)
            y1 = x(2,n1)
            n = n + 1
            cat(n,mcat) = n1
            kc = 0
            nc1 = 0
            nc = ec(n1)
            do while(nc .ne. 0)
               if(nc .gt. 0) goto 30
               nc = -nc
               i = (n1 - e(1,nc))/(e(2,nc) - e(1,nc))
               n2 = e(2-i,nc)
               x2 = x(1,n2)
               y2 = x(2,n2)
               u0 = x0 - x1
               v0 = y0 - y1
               u = x2 - x1
               v = y2 - y1
               xp = v0*u - u0*v
               if(xp .le. 0.) goto 30
               slope = (u0*u + v0*v)/xp
               if(kc .ne. 0 .and. slope .lt. slope0) goto 30
               kc1 = nc1
               kc = nc
               slope0 = slope
               goto 30
 30            i = (n1 - e(1,nc))/(e(2,nc) - e(1,nc))
               nc1 = nc
               nc = cc(i+1,nc)
            enddo
            if(kc .eq. 0) goto 10
            i = (n1 - e(1,kc))/(e(2,kc) - e(1,kc))
            n2 = e(2-i,kc)
            if(kc1 .eq. 0) then
               ec(n1) = kc
            else
               i = (n1 - e(1,kc1))/(e(2,kc1) - e(1,kc1))
               cc(i+1,kc1) = kc
            endif
         enddo
         do while(n .lt. 4)
            cat(n+1,mcat) = cat(n,mcat)
            n = n + 1
         enddo
         ncat = mcat
         goto 10
 20   enddo
      return
      end