c  Program DETOURG
c
c  Contouring of gridded 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,maxcol=2000)
      parameter (MAXIX=1203)
      dimension xyz(1203,3), tetr(2401,3), xpnt(3,3), det(2,3)
      dimension isl(5), dxp(2), dyp(2), xp(4), yp(4), x(100), y(100),
     +          ie(10), ipn(12)
c      dimension row(maxcol,2), flx(100,2)
      dimension row(maxcol,2)
      integer itetr(2401,3), istack(2401), ktetr(50,2), itemp(3,2)
      character ifile*80, label*9, title*56, pgm*8
      dimension post(8)
      character ifile2*80, ifile3*80
        logical    binary, expres
c
      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 ddval/1.0e30/
      data abit1/0.5/
c
      ktmax=0
      print '(a\)',' Input grid file: '
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) title,pgm,nc,nr,nz,xo,dx,yo,dy
      xmin=xo
      ymin=yo
      xmax=xo+(nc-1)*dx
      ymax=yo+(nr-1)*dy
      zmin=1.e+38
      zmax=-1.e+38
      do j=1,nr
      read(10) dlt,(row(i,1),i=1,nc)
      do i=1,nc
      if(row(i,1).gt.ddval) go to 23
      if(row(i,1).gt.zmax) zmax=row(i,1)
      if(row(i,1).lt.zmin) zmin=row(i,1)
   23 continue
      end do
      end do
      rewind 10
      read(10) title
      size=1.0
      print *,'xmin = ',xmin,' xmax = ',xmax
      print *,'ymin = ',ymin,' ymax = ',ymax
      print *,'zmin = ',zmin,' zmax = ',zmax
      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
      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 *
   52 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
c
      ifunc = 0
      ifile2(1:1)=' '
      print '(a\)',' Enter point file (or press ENTER):'
      read(5,100) ifile2
      if(ifile2(1:1).ne.' ') then
        call openf(11,ifile2,istack,maxpst,binary,expres,
     1        nword,maxlog,ntotal)
c  set function switches
	if ( binary  .and.  nword .gt. 3 ) then
	  ifunc = 2
          nchan=nword-2
	endif
	if ( binary  .and.  nword .eq. 3 ) then
	  ifunc = 3
	endif
	if ( expres ) ifunc = 1
      endif
c
        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')
        nl=1
   49 abit1=1.e-10
   55 continue
      xr=(xmax-xmin)
      datax=xr
      yr=(ymax-ymin)
      if(abs(yr).gt.abs(datax)) datax=yr
      call pltset(iplotr,xboard,yboard,isl)
      dxp(1)=0.0
      dxp(2)=xr/datax
      dyp(1)=0.0
      dyp(2)=yr/datax
c      print*,dxp,dyp
      yp(1)=yboard-0.7
      xp(1)=yp(1)*abs(xr/yr)
      if(xp(1).gt.xboard-1.2) then
        xp(1)=xboard-1.2
        yp(1)=xp(1)*abs(yr/xr)
      endif
      if(iplotr.eq.5) then
        print*,'Data area is ',xp(1),' x ',yp(1),' inches.'
      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
      li=len_trim(ifile)
      lt=len_trim(title)
      if(lt.gt.0) then
        ifile = ifile(1:li)//' - '//title
        li=len_trim(ifile)
      endif
      call vchar(.15,.35,ifile,li,3,.1,0.,0.,0.)
c
   44 if(conti.eq.0.) go to 53
      read(10) dlt,(row(i,1),i=1,nc)
      do 200 j=2,nr
        do 1 i=1,3
          itetr(1,i)=i
          tetr(1,i)=xpnt(i,3)
          do 1 jj=1,3
    1   xyz(i,jj)=xpnt(i,jj)
        do 2 i=2,2401
    2   istack(i)=i
        read(10) dlt,(row(i,2),i=1,nc)
        ix=3
        ilast=0
        do 290 i=1,nc-1
          boxmin=min(row(i,1),row(i+1,1),row(i,2),row(i+1,2))
          boxmax=max(row(i,1),row(i+1,1),row(i,2),row(i+1,2))
          cont=contm
          do ic=1,lc
            if(boxmin.le.cont.and.boxmax.ge.cont) go to 299
            cont=cont+conti
          end do
          go to 290
299       continue
          if(ilast.ne.i) then
            if(ix.ge.MAXIX-1) go to 291
            ix=ix+1
            xyz(ix,1)=xo+float(i-1)*dx
            xyz(ix,2)=yo+float(j-2)*dy
            xyz(ix,3)=row(i,1)
            ix=ix+1
            xyz(ix,1)=xo+float(i-1)*dx
            xyz(ix,2)=yo+float(j-1)*dy
            xyz(ix,3)=row(i,2)
          endif
          if(ix.ge.MAXIX-1) go to 291
          ix=ix+1
          xyz(ix,1)=xo+float(i)*dx
          xyz(ix,2)=yo+float(j-2)*dy
          xyz(ix,3)=row(i+1,1)
          ix=ix+1
          xyz(ix,1)=xo+float(i)*dx
          xyz(ix,2)=yo+float(j-1)*dy
          xyz(ix,3)=row(i+1,2)
          ilast=i+1
290     continue
c  normalize the data
291     do ic=4,ix
          xyz(ic,1)=(xyz(ic,1)-xmin)/datax
          xyz(ic,2)=(xyz(ic,2)-ymin)/datax
        end do

        isp=1
        id=2
        do 13 nuc=4,ix
          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
            istack(id)=jt
c  add edges to ktetr but delete if already listed
            do 9 ii=1,3
              l1=itemp(ii,1)
              l2=itemp(ii,2)
              if(km.le.0) go to 8
              kmt=km
              do 7 jj=1,kmt
                if(itetr(jt,l1).ne.ktetr(jj,1)) go to 7
                if(itetr(jt,l2).ne.ktetr(jj,2)) go to 7
                km=km-1
                if(jj.gt.km) go to 9
                do 6 k=jj,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
              ktetr(km,1)=itetr(jt,l1)
              ktetr(km,2)=itetr(jt,l2)
    9       continue
   10     continue
c  form new 3-tuples
          do 12 ii=1,km
            kt = istack(id)
            id = id + 1
c  calculate the circumcircle center and radius
c  squared of points ktetr(i,*) and place in tetr(kt,*)
            do 11 jz = 1,2
              i2 = ktetr(ii,jz)
              det(jz,1) = xyz(i2,1) - xyz(nuc,1)
              det(jz,2) = xyz(i2,2) - xyz(nuc,2)
   11         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
            dd = det(1,1) * det(2,2) - det(1,2) * det(2,1)
            if(dd.eq.0.0) stop 'dd'
            if(kt.lt.1) kt=1
            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
            itetr(kt,1) = ktetr(ii,1)
            itetr(kt,2) = ktetr(ii,2)
   12       itetr(kt,3) = nuc
   13     isp = isp + 2
c         print*,isp
        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            print*,cont,xyz(itetr(jt,1),3),xyz(itetr(jt,2),3)
            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
   14       anum  = cont - xyz(itetr(jt,1),3)
            den = xyz(itetr(jt,3),3) - xyz(itetr(jt,1),3) + abit
c            print*,anum,den
            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
c            print*,anum,den
            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
c            print*,anum,den
            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
   19   continue
        do 200 ic=1,nc
        row(ic,1)=row(ic,2)
  200 continue
c plot lines
      if(nl.ne.0) then
        i=0
   43   read(12,*,end=53) flx,fly
        if(flx.ge.1.e+37) go to 41
        i=i+1
        x(i)=(flx-xmin)/datax
        y(i)=(fly-ymin)/datax
        go to 43
   41   call line(x,y,i,0,1500)
        i=0
        go to 43
      endif
c plot points
   53 continue
      if(ifunc.ne.0) then
        go to (31,30,32) ifunc
   30   read(11,end=31) pgm,(post(i),i=1,nchan)
        x1=(post(1)-xmin)/datax
        y1=(post(2)-ymin)/datax
        call vchar(x1,y1,13,1,1500,.1,0.,0.,0.)
        go to 30
   32   read(11,end=31) x1,y1
        if(x1.gt.ddval) go to 32
        x1=(x1-xmin)/datax
        y1=(y1-ymin)/datax
        call vchar(x1,y1,13,1,1500,.1,0.,0.,0.)
        go to 32
   31   continue
      endif
c plot scale
      if(conti.eq.0.) go to 21
      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
   21 call endpt(ie)
c      print*,abit
      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 )
      dimension  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
        dimension 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
      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

