c
c     simple program to read ascii xyz contour line file
c     (converted from Oasis dxf by NGA conversion program)
c     and find closed contours
c
c     Rick Saltus - 16 Apr 2001 - NPRA project
c
      program findclosures
c
c     xlist,ylist,zlist store first point of each contour line
c     nstatus is 0 if open contour, 1 if closed
c     nstart is starting line of contour in xyz file
c     nend is ending line of contour in xyz file
c     area is area enclosed by closed contour
c
      parameter (nmax=60000)
c     dimension xlist(nmax),ylist(nmax),zlist(nmax),nstatus(nmax)
      dimension nstatus(nmax)
      dimension nstart(nmax),nend(nmax),nline(nmax)
      character*60 buffer
c
c     max distance allowed between points before picking up pen
c
c     distmax=5000
c
c     distance where two points are considered "equal"
c
c     close=1150
c
      open(10,file="input.xyz",form='formatted',status='old')
c
      nclosures=0
      nzeros=0
      ncount=1
      lcount=2
      read(10,'(4x,i10)')line
      read(10,*)x,y,z
      nstart(ncount)=1
      nline(ncount)=line
      xlast=x
      xstart=x
      ylast=y
      ystart=y
      zlast=z
      zstart=z
c
   10 continue
      read(10,'(a60)',end=99)buffer
      lcount=lcount+1
      if (buffer(1:4).eq.'LINE') then
        read(buffer,'(4x,i10)')line
        nend(ncount)=lcount-1
        if ((xstart.eq.xlast).and.(ystart.eq.ylast)) then
          nstatus(ncount)=1
          nclosures=nclosures+1
        else
          nstatus(ncount)=0
        end if
        ncount=ncount+1
        nstart(ncount)=lcount
        nline(ncount)=line
        read(10,*)xstart,ystart,zstart
        lcount=lcount+1
      else
        read(buffer,*)xlast,ylast,zlast
      end if
c
c     s1dist=sqrt((xlist(ncount)-x)**2+(ylist(ncount)-y)**2)
c     dist=sqrt((xlast-x)**2+(ylast-y)**2)
c     if ((dist.gt.distmax).or.(z.ne.zlast)) then
c       s2dist=sqrt((xlist(ncount)-xlast)**2
c    &            +(ylist(ncount)-ylast)**2)
c       if (s2dist.le.close)
c    &    then
c         nstatus(ncount)=1
c         nclosures=nclosures+1
c       else
c         nstatus(ncount)=0
c       end if
c       nend(ncount)=lcount-1
c       if (nstart(ncount).eq.nend(ncount)) then
c         nzeros=nzeros+1
c       end if
c       ncount=ncount+1
c       xlist(ncount)=x
c       ylist(ncount)=y
c       zlist(ncount)=z
c       nstart(ncount)=lcount
c     else if (s1dist.le.close) then
c       print *,s1dist,' = close beg/end'
c       nstatus(ncount)=1
c       nclosures=nclosures+1
c       nend(ncount)=lcount
c       read(10,*,end=99)x,y,z
c       ncount=ncount+1
c       xlist(ncount)=x
c       ylist(ncount)=y
c       zlist(ncount)=z
c       nstart(ncount)=lcount
c     end if
c     xlast=x
c     ylast=y
c     zlast=z
      goto 10
c
   99 continue
      nend(ncount)=lcount
      if ((xstart.eq.xlast).and.(ystart.eq.ylast)) then
        nstatus(ncount)=1
        nclosures=nclosures+1
      else
        nstatus(ncount)=0
      end if
      print *,'Reached end of xyz file - writing report1.xyz'
c     print *,distmax,' = max distance'
      print *,ncount,' contour lines found'
      print *,nclosures,' closed contours found'
c     print *,nzeros,' single point countours found'
      close(10)
c
      open(11,file='report1.xyz',form='formatted',status='unknown')
      do 20 i=1,ncount
        npts=nend(i)-nstart(i)
        write(11,100)i,nline(i),nstart(i),nend(i),
     c        npts,nstatus(i)
  100 format(i4,x,5i6)
   20 continue
      close(11)
c
c     split input file into open & closed contour sets
c
      open(10,file="input.xyz",form='formatted',status='old')
      open(12,file="open.xyz",form='formatted',status='unknown')
      open(13,file="closed.xyz",form='formatted',status='unknown')
      do 30 i=1,ncount
        if (nstatus(i).eq.0) then
          read(10,'(a60)')buffer
          do 25 j=nstart(i)+1,nend(i)
            read(10,*,end=999)x,y,z
c           sdist=sqrt((xlist(i)-x)**2+(ylist(i)-y)**2)
c           write(12,110)x,y,z,nline(i),nstatus(i)
   25     continue
        else 
          read(10,'(a60)')buffer
          do 26 j=nstart(i)+1,nend(i)
            read(10,*,end=999)x,y,z
c           sdist=sqrt((xlist(i)-x)**2+(ylist(i)-y)**2)
            write(13,110)x,y,z,nline(i),nstatus(i)
   26     continue
        end if
   30 continue
  110 format(2f15.3,f12.3,i10,i5)
      close(10)
      close(12)
      close(13)
      print *,'Normal end of program'
      stop
c
  999 continue
      print *,'Abnormal termination - end of file too soon'
      end