c
c     sortclosures.f - input closed.xyz and create closure heirarchy table
c
c     Saltus - 4 Sep 2001 - NPRA project
c
      program sortclosures
c
      Parameter (MAXCLOSE=40000,MAXPTS=1000)
      Dimension xpts(MAXCLOSE,MAXPTS),ypts(MAXCLOSE,MAXPTS),
     & zval(MAXCLOSE),npts(MAXCLOSE),nval(MAXCLOSE),ndir(MAXCLOSE),
     & nparent(MAXCLOSE),xmax(MAXCLOSE),xmin(MAXCLOSE),
     & ymax(MAXCLOSE),ymin(MAXCLOSE)
c
c     xpts, ypts = x,y coords of closed contours
c     zval = z value of contour line
c     npts = # pts in contour
c     nval = contour ID number
c     ndir = +1 for closed high, -1 for closed low, 0 for unknown
c     nparent = 0 for no parent or nval of parent
c
c     print *,'Enter contour increment:'
c     read *,zinc
c
      open(10,file="closed.xyz",form='formatted',status='old')
c
c     set ndir and nparent arrays to zeros
c
      do 5 i=1,MAXCLOSE
        ndir(i)=0
        nparent(i)=0
    5 continue
c
c     read in closed contours from output of findclosures.f
c
      read(10,*)x,y,z,n,junk
      ntest=n
      ncount=1
      npt=1
      nval(ncount)=ncount
      zval(ncount)=z
      xpts(ncount,npt)=x
      ypts(ncount,npt)=y
      xmax(ncount)=x
      xmin(ncount)=x
      ymax(ncount)=y
      ymin(ncount)=y
   10 continue
      read(10,*,end=999)x,y,z,n,junk
      if (n.eq.ntest) then
        npt=npt+1
c       print *,npt,' point input'
        xpts(ncount,npt)=x
        ypts(ncount,npt)=y
        xmax(ncount)=max(x,xmax(ncount))
        xmin(ncount)=min(x,xmin(ncount))
        ymax(ncount)=max(y,ymax(ncount))
        ymin(ncount)=min(y,ymin(ncount))
      else
        npts(ncount)=npt
        ncount=ncount+1
c       print *,ncount,' contour input'
c       print *,n,' orig contour #'
        npt=1
        ntest=n
        nval(ncount)=ncount
        zval(ncount)=z
        xpts(ncount,npt)=x
        ypts(ncount,npt)=y
        xmax(ncount)=x
        xmin(ncount)=x
        ymax(ncount)=y
        ymin(ncount)=y
      end if     
      goto 10
c
  999 continue
      npts(ncount)=npt
      print *,ncount,' contours input to memory...'
      open(14,file="debug.xyz",form='formatted',
     & status='unknown')
c
c     first pass through contours looking for parents and
c     assigning whether the contour is a closed high or low
c
      do 30 i=1,ncount
c
c     find nested contours
c
        do 20 j=1,ncount
          if (j.eq.i) goto 20
          if ((xmax(i).le.xmax(j)).and.(xmin(i).ge.xmin(j)).and.
     &        (ymax(i).le.ymax(j)).and.(ymin(i).ge.ymin(j))) then
c
c     i probably nested inside j
c     test to confirm
c
c           print *,i,' probably nested in',j
            do 15 k=1,npts(i)
              xtest=xpts(i,k)
              ytest=ypts(i,k)
              intest=inside(xtest,ytest,xpts,ypts,npts(j),j,
     &                      MAXCLOSE,MAXPTS)
c     print *,intest,zval(i),zval(j)
              if (intest.eq.0) then
                if (zval(i).lt.zval(j)) then
c
c     i is lower than j
c
                  nparent(i)=j
                  ndir(i)=-1
                else if (zval(i).gt.zval(j)) then
c
c     i is higher than j
c
                  nparent(i)=j
                  ndir(i)=1
                else
c
c     i is same as j
c
                  nparent(i)=j
                  ndir(i)=-1*ndir(j)
                end if
              goto 20
              end if
   15       continue
          end if
   20   continue
   30 continue
c
c     do 30 i=1,ncount
c
c     for each contour, examine contours with zval>zval
c     if the contour polygon is inside one of these, then
c     it is a closed low with that contour as a parent
c     if not, proceed to the next test
c
c     print *,i,' contour test...'
c     xtest=xpts(i,npts(i))
c     ytest=ypts(i,npts(i))
c     itest2=ifix(npts(i)/2)
c     xtest2=xpts(i,itest2)
c     ytest2=ypts(i,itest2)
c     do 20 j=1,ncount
c       if (j.eq.i) goto 20
c       if (zval(j).gt.zval(i)) then
c         print *,j,' inside test...'
c         intest=inside(xtest,ytest,xpts,ypts,npts(j),j,
c    &                 MAXCLOSE,MAXPTS)
c         if (intest.ne.0) then
c         intest1=intest
c         intest=inside(xtest2,ytest2,xpts,ypts,npts(j),j,
c    &                 MAXCLOSE,MAXPTS)
c         if (intest1.ne.intest) then
c           print *,'2nd test diff from first'
c         end if
c         end if
c         print *,'+test, i,j=',i,j,' intest= ',intest,'xtest,ytest=',
c    &            xtest,ytest
c         if (intest.eq.0) then
c           ndir(i)=-1
c           nparent(i)=j
c           goto 30
c         end if
c       end if
c  20 continue
c     print *,' at 20'
c
c     comes here if not a closed low with parent
c     examine contours with zval<zval
c     if the contour polygon is inside one of these, then
c     it is a closed high with that contour as a parent
c     if not, proceed to the next test
c
c     do 25 j=1,ncount
c       if (j.eq.i) goto 25
c       if (zval(j).lt.zval(i)) then
c         print *,j,' -zinc inside test...'
c         intest=inside(xtest,ytest,xpts,ypts,npts(j),j,
c    &                 MAXCLOSE,MAXPTS)
c         print *,'-test, i,j=',i,j,' intest= ',intest,'xtest,ytest=',
c    &            xtest,ytest
c         if (intest.eq.0) then
c           ndir(i)=+1
c           nparent(i)=j
c           goto 30
c         end if
c       end if
c  25 continue
c     print *,' at 25'
c
c     comes here if not a closed high with parent
c     next examine contours with zval=zval
c     if the contour polygon is inside one of these, then
c     it is a closed contour with a parent of opposite sign
c     if not, it is a contour with no parent (i.e. fundamental closed
c     low or high)
c
c     ztest=zval(i)
c     do 27 j=1,ncount
c       if (j.eq.i) goto 27
c       if (zval(j).eq.ztest) then
c         print *,j,' = zval contour test...'
c         intest=inside(xtest,ytest,xpts,ypts,npts(j),j,
c    &                 MAXCLOSE,MAXPTS)
c         if (intest.eq.0) then
c           nparent(i)=j
c           ndir(i)=-1*ndir(j)
c           goto 30
c         end if
c       end if
c  27 continue
c     print *,' at 27'
c     print *,i,' contour has no parent...'
c  30 continue
c
c     print *,' at 30 - completed initial tests'
c
c     make a 2nd pass to resolve two issues:
c
c     1. Set ndir for contours with parent of opposite sign
c        (have parent assigned, but no direction)
c     2. Set ndir for contours with no parent
c
      do 40 i=1,ncount
        if ((ndir(i).eq.0).and.(nparent(i).ne.0)) then
          ndir(i)=-ndir(nparent(i))
        else if (nparent(i).eq.0) then
          do 35 j=1,ncount
            if (j.eq.i) goto 35
            if (nparent(j).eq.i) then
              ndir(i)=ndir(j)
              goto 35
            end if
   35     continue
        end if
   40 continue
c
c     all done with characterization - write out results
c
      open(11,file="closuretable.xyz",form='formatted',
     &     status='unknown')
      open(12,file="chosen.xyz",form='formatted',
     & status='unknown')
      open(13,file="closedstat.xyz",form='formatted',
     & status='unknown')
      npos=0
      nneg=0
      nnull=0
      nhigh=0
      nlow=0
      nundif=0
      do 50 i=1,ncount
        write(11,*)i,npts(i),ndir(i),nparent(i)
        do 43 j=1,npts(i)
          write(13,*)xpts(i,j),ypts(i,j),zval(i),i,ndir(i),
     &               nparent(i)
   43   continue
        if (ndir(i).gt.0) then
          npos=npos+1
        else if (ndir(i).lt.0) then
          nneg=nneg+1
        else
          nnull=nnull+1
        end if
c
        if ((ndir(i).gt.0).and.
     &     ((nparent(i).eq.0).or.(ndir(nparent(i)).lt.0))) then
          nhigh=nhigh+1
          do 45 j=1,npts(i)
            write(12,*)xpts(i,j),ypts(i,j),zval(i),ndir(i),i,npts(i)
   45     continue
        end if
        if ((ndir(i).lt.0).and.
     &     ((nparent(i).eq.0).or.(ndir(nparent(i)).gt.0))) then
           nlow=nlow+1
          do 46 j=1,npts(i)
            write(12,*)xpts(i,j),ypts(i,j),zval(i),ndir(i),i,npts(i)
   46     continue
        end if
        if (ndir(i).eq.0) then
           nundif=nundif+1
          do 47 j=1,npts(i)
            write(12,*)xpts(i,j),ypts(i,j),zval(i),ndir(i),i,npts(i)
   47     continue
        end if
   50 continue
c
c     all done
c
      print *,npos,' highs, ',nneg,' lows,',nnull,' unknowns.'
      print *,nhigh,' closed highs output to chosen.xyz'
      print *,nlow,' closed lows output to chosen.xyz'
      print *,nundif,' closed single-contour closures',
     &  ' output to chosen.xyz'
c
      close(10)
      close(11)
      close(12)
c
      end
c
C+
C      INSIDE - tests whether point is inside polygon
C
C    Call INSIDE (x,y,xbod,ybod,npt,nbd,mbod,mcorn)
C
C    x,y= point to test
C    xbod,ybod= two-dimensional arrays with corners of body
C    npt= number of points in body
C    nbd= body number (first subscript of xbod, ybod)
C    mbod= 1st dimension of xbod,ybod arrays
C    mcorn= 2nd dimension of xbod,ybod arrays
C
C    RETURNS:  0= x,y inside
C              1= x,y outside
C-
C  *********************************************
      Integer Function inside(x,y,xbod,ybod,npt,nbd,mbod,mcorn)
      Dimension xbod(mbod,mcorn),ybod(mbod,mcorn)
c     Double Precision 
C
C     bail for trivial cases npt=1 or 2
C
      if (npt.le.2) then
        inside=1
        return
      end if
C
C     find xmin,xmax,ymin,ymax for polygon
C
      xmin=xbod(nbd,1)
      xmax=xbod(nbd,1)
      ymin=ybod(nbd,1)
      ymax=ybod(nbd,1)
      do 10 i=1,npt
        if (xbod(nbd,i).lt.xmin) xmin=xbod(nbd,i)
        if (xbod(nbd,i).gt.xmax) xmax=xbod(nbd,i)
        if (ybod(nbd,i).lt.ymin) ymin=ybod(nbd,i)
        if (ybod(nbd,i).gt.ymax) ymax=ybod(nbd,i)
   10 continue
c     print *,'xmin,xmax,ymin,ymax=',xmin,xmax,ymin,ymax
C
C     check if test point completely outside polygon box
C
      if ((x.gt.xmax).or.(x.lt.xmin).or.(y.gt.ymax)
     &               .or.(y.lt.ymin)) then
        inside=1
c       print *,'test point outside polygon bounding box'
        return
      end if
C
C     construct radial line from test point to a point
C      outside the polygon min, max box
C
      x2=xmax+2.*abs(xmax)
      y2=ymax+3.*abs(ymax)
c
c
c     print *,'x2,y2=',x2,y2
C
C     check for intersections of radial line with all polygon sides
C     int is number of intersections
C     if int=odd, then pt is inside
C     if int=even, then pt is outside
C
      int=0
      ipar=0
      Do 50 i=1,npt-1
      xp1=xbod(nbd,i)
      yp1=ybod(nbd,i)
      xp2=xbod(nbd,i+1)
      yp2=ybod(nbd,i+1)
      intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c     print *,'i,int=',i,int
      if (intest.ne.0) then
        x2t=xmin-3.*abs(xmin)
        y2t=ymin-2.*abs(ymin)
        intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
        if (intest.ne.0) then
          print *,'Possible missed intersection'
        end if
      end if
   50 Continue
      xp1=xbod(nbd,npt)
      yp1=ybod(nbd,npt)
      xp2=xbod(nbd,1)
      yp2=ybod(nbd,1)
      intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c     print *,'i,int=',i+1,int
      if (intest.ne.0) then
        x2t=xmin-3.*abs(xmin)
        y2t=ymin-2.*abs(ymin)
        intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
        if (intest.ne.0) then
          print *,'Possible missed intersection'
        end if
      end if
C
C     examine int and report results
C
c     print *,'int=',int
      If (mod(int,2).EQ.0) then
        inside=1
c       int1=int
c
c     if first answer is outside, try again (just in case)
c
c       x2=xmin-2.*abs(xmin)
c       y2=ymin-3.*abs(ymin)
c       int=0
c       ipar=0
c       Do 150 i=1,npt-1
c       xp1=xbod(nbd,i)
c       yp1=ybod(nbd,i)
c       xp2=xbod(nbd,i+1)
c       yp2=ybod(nbd,i+1)
c       intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c       print *,'i,int=',i,int
c       if (intest.ne.0) then
c         x2t=xmin-3.*abs(xmin)
c         y2t=ymin-2.*abs(ymin)
c         intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
c         if (intest.ne.0) then
c           print *,'Possible missed intersection'
c         end if
c       end if
c 150   Continue
c       xp1=xbod(nbd,npt)
c       yp1=ybod(nbd,npt)
c       xp2=xbod(nbd,1)
c       yp2=ybod(nbd,1)
c       intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c       print *,'i,int=',i+1,int
c       if (intest.ne.0) then
c         x2t=xmin-3.*abs(xmin)
c         y2t=ymin-2.*abs(ymin)
c         intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
c         if (intest.ne.0) then
c           print *,'Possible missed intersection'
c         end if
c       end if
c       if (mod(int,2).ne.mod(int1,2)) then
c         print *,'Changed result with 2nd insect2 test...'
c       end if
c       if (mod(int,2).EQ.0) then
c         inside=1
c       else
c         inside=0
c       end if
      else
        inside=0
      end if
      Return
      End
C
C     INSECT2  finds intersection of two lines 
C
C  *********************************************
      Integer function insect2(x1,y1,x2,y2,a1,b1,a2,b2,int,ipar)
C
C     if an intersection is found, then increment int
C     if lines are parallel, increment ipar
C
C     insect2=0 normal return
C            =1 intersection falls on a1,b1
C            =2 intersection falls on a2,b2
C
C     conditions 1 & 2 could mean that line is tangent to
C       polygon (no intersection) or that it crosses at 
C       vertex (intersection).  Further testing is required
C       to solve this condition.
C
C     xsect(x1,y1,x2,y2,yloc)=(yloc-y1)*((x2-x1)/(y2-y1))+x1
C 
C     calculate line parameters
C
C     slopes (s) & intercepts (yint)
C
C     ltype = 0 for finite slope
C           = 1 for vertical line
C
      insect2=0
c     print *,'x1,y1,x2,y2=',x1,y1,x2,y2
c     print *,'a1,b1,a2,b2=',a1,b1,a2,b2
      if (x2.eq.x1) then
        ltype1=1
      else
        ltype1=0
        s1=(y2-y1)/(x2-x1)
        yint1=y1-s1*x1
c       print *,'s1,yint1=',s1,yint1
      end if
      if (a2.eq.a1) then
        ltype2=1
      else
        ltype2=0
        s2=(b2-b1)/(a2-a1)
        yint2=b1-s2*a1
c       print *,'s2,yint2=',s2,yint2
      end if
C
C     Set min, max for line segments
C
      if (x1.lt.x2) then
        xmin=x1
        xmax=x2
      else
        xmin=x2
        xmax=x1
      end if
      if (y1.lt.y2) then
        ymin=y1
        ymax=y2
      else
        ymin=y2
        ymax=y1
      end if
      if (a1.lt.a2) then
        amin=a1
        amax=a2
      else
        amin=a2
        amax=a1
      end if
      if (b1.lt.b2) then
        bmin=b1
        bmax=b2
      else
        bmin=b2
        bmax=b1
      end if
C
C     find intersection and test if it is on both segments
C
      if ((ltype1.eq.0).and.(ltype2.eq.0)) then
C
C     both lines not vertical
C
        if (s1.eq.s2) then
C
C     equal slopes = parallel
C
          ipar=ipar+1
        else
C
C     not parallel
C
          xi=(yint2-yint1)/(s1-s2)
          yi=s1*xi+yint1
          if (s1.eq.0) then
C
C     line 1 horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.bmin).and.(yi.lt.bmax)) then
              int=int+1
            end if
          else if (s2.eq.0) then
C
C     line 2 horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.ymin).and.(yi.lt.ymax)) then
              int=int+1
            end if
          else
C
C     neither line horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.ymin).and.(yi.lt.ymax).and.
     &          (yi.gt.bmin).and.(yi.lt.bmax)) then
              int=int+1
            end if
          end if
        end if       
      else if ((ltype1.eq.1).and.(ltype2.eq.0)) then
C
C     line 1 vertical, line 2 not
C
        xi=x1
        yi=s2*xi+yint2
        if (s2.eq.0) then
C
C     line 2 horizontal
C
          if ((xi.gt.amin).and.(xi.lt.amax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax)) then
            int=int+1
          end if
        else
C
C     line 2 not horizontal
C
          if ((xi.gt.amin).and.(xi.lt.amax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        end if
      else if ((ltype1.eq.0).and.(ltype2.eq.1)) then
C
C     line 2 vertical, line 1 not
C
        xi=x2
        yi=s1*xi+yint1
        if (s1.eq.0) then
C
C     line 1 horizontal
C
          if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        else
C
C     line 1 not horizontal
C
          if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        end if
      else
C
C     both lines vertical = parallel
C
        ipar=ipar+1
      end if
C
      Return
      End