c  Program DETOUR
c
c  Contouring of scattered data using Delaunay tessellation
c
c  modified from  D.F. Watson by J. Phillips
c
c  reference: Watson, D.F., 1982, Acord: Automatic contouring of raw data:
c             Computers and Geosciences, v.8, no. 1, p.97-101.
c
      parameter (maxpst=52,npoint=7000,nktetr=500)
      dimension xyz(npoint+3,3), tetr(2*npoint+1,3), xpnt(3,3), det(2,3)
      dimension isl(5), dxp(2), dyp(2), xp(4), yp(4), x(100), y(100),
c     +          ie(10), ipn(12), post(6), flx(100,2)
     +          ie(10), ipn(12), post(6)
      integer*2 itetr(2*npoint+1,3), istack(2*npoint+1),
     +          ktetr(nktetr,2), itemp(3,2)
      character ifile*80, lab*8, label*9, ifile3*80
        logical    binary, expres
      data itemp,xpnt/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./
      data ipn/100,900,300,1100,200,1000,1400,600,400,1200,500,1300/
      data abit1/0.5/,dval/1.e38/,ddval/1.e30/
      print '(a\)',' Input binary post or xyz file name: '
      read(5,100) ifile
  100 format(a)
        call openf(10,ifile,istack,maxpst,binary,expres,
     1        nword,maxlog,ntotal)
c  set function switches
        ifunc = 0
	if ( binary  .and.  nword .gt. 3 ) then
	  ifunc = 2
          nchan=nword-4
	endif
	if ( binary  .and.  nword .eq. 3 ) then
	  ifunc = 3
	endif
        if ( expres ) ifunc = 1
      if(ifunc.eq.2) then
        print '(a\)',' Enter z-field:'
        read(5,*) iz
      endif
      do 1 i=1,3
      itetr(1,i)=i
      tetr(1,i)=xpnt(i,3)
      do 1 j=1,3
    1 xyz(i,j)=xpnt(i,j)
      do 2 i=2,2*npoint+1
    2 istack(i)=i
      size=1.0
c  read the data
      if(ifunc.eq.2) then
222     read(10) lab,xyz(4,1),xyz(4,2),(post(i),i=1,nchan)
        if(post(iz).lt.ddval) then
          xyz(4,3)=post(iz)
        else
          go to 222
        endif
      else if(ifunc.eq.3) then
223     read(10) xyz(4,1),xyz(4,2),xyz(4,3)
        if(xyz(4,3).ge.ddval) go to 223
      endif
      xmin=xyz(4,1)
      xmax=xyz(4,1)
      ymin=xyz(4,2)
      ymax=xyz(4,2)
      zmin=xyz(4,3)
      zmax=xyz(4,3)
      j=5
      do 3 i=5,npoint+3
      if(ifunc.eq.2) then
        read(10,end=4) lab,xyz(j,1),xyz(j,2),(post(k),k=1,nchan)
        if (post(iz).lt.ddval) then
          xyz(j,3)=post(iz)
        else
          go to 3
        endif
      else if(ifunc.eq.3) then
        read(10,end=4) xyz(j,1),xyz(j,2),xyz(j,3)
        if(xyz(j,3).ge.ddval) go to 3
      endif
      if(xyz(j,1).lt.xmin) xmin=xyz(j,1)
      if(xyz(j,1).gt.xmax) xmax=xyz(j,1)
      if(xyz(j,2).lt.ymin) ymin=xyz(j,2)
      if(xyz(j,2).gt.ymax) ymax=xyz(j,2)
      if(xyz(j,3).lt.zmin) zmin=xyz(j,3)
      if(xyz(j,3).gt.zmax) zmax=xyz(j,3)
      j=j+1
    3 continue
      print *,'WARNING - input file contains more than ',npoint,
     +' points'
    4 n=j-1
c eliminate duplicate points
      do 23 i=4,n-1
      do 23 j=i+1,n
      if(xyz(i,1).eq.xyz(j,1).and.xyz(i,2).eq.xyz(j,2)) xyz(j,1)=dval
   23 continue
      m=5
      do 24 i=5,n
      if(xyz(i,1).eq.dval) go to 24
      xyz(m,1)=xyz(i,1)
      xyz(m,2)=xyz(i,2)
      xyz(m,3)=xyz(i,3)
      m=m+1
   24 continue
      n=m-1
      print *,'xmin = ',xmin,' xmax = ',xmax
      print *,'ymin = ',ymin,' ymax = ',ymax
      print *,'zmin = ',zmin,' zmax = ',zmax
c      print ', 'enter contour interval:'
      conti=(zmax-zmin)/12.
   50 print '(a,1p,g15.6,a\)', ' Enter contour interval [',conti,']:'
c      read(5,*) conti
      call getfloat(conti,ierr)
      if(ierr.ne.0) go to 50
      aplotr=10
      print '(a\)',' Enter iplotr (5=HP,8=CGA,9=EGA,10=VGA) [10]:'
c      read(5,*) iplotr
      call getfloat(aplotr,ierr)
      iplotr=aplotr
      if(conti.eq.0.) go to 52
      if(conti.lt.0.0) then
        do 57 i=1,6
        isav=ipn(i)
        ipn(i)=ipn(13-i)
        ipn(13-i)=isav
   57   continue
        conti=-conti
      endif
c      contm=zmin-amod(zmin,conti)
      contm=int(zmin/conti)*conti
      if(contm.lt.zmin) contm=contm+conti
      print *,'minimum contour level = ',contm
      lc=int((zmax-contm)/conti)+1
      print *,'number of contours = ',lc
      print *
c
   52   nl=0
        print '(a\)',' Enter line file (or press ENTER):'
        read(5,100) ifile3
c         print *,ifile
        if(ifile3(1:1).eq.' ') go to 49
        open(12,file=ifile3,form='formatted',status='old')
c        do 40 i=1,100
c        read(12,*,end=48) flx(i,1),flx(i,2)
   40   read(12,*,end=48) flx,fly
c         print *, flx(i,1), flx(i,2)
        if(flx.ge.1.e+37) go to 40
        nl=nl+1
        if(flx.lt.xmin) xmin=flx
        if(flx.gt.xmax) xmax=flx
        if(fly.lt.ymin) ymin=fly
        if(fly.gt.ymax) ymax=fly
c   40   continue
        go to 40
c   48   nl=i-1
   48   rewind(12)
c        close(12)
c      endif
   49 abit1=1.e-10
      if(conti.eq.0.) go to 55
c   51 print '(a,1p,g15.6,a\)',' Enter abit [',abit1,']:'
c      call getfloat(abit1,ierr)
c      if(ierr.ne.0) go to 51
      print*,'computing...'
c      read(5,*) abit1
c
c
c      lc=lc+1
   55 continue
c
c      print *,'Press ENTER to continue'
c      read(5,100) tri
c      print *,'plot triangles (y or n)?:'
c      read(5,100) tri
c      print *,'plot data (y or n)?:'
c      read(5,100) dat
      xr=xmax-xmin
      datax=xr
      yr=ymax-ymin
      if(yr.gt.datax) datax=yr
c  normalize the data
      do 5 i=4,n
      xyz(i,1)=(xyz(i,1)-xmin)/datax
    5 xyz(i,2)=(xyz(i,2)-ymin)/datax
      if(conti.eq.0.) go to 56
      isp=1
      id=2
      do 13 nuc=4,n
      km=0
c  loop thru the established 3-tuples
      do 10 jt=1,isp
c  test if new data point is within the jt circumcircle
      dsq=tetr(jt,3)-(xyz(nuc,1)-tetr(jt,1))**2
      if(dsq.lt.0.0) go to 10
      dsq=dsq-(xyz(nuc,2)-tetr(jt,2))**2
      if(dsq.lt.0.0) go to 10
c  delete this 2-tuple but save its edges
      id=id-1
c      if(id.lt.1) id=1
      istack(id)=jt
c  add edges to ktetr but delete if already listed
      do 9 i=1,3
      l1=itemp(i,1)
      l2=itemp(i,2)
      if(km.le.0) go to 8
      kmt=km
      do 7 j=1,kmt
      if(itetr(jt,l1).ne.ktetr(j,1)) go to 7
      if(itetr(jt,l2).ne.ktetr(j,2)) go to 7
      km=km-1
      if(j.gt.km) go to 9
      do 6 k=j,km
      k1=k+1
      do 6 l=1,2
    6 ktetr(k,l)=ktetr(k1,l)
      go to 9
    7 continue
    8 km=km+1
      if(km.gt.nktetr) stop 'ktetr array dimension exceeded'
      ktetr(km,1)=itetr(jt,l1)
      ktetr(km,2)=itetr(jt,l2)
    9 continue
   10 continue
c  form new 3-tuples
      kmt=km
      do 12 i=1,kmt
      if(i.gt.km) go to 12
      if(id.gt.2*npoint+1) stop 'istack array dimension exceeded'
      kt = istack(id)
      id = id + 1
c  calculate the circumcircle center and radius
c  squared of points ktetr(i,*) and place in tetr(kt,*)
c      print*,i,ktetr(i,1),ktetr(i,2)
      do 11 jz = 1,2
      i2 = ktetr(i,jz)
c      if(i2.lt.1) go to 11
      det(jz,1) = xyz(i2,1) - xyz(nuc,1)
      det(jz,2) = xyz(i2,2) - xyz(nuc,2)
      det(jz,3) = det(jz,1) * (xyz(i2,1) + xyz(nuc,1))/2.0
     +          + det(jz,2) * (xyz(i2,2) + xyz(nuc,2))/2.0
   11 continue
      dd = det(1,1) * det(2,2) - det(1,2) * det(2,1)
      if(abs(dd).lt.1.e-10) then
        print*,'WARNING - Degenerate case encountered.'
        print*,'Contours may be missing, misconected, or incomplete.'
        pause 'Hit <Enter> to plot.'
        go to 56
      endif
c      if(dd.eq.0.0) go to 13
c      if(kt.lt.1) go to 13
      tetr(kt,1) = (det(1,3) * det(2,2) - det(2,3) * det(1,2))/dd
      tetr(kt,2) = (det(1,1) * det(2,3) - det(2,1) * det(1,3))/dd
      tetr(kt,3) = (xyz(nuc,1) - tetr(kt,1))**2
     +           + (xyz(nuc,2) - tetr(kt,2))**2
  115 itetr(kt,1) = ktetr(i,1)
      itetr(kt,2) = ktetr(i,2)
      itetr(kt,3) = nuc
   12 continue
      isp = isp + 2
   13 continue
   56 call pltset(iplotr,xboard,yboard,isl)
      dxp(1)=0.0
      dxp(2)=xr/datax
      dyp(1)=0.0
      dyp(2)=yr/datax
      yp(1)=yboard-0.7
      xp(1)=yp(1)*xr/yr
      if(xp(1).gt.xboard-1.2) then
        xp(1)=xboard-1.2
        yp(1)=xp(1)*yr/xr
      endif
      xp(2)=0
      xp(3)=0.1
      xp(4)=xboard
      yp(2)=0
      yp(3)=0.6
      yp(4)=yboard
      nopts=4
      call scale(dxp,dyp,xp,yp,nopts,ier)
      if(ier.ne.0) stop 'could not scale the plot'
      call neatl
      call vchar(.15,.35,ifile,80,3,.1,0.,0.,0.)
      if(ifunc.eq.2) then
        ifile='z-field = '//char(iz+48)
        x=xp(3)+xp(1)-1.05
        call vchar(x,.35,ifile,11,3,.1,0.,0.,0.)
      endif
      if(nl.ne.0) then
        i=0
c        j=0
c   43   j=j+1
c        if(j.ge.nl-1) go to 44
c        do 42 i=1,nl-1
   43   read(12,*,end=44) flx,fly
        if(flx.ge.1.e+37) go to 41
        i=i+1
        x(i)=(flx-xmin)/datax
        y(i)=(fly-ymin)/datax
c        j=j+1
c   42   continue
c   41   call line(x,y,i-1,0,0)
        go to 43
   41   call line(x,y,i,0,0)
        i=0
        go to 43
      endif
c
   44 if(conti.eq.0.) go to 58
      do 19 jt = 1,isp
      if(itetr(jt,1) .lt. 4 .or. tetr(jt,3) .gt. 1.) go to 19
c  find contour intersections
      abit = 0.0
      if(xyz(itetr(jt,1),3) .eq. xyz(itetr(jt,2),3) .or.
     +   xyz(itetr(jt,1),3) .eq. xyz(itetr(jt,3),3) .or.
     +   xyz(itetr(jt,2),3) .eq. xyz(itetr(jt,3),3)) abit = abit1
      top = amax1(xyz(itetr(jt,1),3),xyz(itetr(jt,2),3),
     +            xyz(itetr(jt,3),3))
      bot = amin1(xyz(itetr(jt,1),3),xyz(itetr(jt,2),3),
     +            xyz(itetr(jt,3),3))
      do 18 jc =  1,lc
      cont=contm+(jc-1)*conti
      ip=1+jc*12/(lc+1)
      jp=1+jc*8/(lc+1)
      if(cont .gt. top .or.  cont .lt. bot) go to 18
c      cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,2),3)
c     +           - xyz(itetr(jt,1),3) + abit)
            anum = cont - xyz(itetr(jt,1),3)
            den =  xyz(itetr(jt,2),3) - xyz(itetr(jt,1),3) + abit
            call fixabit(anum,den,abit,abit1)
            cz = anum/den
      if(cz .le. 0.0 .or. cz .ge. 1.0) go to 15
      x1 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,2),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
      y1 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,2),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
c   14 cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,3),3)
c     +           - xyz(itetr(jt,1),3) + abit)
   14       anum  = cont - xyz(itetr(jt,1),3)
            den = xyz(itetr(jt,3),3) - xyz(itetr(jt,1),3) + abit
            call fixabit(anum,den,abit,abit1)
            cz = anum/den
      if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 16
      x2 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
      y2 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
      go to 17
c   15 cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,3),3)
c     +           - xyz(itetr(jt,1),3) + abit)
   15       anum = cont - xyz(itetr(jt,1),3)
            den = xyz(itetr(jt,3),3) - xyz(itetr(jt,1),3) + abit
            call fixabit(anum,den,abit,abit1)
            cz = anum/den
      if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 18
      x1 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
      y1 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
c   16 cz = (cont - xyz(itetr(jt,2),3))/(xyz(itetr(jt,3),3)
c     +           - xyz(itetr(jt,2),3) + abit)
   16       anum = cont - xyz(itetr(jt,2),3)
            den = xyz(itetr(jt,3),3) - xyz(itetr(jt,2),3) + abit
            call fixabit(anum,den,abit,abit1)
            cz = anum/den
      if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 18
      x2 = (xyz(itetr(jt,2),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,2),1)) * cz) * size
      y2 = (xyz(itetr(jt,2),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,2),2)) * cz) * size
   17 continue
      x(1)=x1
      x(2)=x2
      y(1)=y1
      y(2)=y2
      if(iplotr.eq.5) call newpenhp(jp)
      call line(x,y,2,0,ipn(ip))
   18 continue
c      if(tri.eq.'n') go to 19
c      x(1)=xyz(itetr(jt,1),1)
c      x(2)=xyz(itetr(jt,2),1)
c      y(1)=xyz(itetr(jt,1),2)
c      y(2)=xyz(itetr(jt,2),2)
c      call line(x,y,2,0,801)
c      x(1)=xyz(itetr(jt,1),1)
c      x(2)=xyz(itetr(jt,3),1)
c      y(1)=xyz(itetr(jt,1),2)
c      y(2)=xyz(itetr(jt,3),2)
c      call line(x,y,2,0,801)
c      x(1)=xyz(itetr(jt,2),1)
c      x(2)=xyz(itetr(jt,3),1)
c      y(1)=xyz(itetr(jt,2),2)
c      y(2)=xyz(itetr(jt,3),2)
c      call line(x,y,2,0,801)
   19 continue
c plot scale
      ipold=0
      y=.3
      do 22 i=1,lc
      ip=1+i*12/(lc+1)
      jp=1+i*8/(lc+1)
      cont=contm+(i-1)*conti
      if(ip.ne.ipold) then
      write(label(1:9),'(g9.3)') cont
      y=y+.4
      if(iplotr.eq.5) call newpenhp(jp)
      call vchar(xp(1)+0.4,y,label,9,ipn(ip)+3,.1,0.,0.,0.)
      endif
      ipold=ip
   22 continue
c      if(dat.eq.'n') go to 21
   58 continue
      do 20 i=4,n
      call vchar(xyz(i,1),xyz(i,2),13,1,1500,.1,0.,0.,0.)
   20 continue
   21 call endpt(ie)
      stop
      end
c*******************************************************************************
      subroutine openf( iswt, ifile, ibuff, maxpst, binary, expres,
     1                 nword, maxlog, ntotal )

c  open data file, derive and checkout file structure.
c
c  input file types:
c  xyz      single record xyz data.
c  post     single record containing an id (8 characters), x, y, and
c           1-n dependant-variable channels. Normally n=6, but this
c           program accepts more channels.
c
c  Note: express file options deleted for PC version 1 series.

	parameter  ( maxbuf=1024 )
      integer*2  ibuff(maxbuf)
	character  ifile*(*)
	logical    binary, expres
	nword  = 0
	maxlog = 0
	ntotal = 0
	binary = .false.
	expres = .false.

c  open binary
	open( unit=iswt, file=ifile, status='old', form='unformatted',
     1     mode='read', iostat=ios )
	if ( ios .ne. 0 ) then
        print *,' error opening input data file ( cannot find )'
	  stop
	endif

c  binary or ascii input ?
	read( iswt, err=10, iostat=ios ) itest
10	rewind( iswt )
	if ( ios .eq. 0 ) then
	  binary = .true. 
	  else
	  close( iswt ) 
	  go to 100
	endif

c  file is binary
c  lenrec returns zero if buffer is too short
c      nword = lenrec( iswt, ibuff, maxpst )
      nword = lenrec( iswt, maxpst )
	if ( nword .eq. 1  .or.  nword .eq. 2 ) go to 200

c  simple xyz file.
	if ( nword .eq. 3 ) then
        print *, ' binary xyz data'
	  go to 999
	endif

c  simple post file.
	if ( nword .ge. 4   .and.   nword .le. maxpst ) then
        nchan  = nword - 2
        nnchan=nchan-2
        print *, ' post data,  found', nnchan, ' data channels'
	  go to 999
	endif


c  open for ascii input
100	  open( unit=iswt, file=ifile, status='old', form='formatted', 
     1       mode='read' )
	go to 999

c  record length undefined.
200   print *, ' input record length =', nword
      print *, ' input is neither xyz or post'

999	return
        end
c******************************************************************************
c        function lenrec( iswt, irec, nwork )
        function lenrec( iswt, nwork )
c  brute force recovery of record length
        integer*4 irec(52)

	istep = 10
	if ( nwork .gt. 100 ) istep = 20
	if ( nwork .gt. 500 ) istep = 50
	rewind iswt

	do max = istep, nwork, istep
	  read( iswt, err=10 ) ( irec(i), i=1, max )
	  rewind iswt
	enddo
	lenrec = -1
	go to 99

10	rewind iswt
	do lenrec = max, max - istep, -1
	  read( iswt, err=20 ) ( irec(i), i=1, lenrec )
	  go to 99
20	  rewind iswt
	enddo

99	rewind iswt
	return
        end
c******************************************************************************
         subroutine newpenhp(ipensel)
c
c       written by cindy cooper, january 4, 1982.
c
c       this subroutine is used to select a pen
c       for the hp7580a.  ipensel is the number of
c       pen to be selected.  if ipltr is not equal
c       to the hp7580a, no action is taken.
c
c      include 'pxz99.cmn'
c
c       make sure plot buffer is empty
c        call hpform
c        nxyrs=0
c
c       command the hp7580 to pick up the desired pen,
c       then prepare again for absolute plotting.
c        if (ipltr .eq. 6) then
                write (20,100) ipensel
100             format(' SP',i1,';')
c        endif
        return
        end
c******************************************************************************
      subroutine getfloat(default,ierr)
      character string*80, fmt*40
      ierr=0
      read(*,'(a)') string
      leng=len_trim(string)
      if(leng.ne.0) then
        lenm=leng-index(string(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        write(fmt,'(a,i3,a,i3,a)') '(g',leng,'.',lenm,')'
        read(string(1:leng),fmt,err=990) default
      endif
      go to 999
  990 ierr=1
  999 return
      end
            subroutine fixabit(anum,den,abit,abit1)
            if(anum.eq.0.0) return
131         test = alog10(abs(anum)) - alog10(abs(den))
            if(test.le.38.0) return
              abit = abit*10.
              abit1 = abit
c              print*,abit
              den = den + abit
              go to 131
            end

